OpenFOAM控制方程植入代码讨论

[gooseEast]

大家早上好，我是想通过OpenFOAM研究磁流体的同学，我自己做的求解器（下面详说）收敛性极差，而且感觉算得不太对（本身代码就感觉写得很low）。主要想请教大家看看我将数学方程转化为代码的形式哪里有问题，给我些提高指导意见。

首先需要添加电磁场方程
$\nabla \cdot (\sigma \nabla \phi) = 0$（标量运输方程）
$\nabla \cdot (\nabla \mathbf{A}) = -\mu_0 \mathbf{j}$ （矢量运输方程）
其中$\phi$和$\mathbf{A}$是新添加的场电势和磁矢量势，$\sigma$是电导率，$\mu_0$是一个常数，$\mathbf{j} = - \sigma\nabla \phi$是电流密度。
需要在动量方程添加一个源项洛伦兹力$\mathbf{F}_{Lorentz} = \mathbf{j} \times \mathbf{B}$, 其中$\mathbf{B} = \nabla \times \mathbf{A}$。
需要在能量方程添加源项焦耳热$S=\sigma \mathbf{E}^2$和电子运输焓$\frac{5}{2}\frac{k_B}{e}\mathbf{j}\cdot \nabla T$，其中$\mathbf{E} = -\nabla \phi$，$k_B$是玻尔兹曼常数，$e$是元电荷量（常数），$T$是温度。

我以buoyantSimpleFoam为基础，添加了如下代码来求解磁场方程：

solve(fvm::laplacian(sigma, PotE));
solve
(
    fvm::laplacian(A) == Foam::constant::electromagnetic::mu0*sigma*fvc::grad(PotE)
);

动量方程在buoyantSimpleFoam的UEqn.H上只添加了洛伦兹力源项代码如下：

// Solve the Momentum equation      
MRF.correctBoundaryVelocity(U);      
tmp<fvVectorMatrix> tUEqn
(         
    fvm::div(phi, U)       
  + MRF.DDt(rho, U)       
  + turbulence->divDevRhoReff(U) 
  //添加动量方程源项洛伦兹力      
  - sigma*(-fvc::grad(PotE) ^ fvc::curl(A))      
 ==         
    fvOptions(rho, U)     
);     

fvVectorMatrix& UEqn = tUEqn.ref();      
UEqn.relax();      
fvOptions.constrain(UEqn);    
  
if (simple.momentumPredictor())     
{         
    solve         
    (             
        UEqn          
     ==             
        fvc::reconstruct             
        (                 
            (                   
              - ghf*fvc::snGrad(rho)                   
              - fvc::snGrad(p_rgh)                 
            )*mesh.magSf()             
        )         
    );          
    fvOptions.correct(U);     
}

能量方程就是在原有的方程基础上添加了两个源项代码如下：

{
    volScalarField& he = thermo.he();

    fvScalarMatrix EEqn
    (
        fvm::div(phi, he)
      + (
            he.name() == "e"
          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
        )
      - fvm::laplacian(turbulence->alphaEff(), he)

// 添加能量方程源项焦耳热
      - (fvc::grad(PotE) & fvc::grad(PotE))*sigma

// 添加能量方程源项电子运输焓
      - (
            (2.50*Foam::constant::physicoChemical::k / Foam::constant::electromagnetic::e)
          * sigma
          * (-fvc::grad(PotE) & fvc::grad(he)) 
          / thermo.Cp()
        )

     ==
        rho*(U&g)

      + radiation->Sh(thermo, he)

      + fvOptions(rho, he)
    );

    EEqn.relax();

    fvOptions.constrain(EEqn);

    EEqn.solve();

    fvOptions.correct(he);

    thermo.correct();
    radiation->correct();
}

其中$\nabla T$用$\frac{\nabla he}{thermo.Cp()}$表示了。

谢谢大家的批评指正。

[gooseEast]

另外还想请教大家，在OpenFOAM中可以引入柱坐标系吗，轴向，径向，切向。就像Fluent中计算Axisymmetric Swirl模型，将2D拓展成旋转3D计算域。

[李东岳]

一般收敛性的问题主要在于耦合项。你得动量方程耦合项处理的不是很嘎子。如果是震荡的话，有可能就是我这个链接里面最后一个动图的情况 http://dyfluid.com/class.html 最好是采用这个链接里面说的reconstruct方法来处理 http://dyfluid.com/buoyantBoussinesqPimpleFoam.html

另外你说的收敛性极差是怎么个表现。如果不是震荡可能还是别的原因。

至于球柱坐标，这个openfoam内部自行做变换不需要特殊关注。

[gooseEast]

首先谢谢李老师的回复，我的计算模型长这个样子

最常见的报错信息是这个样子：

--> FOAM FATAL ERROR: (openfoam-2306)
Maximum number of iterations exceeded: 100 when starting from T0:300 old T:18296.6 new T:3283.45 f:1.71386e+07 p:100000 tol:0.03

简言之就是发散了。
上面这个是二维网格，计算结果长这个样子（都给把计算结果对称了）：速度云图

温度云图：

感觉计算结果就是错的，因为在Cathode是wall类型，从Inlet附近到Cathode附近流线应该有弯曲的表现。

阴极我设定的边界条件3500K，它最高温度才3800K，磁场产生的源项作用效果十分小。

谢谢李老师提的reconstruct方法，我研究试试。

[liuaowei]

我最近也是用OpenFOAM算一些关于等离子体的问题，有机会可以讨论一下