Skip to content
  • 最新
  • 版块
  • 东岳流体
  • 随机看[请狂点我]
皮肤
  • Light
  • Cerulean
  • Cosmo
  • Flatly
  • Journal
  • Litera
  • Lumen
  • Lux
  • Materia
  • Minty
  • Morph
  • Pulse
  • Sandstone
  • Simplex
  • Sketchy
  • Spacelab
  • United
  • Yeti
  • Zephyr
  • Dark
  • Cyborg
  • Darkly
  • Quartz
  • Slate
  • Solar
  • Superhero
  • Vapor

  • 默认(不使用皮肤)
  • 不使用皮肤
折叠
CFD中文网

CFD中文网

  1. CFD中文网
  2. OpenFOAM
  3. pisoFoam计算的U,求div(U)结果为何不是严格等于0

pisoFoam计算的U,求div(U)结果为何不是严格等于0

已定时 已固定 已锁定 已移动 OpenFOAM
14 帖子 2 发布者 7.7k 浏览
  • 从旧到新
  • 从新到旧
  • 最多赞同
回复
  • 在新帖中回复
登录后回复
此主题已被删除。只有拥有主题管理权限的用户可以查看。
  • C 离线
    C 离线
    coolhhh 神
    写于 最后由 李东岳 编辑
    #1

    参考李老师的icoFoam解析,PISO算法在根据式(32)计算得到$p^*$,再根据式(31)得到 $\mathbf{U}^{**}$,此时的$\mathbf{U}^{**}$应该满足零散度条件:$\nabla\cdot{\mathbf{U}^{**} }=0$。
    但是在计算of2206自带一个算例:Decay of homogeneous isotropic turbulence,用postProcess -func "div(U)"输出不同时刻的div(U),并求最大值、最小值、均值和标准差,发现0.28,0.66 s的div(U)结果并不是一个很小的值。尝试设置tolerance=1e-10,结果也是基本一样;更改nCorrectors=5,结果也是一样。
    问题:是否理论上$\nabla\cdot\mathbf{U}^{**} =0$,但数值误差影响导致结果不等于0?

    7935a152-692d-4507-b952-3a2268ce274b-image.png


    icoFoam解析

    41907270-05d6-44bb-b3aa-abc50a602faa-image.png

    1 条回复 最后回复
  • 李东岳李 在线
    李东岳李 在线
    李东岳 管理员
    写于 最后由 编辑
    #2

    你测试一下div(phi),这个应该没问题。先试试然后我们再看div(U)的,因为其实守恒的是通量,而不是速度。

    http://dyfluid.com/index.html
    需要帮助debug算例的看这个 https://cfd-china.com/topic/8018

    C 2 条回复 最后回复
  • C 离线
    C 离线
    coolhhh 神
    在 中回复了 李东岳 最后由 编辑
    #3

    谢谢李老师,输出div(phi)就基本等于0了。这样的话如果想生成一个初始零散度U场,对于of来说其实应该考虑U计算得到的phi满足div(phi)=0,而不是div(U)=0

    1 条回复 最后回复
  • C 离线
    C 离线
    coolhhh 神
    在 中回复了 李东岳 最后由 编辑
    #4

    @李东岳 李老师,发现一个问题。

    1. 如果用myphi = fvc::flux(U)计算,然后求div(myphi)和div(U)结果基本是一致的。
    2. 问题出现在PISO算法中,pEqn.H采用phi = phiHbyA - pEqn.flux()得到phi ,U通过U = HbyA - rAU*fvc::grad(p)得到。若接着计算myphi = fvc::flux(U),此时myphi 与phi 结果不一致。
    3. 具体算例结果:
      (1)编译了mypisoFoam,参考createPhi.H创建myphi,pEqn.H中最后增加myphi = fvc::flux(U);,计算输出phi,myphi,U,接着算
    postProcess -func "div(U)"
    postProcess -func "div(phi)"
    postProcess -func "div(myphi)"
    

    (2)取部分结果展示:
    ①phi结果

    internalField   nonuniform List<scalar> 
    774144
    (
    9.2226869e-06
    -1.2808218e-05
    3.6560601e-06
    -1.834312e-06
    -1.5738475e-05
    -5.2343047e-06
    -3.6611478e-06
    

    ②myphi结果

    internalField   nonuniform List<scalar> 
    774144
    (
    9.0864005e-06
    -1.2188111e-05
    3.3814109e-06
    -1.946668e-06
    -1.5267671e-05
    -5.1722677e-06
    -3.3374742e-06
    

    ③div(U)结果

    internalField   nonuniform List<scalar> 
    262144
    (
    -0.40218979
    0.59941437
    0.85901197
    -2.6335161
    -0.55062117
    -0.53919284
    2.4044185
    1.1295023
    -1.7101542
    

    ④div(myphi)结果

    internalField   nonuniform List<scalar> 
    262144
    (
    -0.40218974
    0.59941339
    0.85901118
    -2.6335163
    -0.55062104
    -0.53919262
    2.404419
    1.1295015
    -1.710154
    

    ⑤div(phi)结果

    internalField   nonuniform List<scalar> 
    262144
    (
    -3.1704547e-08
    -6.3409098e-07
    6.3409098e-07
    8.9515991e-16
    -2.1136366e-07
    3.3818186e-07
    1.0568183e-07
    -4.2272732e-07
    -5.7068188e-07
    
    1. 结果分析:
      ①phi与myphi的值看起来差别不大。
      ②但div(phi)结果基本等于0;div(U)和div(myphi)结果一样,但不等于0。
      问题:既然phi定义为 fvc::flux(U),为何of中最后输出phi不用fvc::flux(U),而是用了phi = phiHbyA - pEqn.flux()?此时计算的通量myphi=fvc::flux(U),实际上也没满足div(myphi)=0
    李东岳李 1 条回复 最后回复
  • 李东岳李 在线
    李东岳李 在线
    李东岳 管理员
    在 中回复了 coolhhh 最后由 编辑
    #5

    @coolhhh 在 pisoFoam计算的U,求div(U)结果为何不是严格等于0 中说:

    为何of中最后输出phi不用fvc::flux(U),而是用了phi = phiHbyA - pEqn.flux()?

                    fvScalarMatrix pEqn
                    (
                        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                    );
                    pEqn.setReference(pRefCell, pRefValue);
                    pEqn.solve();
                    phi = phiHbyA - pEqn.flux();
    

    上面这几行代码里面,pEqn实际上是建立在sum(phi)= sum(phiHbyA - pEqn.flux())=0,也就是说,在sum(phi)= sum(phiHbyA - pEqn.flux())=0的时候,才有fvm::laplacian(rAU, p) == fvc::div(phiHbyA)这个形式。所以求解fvm::laplacian(rAU, p) == fvc::div(phiHbyA),sum(phi)= sum(phiHbyA - pEqn.flux())=0。

    http://dyfluid.com/index.html
    需要帮助debug算例的看这个 https://cfd-china.com/topic/8018

    C 1 条回复 最后回复
  • C 离线
    C 离线
    coolhhh 神
    在 中回复了 李东岳 最后由 编辑
    #6

    @李东岳 谢谢李老师,这样能理解,还有个小疑问。我们可以从U计算phi=fvc::flux(U),那么如果已知phi,能不能反算出U?这样的话U与phi就能完全对应上,此时计算就能满足div(U)=div(phi)。但of中后面接着用U = HbyA - rAU*fvc::grad(p)计算U,才导致了这个不对应问题,而U才是我们最后要得到的结果,但这个结果却没能满足零散度。

    1 条回复 最后回复
  • 李东岳李 在线
    李东岳李 在线
    李东岳 管理员
    写于 最后由 李东岳 编辑
    #7

    那么如果已知phi,能不能反算出U?

    就是因为这个搞不出来,所以就只是phi守恒了。各种从phi到U的算法,都涉及到各种插值之类。所以你从守恒的phi,到U(这里面有一次插值),再从U到phi(又来一次插值),误差就出来了。

    有一个算法是reconstruct(phi),不过这个里面也存在数学假设。

    按道理来说守恒的场就是phi场。我们不能说速度是守恒的。我们只能说通量是守恒的。类似只能说能量是守恒的,不能说温度是守恒的。

    http://dyfluid.com/index.html
    需要帮助debug算例的看这个 https://cfd-china.com/topic/8018

    C 2 条回复 最后回复
  • C 离线
    C 离线
    coolhhh 神
    在 中回复了 李东岳 最后由 编辑
    #8

    理解了,谢谢李老师:146:

    1 条回复 最后回复
  • C 离线
    C 离线
    coolhhh 神
    在 中回复了 李东岳 最后由 编辑
    #9

    李老师,我做了简单的对比:

    1. 下面计算的div(U) 与 div(phi) 结果基本相等,区别是小数点后五六位
      ①直接求div(U)
      ②已知U -> phi = fvc::flux(U) -> div(phi)

    2. 接着对比reconstruct(phi)。编译mypisoFoam,把pEqn.H中计算U修改为

    U = fvc::reconstruct(phi);
    // U = HbyA - rAU*fvc::grad(p);
    

    同样计算计算of2206自带一个算例:Decay of homogeneous isotropic turbulence,模拟结果:
    (1)对比div(U) 与 div(phi)的均值和标准差。采用reconstruct(phi),得到的div(U) 的均值和标准差要小一点。这里有个注意点是,of计算的div(phi) ,基本数值都很小,但会有一个比较大的数,下面的计算结果是剔除了这个最大数的计算结果。
    e3b4395e-fa8a-4a42-9696-baa7ee84d18c-image.png

    (2)对比湍流特性。可以看出reconstruct(phi)结果要更差,还是原本of的U = HbyA - rAU*fvc::grad(p)计算准确
    68f537b7-f131-441e-a855-954ad38eddeb-image.png

    1 条回复 最后回复
  • 李东岳李 在线
    李东岳李 在线
    李东岳 管理员
    写于 最后由 编辑
    #10

    对比湍流特性。可以看出reconstruct(phi)结果要更差,还是原本of的U = HbyA - rAU*fvc::grad(p)计算准确

    对比你下面那两个图,我觉得你可以把这部分总结一下写在文章里面,当做appendix。这是2种计算方法,从方程上来说后者是一致的。从物理上来说前者的U看起来更加守恒。不过从图上来看,reconstruct耗散掉了一部分没有后者好。

    你这个做的很细,值得一写。

    http://dyfluid.com/index.html
    需要帮助debug算例的看这个 https://cfd-china.com/topic/8018

    C 1 条回复 最后回复
  • C 离线
    C 离线
    coolhhh 神
    在 中回复了 李东岳 最后由 编辑
    #11

    谢谢李老师,文章已经发表了,只是好奇做个对比。在想的问题是,一直以来各种湍流生成方法都在追求实现生成零散度的湍流场U,我们直接构造的是0时刻的U,而不是phi。

    根据测试结果看,如果严格让div(U)=0, 此时U -> phi = fvc::flux(U) -> div(phi)得到的div(phi)也基本等于0。但这种情况下0时刻的湍流特性会变化的比较大,结果会变差。

    of在后面时间步计算输出的U实际上不满足div(U)=0,所以想法是0时刻提供的初始U也不需要强制要求div(U)=0,这样得到的结果可能才是更合理的。

    1 条回复 最后回复
  • 李东岳李 在线
    李东岳李 在线
    李东岳 管理员
    写于 最后由 编辑
    #12

    of在后面时间步计算输出的U实际上不满足div(U)=0,所以想法是0时刻提供的初始U也不需要强制要求div(U)=0,这样得到的结果可能才是更合理的。

    你目前0时间步的U与phi,是前者散度为0,还是后者。

    不过我感觉0时间步的问题不大,稍微步进一个时间步,后面就都是div(phi)=0了。

    厉害啊,小伙子,我记得你好像做建筑领域的,现在在做DNS么

    http://dyfluid.com/index.html
    需要帮助debug算例的看这个 https://cfd-china.com/topic/8018

    C 1 条回复 最后回复
  • C 离线
    C 离线
    coolhhh 神
    在 中回复了 李东岳 最后由 编辑
    #13

    0时间步的U,是用湍流生成方法直接生成的,可以计算得到div(U)。然后就可以计算div(phi): U -> phi = fvc::flux(U) -> div(phi)。测试结果是,这种方式计算的div(U)和div(phi)结果基本一样的,区别是小数点后五六位。此时div(U)=div(phi)=0

    稍微步进一个时间步,of输出的U和phi,就会出现div(U)不等于0,div(phi)=0

    是做建筑领域的,DNS搞不了,还是做的LES

    1 条回复 最后回复
  • 李东岳李 在线
    李东岳李 在线
    李东岳 管理员
    写于 最后由 编辑
    #14

    0时间步的U,是用湍流生成方法直接生成的,可以计算得到div(U)。然后就可以计算div(phi): U -> phi = fvc::flux(U) -> div(phi)。测试结果是,这种方式计算的div(U)和div(phi)结果基本一样的,区别是小数点后五六位。此时div(U)=div(phi)=0

    那这个算法很屌。不错。

    http://dyfluid.com/index.html
    需要帮助debug算例的看这个 https://cfd-china.com/topic/8018

    1 条回复 最后回复
  • C coolhhh 被引用 于这个主题
  • C coolhhh 被引用 于这个主题

  • 登录

  • 登录或注册以进行搜索。
  • 第一个帖子
    最后一个帖子
0
  • 最新
  • 版块
  • 东岳流体
  • 随机看[请狂点我]