Navigation

    CFD中文网

    CFD中文网

    • Login
    • Search
    • 最新

    HbyA,phiHbyA,fvc::div(phiHbyA)计算错误问题

    OpenFOAM
    2
    16
    5461
    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.
    • yhdthu
      yhdthu 讲师 last edited by yhdthu

      最近做pEqn.H中的演算,输出了HbyA,phiHbyA和fvc::div(phiHbyA),通过手动计算进行验证

      首先测试了一个简单的331的立方体模型,由HbyA计算phiHbyA,再进行求和并除以体积求div,计算结果和输出的结果一致,没问题。

      换了个模型,吊诡的事情出现了,由HbyA计算phiHbyA没问题,但是由phiHbyA进行求和算div居然是错的?!我很奇怪难道连加法也能算错?

      以下是我的一个结果,测试的cell编号是47794:

      手动计算的div:-6.3169e-3
      输出的div:1.90343e-6

      且不说正负号,cell体积的量级就是10^-8,手动求和phiHbyA的量级是10^-10,check了一下代码,也没有其它干扰了,实在想不通内部是怎么算的

      长风破浪会有时,直挂云帆济沧海

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

        以下是我的一个结果,测试的cell编号是47794

        不太清楚你怎么计算的。你要测试什么?计算$\nabla\cdot\mathbf{U}=0$?

        phiHbyA进行求和算div居然是错的?

        在OpenFOAM里面,方程$\nabla\cdot\mathbf{U}=0$的植入就是fvc::div(phi),对应:fvc::surfaceIntegrate(phi),对应面加和。

        我对这个也感兴趣,能否提供更详细信息。

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

        yhdthu 1 Reply Last reply Reply Quote
        • yhdthu
          yhdthu 讲师 @李东岳 last edited by yhdthu

          @李东岳 前辈好,是这样,我打算测试压力泊松方程源项对计算结果的影响

          fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
          

          以上是压力方程,我想看看**源项fvc::div(phiHbyA)**的计算结果和我自己手动计算的是否一样

          所以下一步我需要求phiHbyA,我看了fvc::surfaceIntegrate的源码,fvc::div(phiHbyA)的计算就是将一个cell中不同face的通量phiHbyA相加,之后除以体积mesh.V()

          关于phiHbyA的计算,查看代码如下,

          surfaceScalarField phiHbyA
          (
              "phiHbyA",
              fvc::flux(HbyA)
            + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
          );
          

          有两个部分组成,一个是通过volVectorField类型的HbyA中心插值到face上再与面矢量做点积,第二个修正项是通过上一时间步对本时间的修正,具体算法代码中有

          目前采用的都是2D的算例,一个cell有4个neighbour,具体细节如下:

          对于简单算例,几何模型如下,是一个3 * 3 * 1的一个立方体,选择4号单元进行验证,具体数据如图

          0_1500374893743_12.jpeg
          0_1500374162407_4.jpeg

          由HbyA手动计算得得phiHbyA与输出的结果一致,再由phiHbyA计算fvc::div(phiHbyA)也没问题;

          对于网格量比较大的复杂问题

          0_1500373676075_3.jpeg
          0_1500373688997_1.jpeg

          由HbyA手动计算得得phiHbyA与输出的结果一致但是,由phiHbyA求和除以mesh.V()所得的结果,与输出的fvc::div(phiHbyA)相差了好几个量级!也就是之前列出来的结果了。

          求前辈指教,我现在不懂为啥这个求和还能求错了?或者是求这个fvc::div(phiHbyA)时根本用的就不是phiHbyA求和来的?

          长风破浪会有时,直挂云帆济沧海

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

            @yhdthu

            #include "fvCFD.H"
            
            // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
            
            int main(int argc, char *argv[])
            {
                #include "setRootCase.H"
                #include "createTime.H"
                #include "createMesh.H"
            
                volVectorField U
                (
                    IOobject
                    (
                        "U",
                        runTime.timeName(),
                        mesh,
                        IOobject::MUST_READ,
                        IOobject::AUTO_WRITE
                    ),
                    mesh
                );
            
            	surfaceScalarField phi
            	(
                    IOobject
                    (
                        "phi",
                        runTime.timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                    ),
                   	fvc::interpolate(U) & mesh.Sf()
                );
            
            	Info << "fvc::interpolate(U) " << fvc::interpolate(U) << nl;
            
            	Info << "phi is " << phi << nl;
            
            	Info << "mesh surface vector is " << mesh.Sf() << nl;
            
            	Info << "fvc::div(phi) is " << fvc::div(phi) << nl;
            
                return 0;
            }
            

            请运行上述代码,并把U初始化为:

            dimensions      [0 1 -1 0 0 0 0];
            
            internalField   nonuniform List<vector> 
            9
            (
            (0 1 0)
            (0 2 0)
            (0 3 0)
            (0 4 0)
            (0 5 0)
            (0 6 0)
            (0 7 0)
            (0 8 0)
            (0 9 0)
            )
            ;
            
            
            //internalField   uniform (0 1 0);
            
            boundaryField
            {
                wall
                {
                    type            fixedValue;
                    value           uniform (0 1 0);
                }
                defaultFaces
                {
                    type            empty;
                }
            }
            

            我没看出计算错误的问题。fvc::div就是对速度做体积分,通过加和的方式计算。

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

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

              Yes, 我觉得这个算例非常好,因为OKSS2上来就是高斯积分计算,也算是我对你这个感兴趣的动力,我会把这个放在课上分析。

              由phiHbyA求和除以mesh.V()所得的结果,与输出的fvc::div(phiHbyA)相差了好几个量级!

              用上面那个简单试试,可能你那个太复杂了,小数点很多,数也不工整。客观来讲,不会出现这个现象。

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

              yhdthu 1 Reply Last reply Reply Quote
              • yhdthu
                yhdthu 讲师 @李东岳 last edited by

                @李东岳 前辈 您说的这个算例我没异议,这个就和我的那个简单例子一样的,我目前已经掌握HbyA到fvc::div(phiHbyA)的计算方法。

                但是,问题来了,我用相同的计算方法去计算我的第二个例子的某个cell,发现求和的时候出现了错误。

                即,我手动计算的fvc::div(phiHbyA)为:-6.3169e-3
                但是输出的结果却是1.90343e-6

                如果只看phiHbyA求和,其结果是10^-10量级,cell的体积是10^-7量级,两个相除就是10^-3次方量级了;

                输出的fvc::div(phiHbyA)的结果居然是10^-6量级,这可是除完体积之后的值,也就是说之前的phiHbyA的求和量级是10^-13量级

                前辈,这就是我想表达的意思,为什么从phiHbyA就算不到fvc::div(phiHbyA)了??

                长风破浪会有时,直挂云帆济沧海

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

                  所以就是说:你实际上已经知道手动计算的方式,但是你算的和代码结果不一样。如果可以的话,请把4号网格4个面的phiHbyA输出,我看看你是怎样 手动计算的fvc::div(phiHbyA) 的,就是说:把你计算-6.3169e-3的过程写出来。

                  输出的fvc::div(phiHbyA)的结果居然是10^-6量级

                  这个非常正常,匹配你的压力方程残差。
                  你可以看出来,这只不过是加减法问题,可能你手动计算存在一个非常小的bug。

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

                  yhdthu 2 Replies Last reply Reply Quote
                  • yhdthu
                    yhdthu 讲师 @李东岳 last edited by

                    @李东岳 前辈,我的意思是:

                    我计算那个方形的4号cell全都没问题,以上第一个方形网格的算例给出了我所查找的所有参数;

                    有问题的是,用相同的计算方法,去计算第二个算例的47794号cell,从phiHbyA到fvc::div(phiHbyA)的计算出了问题。

                    我明白量级什么的无所谓,最重要的是数得算的对,但是现在手动算的和程序算的根本就对不上了,而且差了n个量级,这就说不过去了:crying:

                    长风破浪会有时,直挂云帆济沧海

                    1 Reply Last reply Reply Quote
                    • yhdthu
                      yhdthu 讲师 @李东岳 last edited by yhdthu

                      @李东岳 前辈 具体的计算过程如下,这里只说有错误的第二个算例(cell_47794)

                      需要进行计算的场输出为:
                      HbyAVol = HbyA;
                      phiHbyASuf = phiHbyA;
                      sourcePoissionU = fvc::div(phiHbyA);

                      根据输出的phiHbyA,4个面的值如下:

                      E:0.000125002
                      W:0.000124669
                      N:-1.94136e-05
                      S:-1.908e-05
                      F 和 B 为empty,0

                      加和:E - W + N - S = 6e-10

                      mesh.V() = 9.49835e-8

                      fvc:div(phiHbyA) = (E - W + N - S) / mesh.V() = -0.0063169

                      但是,输出的fvc::div(phiHbyA) = 1.90343e-6

                      长风破浪会有时,直挂云帆济沧海

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

                        E:0.000125002
                        W:0.000124669
                        N:-1.94136e-05
                        S:-1.908e-05
                        

                        这个找的正确么

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

                        yhdthu 1 Reply Last reply Reply Quote
                        • yhdthu
                          yhdthu 讲师 @李东岳 last edited by

                          @李东岳 前辈,这个百分百没问题,我是用paraview一个一个面点出来看的

                          长风破浪会有时,直挂云帆济沧海

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

                            你是怎么在paraview查看face的值的?

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

                            yhdthu 1 Reply Last reply Reply Quote
                            • yhdthu
                              yhdthu 讲师 @李东岳 last edited by

                              @李东岳 通过paraview找到cell编号,查找polymesh里面的owner和neighbor文件定位行号,打开surfacefield对应找就好啦~

                              长风破浪会有时,直挂云帆济沧海

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

                                方法是对的,不过为什么错了?我也很奇怪。

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

                                yhdthu 2 Replies Last reply Reply Quote
                                • yhdthu
                                  yhdthu 讲师 @李东岳 last edited by

                                  @李东岳 更奇怪的是,程序前半段计算phiHbyA都是对的,但是之后求div的时候(加法)却错了,难道加法也能算错:upset:

                                  长风破浪会有时,直挂云帆济沧海

                                  1 Reply Last reply Reply Quote
                                  • yhdthu
                                    yhdthu 讲师 @李东岳 last edited by

                                    @李东岳 前辈,解决了,原来是输出精度不够导致的,我的锅:lol:

                                    长风破浪会有时,直挂云帆济沧海

                                    1 Reply Last reply Reply Quote
                                    • First post
                                      Last post

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