sHM后如何在临近边界的区域设置高粘度

林之流风

大家好~我用snappyHexMesh剖分出了一套微观尺度岩石内部孔隙的网格，在下一步的流动模拟时，想在每个通道临近边界的区域设置一个高粘度（通道中心是另一个粘度），不知如何实现，希望大神们详细指点，非常感谢！

林之流风

在帖子在openfoam对边界层的处理，为什么仅仅修正边界，而不是跟边界相邻的那一层网格？得知跟边界相邻的那一层网格与patchInternalField 有关，能否在求解器程序中，将边界相邻的那一层网格（或区域）的粘度固定为一个值。但又不知如何写起，恳请大神指点，谢谢！

东岳

你的思路正确，你先研究下或者看看其他大佬可以可以给你写个代码我先处理点其他事情

林之流风

非常感谢东岳老师留意，主要问题有这几点：

想加一个类似于这样的代码，但不知是否正确，没有找到类似的例子

nu_.patchInternalField() = 1000 ；

patchInternalField 指的是跟边界相邻的区域，如何设置这个区域的厚度？比如跟边界相邻的两层网格、三层网格？
以simpleFoam求解器为例，这一行代码应该加到createFields.H里，还是应该加到它所调用的粘度模型如Newtonian.C里？

林之流风

@东岳上一条回复中的代码我更改一下```

nu_.patchInternalField() =  scalar(1000)  ；

想加一个类似于这样的代码，不知是否正确

东岳

是的，这种方法虽然可行，但总感觉你这个不符合物理，具体为什么这么做？

微观尺度岩石内部孔隙的网格，在下一步的流动模拟时，想在每个通道临近边界的区域设置一个高粘度

林之流风

@东岳东岳老师您好，我研究的空间尺度是微纳米尺度，想这么做是因为有文章在微观尺度上、通过分子模拟等指出，流体在靠近壁面的粘度与通道中心是不同的，这个临界厚度也就0.7nm左右，于是我就在每个通道临近边界的0.7nm厚的区域设置一个额外的粘度。
有一些困惑，一是不知道上面我写的代码对不对，应该放在哪个位置；二是如果想仅在一个名为“fixedWalls”的path的临近区域设置一个额外的粘度，该如何编写？
感谢老师

林之流风

@东岳老师，这个方法似乎有一些不妥，我修改了simpleFoam求解器调用的Newtonian粘度模型（暂且新命名为“twonuNewtonian”），添加了代码

nu_.patchInternalField() = scalar(1000);

但是wmake libso时报错如下

error: ‘Foam::volScalarField’ has no member named ‘patchInternalField’

请问老师应该如何修改初次接触编程，网上看了一些关于patchInternalField的内容，仍然一知半解，它的本质指的是什么，不知道能不能这样对其赋值。

twonuNewtonian.C 部分代码如下：

bool Foam::viscosityModels::twonuNewtonian::read
(
    const dictionary& viscosityProperties
)
{
    viscosityModel::read(viscosityProperties);

    viscosityProperties_.lookup("nu") >> nu0_;
    nu_ = nu0_;
    nu_.patchInternalField() = scalar(1000);  //此处一行是所添加的
    return true;
}

东岳

你做的挺好的，

const fvPatchList& patches = nu_.mesh().boundary();
forAll(patches, patchi)
{ 
    scalarField one(1.0, nu_.mesh().boundary()[patchi].size());
    nu_.boundaryFieldRef()[patchi].patchInternalField() = one;
}

东岳

我见过一些相关的帖子：

如果设置高粘度，会导致雷诺数非常小，流动非常慢，粘度过大会导致非对角占优引起发散。具体多大也不好说。更好的方式是那种直接把网格扣除掉，设置固定值也不太好，最好直接操作矩阵，类似p.setReference()

你们需要自己测试下，我只是感觉可能会有问题

林之流风

@东岳非常感谢老师深入讲解！这个方法似乎还有一些问题，我添加了您修改的代码如下所示，成功编译。

bool Foam::viscosityModels::twonuNewtonian::read
(
    const dictionary& viscosityProperties
)
{
    viscosityModel::read(viscosityProperties);

    viscosityProperties_.lookup("nu") >> nu0_;
    nu_ = nu0_;
    const fvPatchList& patches = nu_.mesh().boundary();
    forAll(patches, patchi)
    { 
    scalarField one(1000.0, nu_.mesh().boundary()[patchi].size());
    nu_.boundaryFieldRef()[patchi].patchInternalField() = one;
    }
    return true;
}

不过在计算后，发现与更改粘度模型前的计算结果对比发现没有区别，是不是边界条件覆盖了它的作用？

从代码上来看，是对边界临近的一层网格赋予一个粘度值，能否控制对距离边界的N层内网格赋一个粘度，这样边界条件就不会对它影响？
从代码上来看，是对所有边界临近的一层网格赋予一个粘度值，可是只想对某一特定的边界赋值，我在这个网页 openfoam寻找与某一边界相邻的一层网格得知可以通过以下方法定位到边界的ID，如何加到上述代码里来实现呢？

label patchID = mesh.boundaryMesh().findPatchID("cylinder");

老师您曾在这个帖子回复过我用坐标划分区域设置不同粘度，不过经过sHM后剖分的网格，流动通道很多而且壁面都是不规则的，难以通过用坐标来划分，所以想通过边界附近的网格入手去设置。
问题有点多，恳请老师指点迷津，麻烦老师了

林之流风

@东岳老师您好，我尝试了您给出的这行代码

 scalarField one(1.0, nu_.mesh().boundary()[patchi].size());
 nu_.boundaryFieldRef()[patchi].patchInternalField() = one;

成功编译，但是运行时报错

Attempted assignment to a const reference to an object of type N4Foam5FieldIdEE

    From function void Foam::tmp<T>::operator=(const Foam::tmp<T>&) [with T = Foam::Field<double>]

是不是不能这样赋值给patchInternalField()？应该怎么修改呢?

东岳

好像你放的位置不对，对私有成员函数的更改，不能放到const里面

林之流风

@东岳老师，应该放在什么位置呢，我放在了crateFields.H 里面新定义的量下面

Info<< "Reading field nuk\n" << endl;
volScalarField nuk
(
    IOobject
    (
        "nuk",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    ),
    mesh
);
  
   const fvPatchList& patches = nuk.mesh().boundary();
   forAll(patches, patchi)
   {
        scalarField one(3.0, nuk.mesh().boundary()[patchi].size()); 
        nuk.boundaryFieldRef()[patchi].patchInternalField() = one;
   }

还是这种问题，这是为什么呢？

东岳

这个应该放置在粘度模型里面，比如src里面的laminar.C里面

马乔

dimensionedScalar nu
(
    "nu",
    dimViscosity,
    transportProperties
);

wordList patchTypes(mesh.boundaryMesh().size(),"fixedValue");
volScalarField nu0
(
    IOobject
    (
        "nu0",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    nu,
    patchTypes
);

const scalar nuWanted(1000.0);
labelUList neiCells(nu0.boundaryFieldRef()[0].patch().faceCells());
forAll(neiCells, cellI)
{
    nu0.primitiveFieldRef()[neiCells[cellI]]=nuWanted;
}