Navigation

    CFD中文网

    CFD中文网

    • Login
    • Search
    • 最新

    fvc::div的前后量纲关系到底是啥?

    OpenFOAM
    4
    14
    12621
    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.
    • 程
      程迪 last edited by CFD中文网

      如题,参考rhoCentralFoam.C代码178-179行。

              // --- Solve density
              solve(fvm::ddt(rho) + fvc::div(phi));
      

      phi变量根据前面的定义是

      phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg; //159行:[phi]= 密度*[aphiv_pos]
      surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);//136行:[aphiv_pos]= [phiv_pos ]
      surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());//93行:[phiv_pos]=速度*面积
      

      phi的量纲是M/T,也就是质量流量(mass flow rate, [M/T]),而不是单位面积的质量流率(mass flux, [M/L^2/T])

      而密度方程的意思应该是

      \begin{equation}
      \frac{\partial \rho}{\partial t}=\nabla\cdot\phi
      \end{equation}
      明显二者量纲是不一致的呀。

      github: chengdi123000
      网站:chengdi123000.github.io
      本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。

      程 1 Reply Last reply Reply Quote
      • 程
        程迪 @程迪 last edited by 程迪

        @程迪

        莫非OF中的div和数学$\nabla\cdot$不一样,是这个意思:

        \begin{equation}
        \text{fvc::div}(\phi_f,\alpha_f) == \frac{\sum_f{\phi_f\alpha_f}}{V}
        \end{equation}

        github: chengdi123000
        网站:chengdi123000.github.io
        本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。

        程 1 Reply Last reply Reply Quote
        • 程
          程迪 @程迪 last edited by

          @程迪

          还有一种可能是非定常项是乘以了体积V的。

          github: chengdi123000
          网站:chengdi123000.github.io
          本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。

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

            @程迪

            \begin{equation}
            \int\int\frac{\partial\mathbf{U}}{\partial t}\mathrm{d}V\mathrm{d}t=(\mathbf{U}_\mathrm{P}^{n+1}-\mathbf{U} _\mathrm{P})V _\mathrm{P}
            \end{equation}

            没有细看你一楼的alphav什么的,不过从方程离散来看,你最后的解释是正确的。

            公式显示有问题?我处理一下,

            我处理好了:sunglasses:

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

            程 2 Replies Last reply Reply Quote
            • 程
              程迪 @CFD中文网 last edited by

              @cfd-china
              我觉得你是错的,应该是
              以Euler格式为例:

              /*src\finiteVolume\finiteVolume\ddtSchemes\EulerDdtScheme\EulerDdtScheme.C*/
              //line: 102
              
              template<class Type>
              tmp<GeometricField<Type, fvPatchField, volMesh>>
              EulerDdtScheme<Type>::fvcDdt
              (
                  const GeometricField<Type, fvPatchField, volMesh>& vf
              )
              {
                  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
              
                  IOobject ddtIOobject
                  (
                      "ddt("+vf.name()+')',
                      mesh().time().timeName(),
                      mesh()
                  );
              
                  if (mesh().moving())
                  {
                      return tmp<GeometricField<Type, fvPatchField, volMesh>>
                      (
                          new GeometricField<Type, fvPatchField, volMesh>
                          (
                              ddtIOobject,
                              rDeltaT*
                              (
                                  vf()
                                - vf.oldTime()()*mesh().Vsc0()/mesh().Vsc() //动网格也是同理
                              ),
                              rDeltaT.value()*
                              (
                                  vf.boundaryField() - vf.oldTime().boundaryField()
                              )
                          )
                      );
                  }
                  else
                  {
                      return tmp<GeometricField<Type, fvPatchField, volMesh>>
                      (
                          new GeometricField<Type, fvPatchField, volMesh>
                          (
                              ddtIOobject,
                              rDeltaT*(vf - vf.oldTime()) //显然如果vf的量纲是D,那么ddt(vf)的量纲是D/T
                          )
                      );
                  }
              }
              
              

              翻译成公式是:
              \begin{equation}
              \text{fvc::ddt}(\alpha)=\frac{1}{\Delta t}(\alpha^n-\alpha^{n-1})
              \end{equation}
              这样OpenFOAM的代码中的公式量纲才是正确的。

              github: chengdi123000
              网站:chengdi123000.github.io
              本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。

              1 Reply Last reply Reply Quote
              • 程
                程迪 @CFD中文网 last edited by

                @cfd-china

                确认了,OpenFOAM嘴上说的和手底下做的不一致。

                /*src/finiteVolume/finiteVolume/fvc/fvcSurfaceIntegrate.H*/
                //line:79
                
                template<class Type>
                tmp<GeometricField<Type, fvPatchField, volMesh>>
                surfaceIntegrate
                (
                    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
                )
                {
                    const fvMesh& mesh = ssf.mesh();
                
                    tmp<GeometricField<Type, fvPatchField, volMesh>> tvf
                    (
                        new GeometricField<Type, fvPatchField, volMesh>
                        (
                            IOobject
                            (
                                "surfaceIntegrate("+ssf.name()+')',
                                ssf.instance(),
                                mesh,
                                IOobject::NO_READ,
                                IOobject::NO_WRITE
                            ),
                            mesh,
                            dimensioned<Type>
                            (
                                "0",
                                ssf.dimensions()/dimVol, //**HERE!!!!!!!!**
                                Zero
                            ),
                            extrapolatedCalculatedFvPatchField<Type>::typeName
                        )
                    );
                    GeometricField<Type, fvPatchField, volMesh>& vf = tvf.ref();
                
                    surfaceIntegrate(vf.primitiveFieldRef(), ssf);
                    vf.correctBoundaryConditions();
                
                    return tvf;
                }
                
                

                github: chengdi123000
                网站:chengdi123000.github.io
                本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。

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

                  嗯,没有显性处理时间项,而是在其他项除掉体积。

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

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

                    @cfd-china 同样的,时间项除掉了时间,也就是本贴公式中的公式(4):joking: 。
                    @李东岳 或许可以更新一下icoFoam解析,目前OpenFOAM这种颠倒的操作,可能会引起困惑。

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

                    程 2 Replies Last reply Reply Quote
                    • 程
                      程迪 @CFD中文网 last edited by

                      @cfd-china

                      今儿又发现一个吊诡的地方:fvm::div(phi,X)的量纲是什么?似乎TNND和fvc::div又不一样...

                      /*
                      OpenFOAM-dev/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C
                      */
                      //line: 90
                      
                          tmp<fvMatrix<Type>> tfvm
                          (
                              new fvMatrix<Type>
                              (
                                  vf,
                                  faceFlux.dimensions()*vf.dimensions()
                              )
                          );
                      

                      有人能解释一下么?

                      github: chengdi123000
                      网站:chengdi123000.github.io
                      本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。

                      1 Reply Last reply Reply Quote
                      • 程
                        程迪 @CFD中文网 last edited by

                        @cfd-china

                        关于OF中icoFoam的量纲测试

                        // 在合适的地方插入以下代码
                        fvVectorMatrix UEqn
                        (
                            fvm::ddt(U)
                            +fvm::div(phi,U)
                            -fvm::laplacian(nu,U)
                        )
                        Info << "Dimension of UEqn is "<<UEqn.dimensions()<<endl;
                        

                        输出显示量纲为

                        Dimension of UEqn is [0 4 -2 0 0 0 0]

                        相比之下,速度U的量纲为[0 1 -1 0 0 0 0],体积流量phi的量纲是[0 3 -1 0 0 0 0]
                        所以实际的情况似乎又是fvm::ddt项乘以了体积,fvm::laplacian项乘以了体积,fvm::div项没有乘也没有除以体积。

                        可能我前面的代码找到的地方不是fvm::ddt最终实现的地方吧。

                        注:icoFoam中的压力p是kinematic pressure. 量纲是压力除以密度[0 2 -2 0 0 0 0]

                        github: chengdi123000
                        网站:chengdi123000.github.io
                        本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。

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

                          @程迪 fvm下面的operator返回的是fvMatrix,而fvc下面的operator返回的是GeometricField。

                          fvMatrix是将物理量对单元体积分后离散得到的矩阵,由于对体积积分,所以需要将量纲乘以体积。

                          而之所以fvMatrix和GeometricField两种不同的类能相加相减等操作,是由于操作符已被重载。

                          程 1 Reply Last reply Reply Quote
                          • 程
                            程迪 @wwzhao last edited by

                            @wwzhao

                            不太理解,你的意思是对于fvm算子生成的fvMatrix和fvc生成的volScalarField相加生成fvMatrix的时候,volScalarField会被自动乘以网格体积?

                            我觉得也有这种可能,我对icoFoam里面fvMatrix::A()和H()也输出dimensions看了一下。

                            rAU=1.0/UEqn.A()的量纲是[0 0 1 0 0 0 0],所以UEqn.A()的量纲是[0 0 -1 0 0 0 0],暂,感觉A是除以了体积的对角系数,这样$AV_cU$的量纲就能和UEqn的量纲对上了。HbyA的量纲是[0 1 -1 0 0 0 0],和速度一样。所以H的量纲似乎应该是[0 1 -2 0 0 0 0],也是除以过体积的非对角项之和,然后phiHbyA的量纲是[0 3 -1 0 0 0 0]。

                            比较诡异的是H的量纲和UEqn的量纲是不一样的,OF这事儿办得忒不地道了。

                            github: chengdi123000
                            网站:chengdi123000.github.io
                            本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。

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

                              @cfd-china 多谢,有空更新一下。

                              @程迪 OpenFOAM这方面操作很绕。你输出的是矩阵的单位。以时间项为例:EulerDdtScheme.C

                              tmp<fvMatrix<Type> > tfvm
                                  (
                                      new fvMatrix<Type>
                                      (
                                          vf,
                                          vf.dimensions()*dimVol/dimTime//fvMatrix单位设定
                                      )
                                  );
                              

                              比如我们速度fvm:ddt(U),离散后矩阵单位后:速度单位/时间单位×体积单位为:Dimension of UEqn is [0 4 -2 0 0 0 0]

                              CFD高性能服务器 http://dyfluid.com/servers.html

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

                                @程迪

                                对于fvm算子生成的fvMatrix和fvc生成的volScalarField相加生成fvMatrix的时候,volScalarField会被自动乘以网格体积?

                                没错,参考 https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C#L1691

                                template<class Type>
                                Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
                                (
                                    const fvMatrix<Type>& A,
                                    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
                                )
                                {
                                    checkMethod(A, tsu(), "+");
                                    tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
                                    tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
                                    tsu.clear();
                                    return tC;
                                }
                                

                                感觉A是除以了体积的对角系数

                                没错,参考 https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C#L724

                                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;
                                }
                                
                                1 Reply Last reply Reply Quote
                                • First post
                                  Last post

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