气液流动过程中温度直接求解问题

[李东岳]

你这个温度跟压力是耦合在一起的么

[upc_ngh]

@李东岳 对的，东岳老师，我是先求解动量方程，压力修正方程之后，最后一步再求解温度场，压力这个时候已经是已知的吧，然后就各种发散

[upc_ngh]

@李东岳 东岳老师，上次温度方程求解错误主要是由于两个原因，一个是边界条件给错了，还一个是温度方程和动量、连续相方程耦合求解了，经过更改已经可以正常运行。但目前在验证算例时发现和实验值存在较大的误差。
验证算例如下： 一根长为 12 m， 直径为 1 m的垂直管， 管道上端为入口， 管道下端为出口，出口压力是 1e5Pa， 气液相分别为空气和水，温度均为 25 ℃ ．入口处 x = 0， 出口处 x = 12 m．初始状态如图 1(a)所示， 气相相含率 0. 2， 气相速度 Vg=0；液相相含率 0. 8， 液相速度 Vl=10 m/s。使用双欧拉求解。耦合求解动量及连续性方程（参考瞬态可压缩压力泊松方程），收敛后求解相率方程，进行下一时间步。图一是示意图，图二是解析解的气相分数及液相流速，图三为求解出的压力、持气率、液相流速，可以看出波动很大，尤其持气率波动尤其大，已经排除网格密度的原因。能力有限，查找一周的错误了，实在找不到什么别的原因了，还望东岳老师和前辈指点一二。

气液相动量方程：

fvVectorMatrix UEqn1
           (
                fvm::ddt(alpharho1, U1) + fvm::div(alphaphi1, U1)
                == -pGrad1-Sg
           );
           UEqn1.solve();

           fvVectorMatrix UEqn2
           (
                fvm::ddt(alpharho2, U2) + fvm::div(alphaphi2, U2)
                   == -pGrad2-Sl
            );
           UEqn2.solve();

连续性方程

psi1 = 1/Cg/Cg;
           psi2 = rho2/p;

           volScalarField rAU1(1.0/UEqn1.A());
           surfaceScalarField rhorAUf1("rhorAUf1", fvc::interpolate(alpharho1*rAU1/rhog));
           volVectorField HbyA1(constrainHbyA(rAU1*UEqn1.H(), U1, p));

           surfaceScalarField phiHbyA1
             (
                    "phiHbyA1",
                 (
                   fvc::interpolate(psi1*alpha1/rhog)
                   *fvc::flux(HbyA1)
                 ) 
              );

           volScalarField alphapsi1("alphapsi1",alpha1*psi1/rhog);

   
           volScalarField rAU2(1.0/UEqn2.A());

           surfaceScalarField rhorAUf2("rhorAUf2", fvc::interpolate(alpharho2*rAU2/rho2));

           volVectorField HbyA2(constrainHbyA(rAU2*UEqn2.H(), U2, p));

           surfaceScalarField phiHbyA2
             (
                    "phiHbyA2",
                  (
                      fvc::interpolate(psi2*alpha2/rho2)
                      *fvc::flux(HbyA2)
                  ) 
              );

           volScalarField alphapsi2("alphapsi2",alpha2*psi2/rho2);

  
   
fvScalarMatrix pEqn
        (
            fvm::ddt(alphapsi1, p)
          + fvm::ddt(alphapsi2, p)
          + fvm::div(phiHbyA1, p)
          + fvm::div(phiHbyA2, p)
          - fvm::laplacian(rhorAUf1, p)
          - fvm::laplacian(rhorAUf2, p)
        );
  pEqn.solve();
  rho1 = p*psi1;
  rho1.correctBoundaryConditions();

surfaceScalarField phi1 = fvc::interpolate(rho1*U1) & mesh.Sf();
   surfaceScalarField phi2 = fvc::interpolate(rho2*U2) & mesh.Sf();  

   fvScalarMatrix alpharhoEqn1
    (
        fvm::ddt(rho1,alpha1)
      + fvm::div(phi1,alpha1)
    );
   
fvScalarMatrix alpharhoEqn2
    (
        fvm::ddt(rho2,alpha2)
      + fvm::div(phi2,alpha2)
    );

 
   alpharhoEqn1.solve();
   alpharhoEqn2.solve();
  
  
   alpha1 = alpha1/(alpha1+alpha2);
   alpha2 = 1 - alpha1;
   alpha2.correctBoundaryConditions();

[李东岳]

前几天谈融资去了，今天刚回来，你这个我把这几天积压的事情弄完 给你看一下

[upc_ngh]

@李东岳 哇，谢谢东岳老师，非常感谢感谢

[upc_ngh]

@李东岳 东岳老师，上次那个程序您帮忙看看，有结果了吗

[李东岳]

你是不是之前来参加过OKS课？

你这个看起来简单一句两句说不清。我今天面试了一天。周四晚上9.30如果你方便，我安排个腾讯会议，你给我介绍一下你的问题。

[upc_ngh]

@李东岳 对的李老师，我去年参加过OKS课程，周四晚上有时间的，那就提前谢谢李老师了！

[李东岳]

@upc_ngh

李东岳 邀请您参加腾讯会议
会议主题：李东岳预定的会议
会议时间：2021/12/30 21:30-22:00 (GMT+08:00) 中国标准时间 - 北京

点击链接入会，或添加至会议列表：
https://meeting.tencent.com/dm/fevw9LEt0xWa

#腾讯会议：381-992-089

复制该信息，打开手机腾讯会议即可参与

[upc_ngh]

@李东岳 收到，谢谢东岳老师！