Navigation

    CFD中文网

    CFD中文网

    • Login
    • Search
    • 最新

    关于 phi 和 fvVectorMatrix 的两个问题

    OpenFOAM
    6
    23
    14423
    Loading More Posts
    • Oldest to Newest
    • Newest to Oldest
    • Most Votes
    Reply
    • Reply as topic
    Log in to reply
    This topic has been deleted. Only users with topic management privileges can see it.
    • M
      mengweilm425 last edited by 李东岳

      刚开始学习openfoam,遇到两个问题, 前来请教。

      我在cavity的算例里面建立了4*4的网格,如下

      vertices
      (
      (0 0 0)
      (1 0 0)
      (1 1 0)
      (0 1 0)
      (0 0 0.1)
      (1 0 0.1)
      (1 1 0.1)
      (0 1 0.1)
      );
      
      blocks
      (
      hex (0 1 2 3 4 5 6 7) (4 4 1) simpleGrading (1 1 1)
      );
      

      然后输出了flux phi的值

      Info<< " phi " << phi << nl <<endl;
      
      Info<< " phi size " << phi.size() << nl <<endl;
      

      得到如下结果:

      phi dimensions [0 3 -1 0 0 0 0];
      
      internalField nonuniform List<scalar>
      24
      (
      -2.28233e-06
      2.28233e-06
      -3.2402e-06
      9.57866e-07
      -2.28434e-06
      -9.55866e-07
      -2.28433e-06
      -3.94416e-06
      6.22649e-06
      -5.34147e-06
      2.35517e-06
      -3.95016e-06
      -2.34716e-06
      -6.2345e-06
      -6.72257e-06
      1.29491e-05
      -7.99247e-06
      3.62507e-06
      -6.71657e-06
      -3.62306e-06
      -1.29511e-05
      1.29491e-05
      1.65741e-05
      1.29511e-05
      )
      ;
      
      boundaryField
      {
      movingWall
      {
      type calculated;
      value uniform 0;
      }
      fixedWalls
      {
      type calculated;
      value uniform 0;
      }
      frontAndBack
      {
      type empty;
      value nonuniform 0();
      }
      }
      
      
      phi size 24
      
      

      我的第一个问题是, 为什么phi只在internal face上面赋值了, 边界上的face不用计算吗?在 fvm::div(phi, U) 也用到了phi值, 边界上的flux是怎么考虑的?

      第二个问题是, 我声明了一个新的fvVectorMatrix

      fvVectorMatrix UEqn_time_term
      (
      fvm::ddt(U)
      );
      

      查看源项的值

      Info<< " UEqn_time_term.s " << UEqn_time_term.source() << nl << endl;
      

      输出结果:

      UEqn_time_term.s
      16
      (
      (-5.83396e-06 5.74984e-06 0)
      (-1.42259e-05 2.43986e-06 0)
      (-1.42171e-05 -2.34961e-06 0)
      (-5.92331e-06 -5.6699e-06 0)
      (-1.00527e-05 2.19123e-05 0)
      (-2.40783e-05 8.37173e-06 0)
      (-2.41399e-05 -8.28764e-06 0)
      (-1.01346e-05 -2.19111e-05 0)
      (-1.39378e-05 4.90818e-05 0)
      (-3.80295e-05 1.57254e-05 0)
      (-3.80131e-05 -1.56865e-05 0)
      (-1.39231e-05 -4.90753e-05 0)
      (9.96865e-05 3.28844e-05 0)
      (7.57897e-05 9.67425e-06 0)
      (7.57862e-05 -9.65121e-06 0)
      (9.97123e-05 -3.2868e-05 0)
      )
      

      这个fvmatrix的源项不是应该为零吗?为什么会有值呢?

      非常感谢!!

      1 Reply Last reply Reply Quote
      • 李东岳
        李东岳 管理员 last edited by 李东岳

        Hi,

        我的第一个问题是, 为什么phi只在internal face上面赋值了, 边界上的face不用计算吗?在 fvm::div(phi, U) 也用到了phi值, 边界上的flux是怎么考虑的?

        phi是通量值,只在面上有值,也就是说只有面有通量。看这个:http://www.dyfluid.com/flux.html 边界处的通量计算而来。要么给定,要么通过边界条件计算。

        这个fvmatrix的源项不是应该为零吗?为什么会有值呢?

        不能是0. 是零的话:AU=0。你的速度都是零。具体的在离散过程中:
        0_1474522774670_捕获.JPG
        对时间的离散调用了显性的值。

        线上CFD课程 7月1日报名截止 http://dyfluid.com/class.html
        CFD高性能服务器 http://dyfluid.com/servers.html

        M 1 Reply Last reply Reply Quote
        • M
          mengweilm425 last edited by CFD中文网

          东岳老师,非常感谢您的回复!
          对于第一问题,我又进行了尝试,下面是我的理解,如果有不对的地方,还请您指正。
          边界上的phi可以由边界条件得到,在cavity这个算例中,边界上没有法向的速度,因而phi全部为‘0’。而在internal face上的phi,可以由已知的速度场进行差值得到face上的速度,再乘以该face的面积和法向,就得到了phi。

          对于第二个问题,我又输出了U的值。

          Info<< "   U   " << U << nl << endl;
          

          得到

             U   dimensions      [0 1 -1 0 0 0 0];
          
          internalField   nonuniform List<vector>
          16
          (
          (-0.00109257 0.00108567 0)
          (-0.00260853 0.000420926 0)
          (-0.00260827 -0.000416558 0)
          (-0.0010966 -0.00108186 0)
          (-0.00181643 0.00400068 0)
          (-0.00421015 0.00141225 0)
          (-0.00421276 -0.00140849 0)
          (-0.00182014 -0.00400052 0)
          (-0.0033937 0.00937929 0)
          (-0.00774594 0.0025887 0)
          (-0.00774527 -0.00258687 0)
          (-0.00339316 -0.00937899 0)
          (0.0210076 0.00646146 0)
          (0.0145416 0.00159495 0)
          (0.0145415 -0.00159377 0)
          (0.0210087 -0.00646063 0)
          )
          ;
          
          boundaryField
          {
              movingWall
              {
                  type            fixedValue;
                  value           uniform (1 0 0);
              }
              fixedWalls
              {
                  type            noSlip;
              }
              frontAndBack
              {
                  type            empty;
              }
          }
          

          发现 UEqn_time_term.source() 里面的值对应为速度的值除以时间步长再乘以对应cell的体积。(时间步长为0.001,体积为0.01乘以0.025乘以0.025)

          1 Reply Last reply Reply Quote
          • C
            CFD中文网 last edited by

            0_1474590590150_捕获.JPG
            这一项吧?对应icoFoam解析公式6的第一项。
            这个研究思路很好。小伙子我看好你。:sunglasses: 哈哈

            边界上的phi可以由边界条件得到,在cavity这个算例中,边界上没有法向的速度,因而phi全部为‘0’。而在internal face上的phi,可以由已知的速度场进行差值得到face上的速度,再乘以该face的面积和法向,就得到了phi。

            是的

            CFD中国标准用户测试帐号
            目前由徐笑笑登录

            M 1 Reply Last reply Reply Quote
            • M
              mengweilm425 @李东岳 last edited by mengweilm425

              @李东岳

              东岳老师,我在研究完icoFoam中的 ‘fvm::ddt(U)’ 和 ‘fvm::div(phi, U)’ 之后又进而研究了‘fvm::laplacian(nu, U)’,有如下几个问题没有弄明白。

              为了看的清楚,附上网格编号的图。

              0_1474592981661_grid.JPG

              我分别定义了不同的 fvVectorMatrix

                      fvVectorMatrix UEqn_time_term
                      (
                         fvm::ddt(U)
                      );
              
                      fvVectorMatrix UEqn_advection_term
                      (
                         fvm::div(phi, U)
                      );
              
                      fvVectorMatrix UEqn_diffusion
                      (
                         fvm::laplacian(nu, U)
                      );
              

              对他们进行输出之后发现 UEqn_time_term 和 UEqn_advection_term 的diag() 和A() 每个元素都是相差一个cell 的体积

              UEqn.A() = UEqn.diag()/mesh.V()
              

              但是UEqn_diffusion的diag() 和A()只有在非边界的单元才满足上述关系。如下

                UEqn_diffusion.diag
              16
              (
              -0.0002
              -0.0003
              -0.0003
              -0.0002
              -0.0003
              -0.0004 \\这个满足上面的关系
              -0.0004 \\这个满足上面的关系
              -0.0003
              -0.0003
              -0.0004 \\这个满足上面的关系
              -0.0004 \\这个满足上面的关系
              -0.0003 
              -0.0002
              -0.0003
              -0.0003
              -0.0002
              )
              
              
                UEqn_diffusion.A  dimensions      [0 0 -1 0 0 0 0];
              
              internalField   nonuniform List<scalar>
              16
              (
              -96
              -80
              -80
              -96
              -80
              -64
              -64
              -80
              -80
              -64
              -64
              -80
              -96
              -80
              -80
              -96
              )
              ;
              
              boundaryField
              {
                  movingWall
                  {
                      type            extrapolatedCalculated;
                      value           nonuniform List<scalar> 4(-96 -80 -80 -96);
                  }
                  fixedWalls
                  {
                      type            extrapolatedCalculated;
                      value           nonuniform List<scalar>
              12
              (
              -96
              -80
              -80
              -96
              -96
              -80
              -80
              -96
              -96
              -80
              -80
              -96
              )
              ;
                  }
                  frontAndBack
                  {
                      type            empty;
                  }
              }
              

              我的第一个问题是 在UEqn_diffusion中.A()是怎么计算的?

              我在输出UEqn_diffusion的时候也发现它与其他的输出结果不同,在网上搜了半天也没有得到答案。

              Info<< "  UEqn_diffusion  " << UEqn_diffusion << nl << endl;
              

              输出结果如下:

                UEqn_diffusion  false true true
              
              16 \\这个应该是diag()
              (
              -0.0002
              -0.0003
              -0.0003
              -0.0002
              -0.0003
              -0.0004
              -0.0004
              -0.0003
              -0.0003
              -0.0004
              -0.0004
              -0.0003
              -0.0002
              -0.0003
              -0.0003
              -0.0002
              )
              
              24 \\这个应该是upper
              (
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              0.0001
              )
              
              [0 4 -2 0 0 0 0]
              16{(0 0 0)} \\我猜这个是source,这里面没有源项,4by4的网格,应该都是零。如果不对,请您指正。
              
              3 \\从这个开始往下我不知道是什么了,其他的matrix这些项都是零,请问这些值都代表了什么?
              (
              4((-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002))
              
              12
              (
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              (-0.0002 -0.0002 -0.0002)
              )
              
              0()
              )
              
              
              3
              (
              4((-0.0002 -0 -0) (-0.0002 -0 -0) (-0.0002 -0 -0) (-0.0002 -0 -0))
              12{(-0 -0 -0)}
              0()
              )
              

              第二个问题 写在了代码的注释里,主要是不明白这个输出的结果里面除了diag(), upper() 和source()之外的那些项是什么?

              我在推倒UEqn_diffusion 的系数矩阵的时候发现,由边界条件得到边界上的速度值,那么最上面的滑动边界上的cell (编号为P)

              \begin{equation}
              S_{e}\nu_{e}\frac{u_{E}-u_{P}}{\Delta x} + S_{w}\nu_{w}\frac{u_{W}-u_{P}}{\Delta x} + S_{n}\nu_{n}\frac{u_{N}-u_{P}}{\Delta y} + S_{s}\nu_{s}\frac{u_{S}-u_{P}}{\Delta y}
              \end{equation}

              里面的$u_{N}$可以根据边界条件差值

              \begin{equation}
              u_{N} = 2 u_{n} - u_{P}
              \end{equation}

              其中$u_{n}$为上面滑动边界已知的速度。

              但是这个边界条件并没有在系数矩阵的diag()中体现出来,这也带来了我的第三个问题,这个边界条件在 UEqn_diffusion 是怎么体现出来的?是否与第二个问题有关?

              第四个问题就是也是又第三个问题引出来的,我试着求解了一下UEqn_diffusion,并输出了一下U的结果

              solve(UEqn_diffusion);
              
              Info<< "   U solution  " << U << nl << endl;
              

              但这次与前面的 UEqn_time_term 和 UEqn_advection_term 不同, 没能找到其结果和UEqn.A() 与UEqn.H() 的关系。

              我的问题有点儿多,我也会继续查找相关的信息,期待您的解答!

              非常感谢!

              M 1 Reply Last reply Reply Quote
              • M
                mengweilm425 @CFD中文网 last edited by

                @cfd-china 谢谢您!我是刚开始学Openfoam,希望能多多向前辈们请教!

                1 Reply Last reply Reply Quote
                • M
                  mengweilm425 @mengweilm425 last edited by CFD中文网

                  @mengweilm425

                  补充一下上面的帖子

                  solve(UEqn_diffusion);
                  
                  Info<< "   U solution  " << U << nl << endl;
                  

                  结果如下:

                  smoothSolver:  Solving for Ux, Initial residual = 0.921181, Final residual = 9.32208e-06, No Iterations 11
                  smoothSolver:  Solving for Uy, Initial residual = 0.999778, Final residual = 8.66271e-06, No Iterations 8
                     U solution  dimensions      [0 1 -1 0 0 0 0];
                  
                  internalField   nonuniform List<vector>
                  16
                  (
                  (0.0189067 -1.2469e-08 0)
                  (0.0430646 -3.74215e-08 0)
                  (0.0430639 -4.90001e-08 0)
                  (0.0189057 -2.83642e-08 0)
                  (0.0703756 -3.73926e-08 0)
                  (0.153354 -1.08185e-07 0)
                  (0.153353 -1.33266e-07 0)
                  (0.070374 -6.9139e-08 0)
                  (0.179619 -4.88555e-08 0)
                  (0.34663 -1.32978e-07 0)
                  (0.34663 -1.48122e-07 0)
                  (0.179618 -6.64227e-08 0)
                  (0.481091 -2.81724e-08 0)
                  (0.706929 -6.88449e-08 0)
                  (0.706929 -6.63391e-08 0)
                  (0.481091 -2.54916e-08 0)
                  )
                  ;
                  
                  boundaryField
                  {
                      movingWall
                      {
                          type            fixedValue;
                          value           uniform (1 0 0);
                      }
                      fixedWalls
                      {
                          type            noSlip;
                      }
                      frontAndBack
                      {
                          type            empty;
                      }
                  }
                  

                  但这次与前面的 UEqn_time_term 和 UEqn_advection_term 不同, 没能找到其结果和UEqn.A() 与UEqn.H() 的关系。感觉openfoam求解了这个方程,从UEqn_diffusion的系数和source()项为零来看,求解结果应该是零,但是显然是我想错了。如果是$u_{P}$为未知,$u_{W}$等相邻的cell的速度用的是上一个时间步的值,那么线性方程组只存在对角线的值,不需要迭代,直接就能出结果。在这个问题上没有想明白,还希望您能指点一下这个方程是怎么求解的。

                  1 Reply Last reply Reply Quote
                  • 李东岳
                    李东岳 管理员 last edited by

                    不好意思信息量太大。

                    如果对时间离散fvm::ddt(T),其离散后的为对角阵。因此如果源项为0,AT=0则T也为零。并且不需要迭代求解。对于对流项fvm::div(phi,T),离散后非对角阵因此不一定T为0,且需要迭代求解。

                    线上CFD课程 7月1日报名截止 http://dyfluid.com/class.html
                    CFD高性能服务器 http://dyfluid.com/servers.html

                    W 1 Reply Last reply Reply Quote
                    • W
                      wwzhao 教授 @李东岳 last edited by

                      @李东岳 fvMatrix的source并不是线性代数方程组Ax=b中的b,而是指加在网格单元中心的力源项。

                      李东岳 1 Reply Last reply Reply Quote
                      • 李东岳
                        李东岳 管理员 @wwzhao last edited by

                        @wwzhao
                        即使没有在方程等号右边添加源项的fvMatrix,比如:

                            fvScalarMatrix TEqn
                            (
                                fvm::ddt(T) 
                            );
                        

                        这个矩阵系统依然具有源项。显性离散的值全部进入了fvMatrix的source()。

                        P.S. 这个问题@mengweilm425 不是在一楼测试过了么?
                        用一个网格单元测试:

                        #include "fvCFD.H"
                        
                        // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
                        
                        int main(int argc, char *argv[])
                        {
                            #include "setRootCase.H"
                        
                            #include "createTime.H"
                            #include "createMesh.H"
                        
                            volScalarField T
                            (
                                IOobject
                                (
                                    "T",
                                    runTime.timeName(),
                                    mesh,
                                    IOobject::NO_READ,
                                    IOobject::NO_WRITE
                                ),
                                mesh,
                                1
                            );
                        
                            fvScalarMatrix TEqn
                            (
                                fvm::ddt(T) 
                            );
                        
                            Info<< "source" << TEqn.source() << endl;
                        }
                        

                        线上CFD课程 7月1日报名截止 http://dyfluid.com/class.html
                        CFD高性能服务器 http://dyfluid.com/servers.html

                        W 1 Reply Last reply Reply Quote
                        • 李东岳
                          李东岳 管理员 last edited by 李东岳

                          针对@mengweilm425 提到的A()和diag()。区别是A()=diag()/V;可测试:

                          #include "fvCFD.H"
                          
                          // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
                          
                          int main(int argc, char *argv[])
                          {
                              #include "setRootCase.H"
                          
                              #include "createTime.H"
                              #include "createMesh.H"
                          
                              volScalarField T
                              (
                                  IOobject
                                  (
                                      "T",
                                      runTime.timeName(),
                                      mesh,
                                      IOobject::NO_READ,
                                      IOobject::NO_WRITE
                                  ),
                                  mesh,
                                  1
                              );
                          
                              fvScalarMatrix TEqn
                              (
                                  fvm::ddt(T) 
                              );
                          
                              Info<< "source" << TEqn.source() << endl;
                              Info<< "source" << TEqn.A() << endl;
                              Info<< "source" << TEqn.diag() << endl;//diag() equals TEqn.A() multiply mesh.V()
                          }
                          

                          源代码:

                          template<class Type>
                          Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
                          {
                              tmp<volScalarField> tAphi
                              (
                                  new volScalarField
                                  (
                                      IOobject
                                      (
                                          "A("+psi_.name()+')',
                                          psi_.instance(),
                                          psi_.mesh(),
                                          IOobject::NO_READ,
                                          IOobject::NO_WRITE
                                      ),
                                      psi_.mesh(),
                                      dimensions_/psi_.dimensions()/dimVol,
                                      extrapolatedCalculatedFvPatchScalarField::typeName
                                  )
                              );
                          
                              tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V();
                              tAphi.ref().correctBoundaryConditions();
                          
                              return tAphi;
                          }
                          
                          template<class Type>
                          Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
                          {
                              tmp<scalarField> tdiag(new scalarField(diag()));
                              addCmptAvBoundaryDiag(tdiag.ref());
                              return tdiag;
                          }
                          

                          线上CFD课程 7月1日报名截止 http://dyfluid.com/class.html
                          CFD高性能服务器 http://dyfluid.com/servers.html

                          W 1 Reply Last reply Reply Quote
                          • W
                            wwzhao 教授 @李东岳 last edited by

                            @李东岳 我是这么理解的,fvm::ddt 反映的是上一时刻(或迭代步)对当前时刻(或迭代步)的物理量的影响,OpenFOAM把这种影响当作源项处理,因此将离散后的系数放在source()里面。

                            1 Reply Last reply Reply Quote
                            • W
                              wwzhao 教授 @李东岳 last edited by 李东岳

                              lduMatrix::diag()和fvMatrix<Type>::D()并不等价。相比lduMatrix::diag(),fvMatrix<Type>::D()考虑了边界对系数矩阵的影响,具体代码对应这句:addCmptAvBoundaryDiag(tdiag.ref());。

                              1 Reply Last reply Reply Quote
                              • C
                                CFD中文网 last edited by CFD中文网

                                有关diag()和D()严格来说是的。你可以测试一下然后给楼主提供更详细的信息。:sunglasses:

                                CFD中国标准用户测试帐号
                                目前由徐笑笑登录

                                W 1 Reply Last reply Reply Quote
                                • W
                                  wwzhao 教授 @CFD中文网 last edited by wwzhao

                                  @cfd-china 离散的具体实现可参阅相关源码。

                                  以fvm::ddt(T)为例,调用的是

                                  template<class Type>
                                  tmp<fvMatrix<Type> >
                                  ddt
                                  (
                                      const GeometricField<Type, fvPatchField, volMesh>& vf
                                  )
                                  {
                                      return fv::ddtScheme<Type>::New
                                      (
                                          vf.mesh(),
                                          vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
                                      )().fvmDdt(vf);
                                  }
                                  

                                  这里调用了fv::ddtScheme<Type>::New,这个函数是一个runtime construction的实现,即在运行时读取字典文件中的离散格式并构造相应的对象:

                                  template<class Type>
                                  tmp<ddtScheme<Type> > ddtScheme<Type>::New
                                  (
                                      const fvMesh& mesh,
                                      Istream& schemeData
                                  )
                                  {
                                      if (fv::debug)
                                      {
                                          Info<< "ddtScheme<Type>::New(const fvMesh&, Istream&) : "
                                                 "constructing ddtScheme<Type>"
                                              << endl;
                                      }
                                  
                                      if (schemeData.eof())
                                      {
                                          FatalIOErrorIn
                                          (
                                              "ddtScheme<Type>::New(const fvMesh&, Istream&)",
                                              schemeData
                                          )   << "Ddt scheme not specified" << endl << endl
                                              << "Valid ddt schemes are :" << endl
                                              << IstreamConstructorTablePtr_->sortedToc()
                                              << exit(FatalIOError);
                                      }
                                  
                                      const word schemeName(schemeData);
                                  
                                      typename IstreamConstructorTable::iterator cstrIter =
                                          IstreamConstructorTablePtr_->find(schemeName);
                                  
                                      if (cstrIter == IstreamConstructorTablePtr_->end())
                                      {
                                          FatalIOErrorIn
                                          (
                                              "ddtScheme<Type>::New(const fvMesh&, Istream&)",
                                              schemeData
                                          )   << "Unknown ddt scheme " << schemeName << nl << nl
                                              << "Valid ddt schemes are :" << endl
                                              << IstreamConstructorTablePtr_->sortedToc()
                                              << exit(FatalIOError);
                                      }
                                  
                                      return cstrIter()(mesh, schemeData);
                                  }
                                  

                                  以常见的Euler格式为例,将会构造EulerDdtScheme类型的对象,并用其fvmDdt()方法返回的fvMatrix,看fvmDdt()的实现:

                                  template<class Type>
                                  tmp<fvMatrix<Type> >
                                  EulerDdtScheme<Type>::fvmDdt
                                  (
                                      const GeometricField<Type, fvPatchField, volMesh>& vf
                                  )
                                  {
                                      tmp<fvMatrix<Type> > tfvm
                                      (
                                          new fvMatrix<Type>
                                          (
                                              vf,
                                              vf.dimensions()*dimVol/dimTime
                                          )
                                      );
                                  
                                      fvMatrix<Type>& fvm = tfvm();
                                  
                                      scalar rDeltaT = 1.0/mesh().time().deltaTValue();
                                  
                                      fvm.diag() = rDeltaT*mesh().Vsc();
                                  
                                      if (mesh().moving())
                                      {
                                          fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc0();
                                      }
                                      else
                                      {
                                          fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc();
                                      }
                                  
                                      return tfvm;
                                  }
                                  

                                  可以看到fvm::ddt(T)产生的fvm.diag()为时间倒数乘以网格单元体积,其产生的源项fvm.source()为时间倒数乘以上一时刻的T乘以网格单元体积。

                                  对于fvm::div及fvm::laplacian可以做类似分析。

                                  插一句,fvm::div离散不会产生source(),fvm::laplacian会产生source()。

                                  M 1 Reply Last reply Reply Quote
                                  • M
                                    mengweilm425 @wwzhao last edited by

                                    @wwzhao @李东岳 @cfd-china

                                    非常感谢大家的回复!我现在的阶段还看不好代码,但是我又对fvMatrix做了一些测试。

                                    还是如下的fvMatrix

                                        fvVectorMatrix UEqn_diffusion
                                        (
                                           fvm::laplacian(nu, U)
                                        );
                                    

                                    我发现在解这个方程的时候

                                    \begin{equation}
                                    [A][x]=[b]
                                    \end{equation}

                                    矩阵[A] 实际上是 UEqn.D() + UEqn.upper() + UEqn.lower(), 而没有使用UEqn.diag()。UEqn.D()在UEqn.diag()的基础上加上了边界条件对当前网格的影响。上面方程中的[b],也不是UEqn.source(),而是和边界条件中已知的速度相关的值。我直接求解得到的结果和openfoam算出来的结果是一样的。由于没有真正的去看代码,不知道Openfoam里面有没有相关的function可以直接输出上述$[A][x]=[b]$方程组中的[A]和[b],还请各位大神指点。

                                    再次感谢大家的帮助!

                                    W 1 Reply Last reply Reply Quote
                                    • W
                                      wwzhao 教授 @mengweilm425 last edited by

                                      @mengweilm425 $[A][x]=[b]$ 中的 [A] 可以看成是 internal field 离散后的系数矩阵和 boundary conditions 对系数矩阵的影响两部分之和,其中 internal field 离散之后的系数矩阵表示为 upper(),diag() 和 lower()。boundary condition 对系数矩阵的影响则放在 internalCoeffs_ 里面。[b] 也一样,可以分为 internal field 产生的源项 source(),及 boundary conditions 对源项的影响 boundaryCoeffs_。

                                      OpenFOAM中没有直接输出 [A] 和 [b] 的函数,关于如何将 boundary conditions 的影响加到 fvMatrix 里面,你可以参考 addBoundaryDiag,addCmptAvBoundaryDiag 及 addBoundarySource 这几个函数。

                                      @wwzhao 在 关于 phi 和 fvVectorMatrix 的两个问题 中说:

                                      @李东岳 fvMatrix的source并不是线性代数方程组Ax=b中的b,而是指加在网格单元中心的力源项。

                                      上面的帖子描述不太准确,确切说 [b] 是 internal field 产生的源项 source() + boundary field 产生的源项 boundaryCoeffs_。

                                      M 1 Reply Last reply Reply Quote
                                      • M
                                        mengweilm425 @wwzhao last edited by

                                        @wwzhao 谢谢您的回复,这就明朗多了。我再好好看看这部分代码。:happy:

                                        V 1 Reply Last reply Reply Quote
                                        • V
                                          vivian @mengweilm425 last edited by

                                          @mengweilm425 您好 看了你的研究过程,觉得很好。我也想这样分步骤运行各个程序,以分清各个量含义。请问您是怎么做到的。我先编译了自己求解器,在文件中加入自己想输出的量,希望在输出中出现想看到的量。结果并没有看到想看的内容。请问您是如何。。

                                          C M 2 Replies Last reply Reply Quote
                                          • C
                                            CFD中文网 @vivian last edited by

                                            @vivian
                                            我猜测你的求解器名字没有替换?如果编译后的是icoFoam2,记得输入icoFoam2而不是原来的icoFoam..

                                            CFD中国标准用户测试帐号
                                            目前由徐笑笑登录

                                            V 1 Reply Last reply Reply Quote
                                            • V
                                              vivian @CFD中文网 last edited by

                                              @cfd-china 您好,谢谢您的回复。是替换了的,是我给他速度赋值赋错啦。。

                                              1 Reply Last reply Reply Quote
                                              • M
                                                mengweilm425 @vivian last edited by

                                                @vivian @cfd-china 看来我来晚了:happy:

                                                1 Reply Last reply Reply Quote
                                                • Wayne
                                                  Wayne last edited by

                                                  厉害了我的哥哥姐姐们

                                                  Blog: http://blog.sina.com.cn/multiphyzks
                                                  RG:https://www.researchgate.net/profile/Yan_Wang154

                                                  1 Reply Last reply Reply Quote
                                                  • First post
                                                    Last post

                                                  CFD中文网 | 东岳流体 | 京ICP备15017992号-2
                                                  论坛登录问题反馈可联系 li.dy@dyfluid.com