只给热传导单独的τ,粘性加热却被锁在Pr=1 — 双分布LB中的准平衡g*

Pr变了,温度为什么没变#

设想一个Couette通道。下壁静止,上壁以速度$U$滑动。两壁温度相同。中间依然会升温,因为粘性把动能转化为热(粘性加热,viscous heating)。中心的温升与Prandtl数(Pr,动量扩散与热扩散之比)严格成正比。

本文讨论如何在格子Boltzmann(LB,Lattice Boltzmann)热模型中改变Pr。把分布函数分成两组,各给一个弛豫时间,Pr看起来就能随意调节。真正运行后会发现,热传导变了,粘性加热却一直停在Pr=1的值。读者将了解其原因,并得到只需一行的解决办法:准平衡(quasi-equilibrium)$g^*$。

拖动 Pr 滑块,应力时钟和热传导时钟的间隔会变,粘性功通量贴在哪条曲线上也会一起变。

把Pr滑块放在0.71,橙色(应力)和绿色虚线(g*模型的粘性功通量)始终重合。红色(简单双τ模型)跟着紫色(热传导)的时钟走,与应力分道扬镳。把Pr移到1,四条线合并成两条。

一组分布函数只认得Pr=1#

BGK碰撞算子以同一个速率$\omega$把所有非平衡矩拉回平衡。应力和热通量按同一个时钟弛豫。因此粘性系数和导热系数被绑在同一个$\omega$上。连续BGK方程给出Pr=1,原因就在这里。

空气的Pr为0.71,水约为7。锁定在Pr=1的模型连边界层厚度之比都对不上。补救办法有两条。一条是在碰撞算子内部只让热通量矩以不同速率弛豫(MRT、Shakhov型)。另一条是把能量交给一个独立的分布函数。本文走第二条路。

两组:f负责质量和动量,g负责总能#

在双分布模型中,$f_i$输运密度和动量,$g_i$输运总能$E=\rho(c_vT+u^2/2)$($i=0,\dots,Q-1$)。动理学方程如下。

$\omega_\mu$决定粘性(运动粘度$\nu$),$\omega_\kappa$决定热传导(热扩散率$\alpha$)。$g_i^{eq}$是局部平衡,$g_i^*$是准平衡分布。

两个系数与扩散系数的关系如下。

$c_s^2=1/3$是格子声速的平方。记$\tau_\mu=1/\omega_\mu$、$\tau_\kappa=1/\omega_\kappa$,则$\mathrm{Pr}=(\tau_\mu-\tfrac12)/(\tau_\kappa-\tfrac12)$。

代码把 $\tau_\mu$ 固定为 1。于是 $\tau_\mu-\tfrac12=\tfrac12$,$\tau_\kappa=\tfrac12+(\tfrac12)/\mathrm{Pr}$。表里四个值就是这样来的。Pr=0.5 时 $\tau_\kappa=1.5$。Pr=0.71 时加上 $0.5/0.71\approx0.704$,$\tau_\kappa\approx1.204$。Pr=1 时 $\tau_\kappa=1$。Pr=2 时 $\tau_\kappa=0.75$。因为 $\omega=1/\tau$,Pr=2 的热弛豫是 $\omega_\kappa\approx1.333$,黏性是 $\omega_\mu=1$。粘性功应该留在慢的那一侧。简单双 τ 却把这项放到更快的 $\omega_\kappa$ 上。只有 Pr=1 那一行两个 $\omega$ 相等,两列都是 1.000。只看这一行,错误会一直安静。

粘性功藏在能量通量里

在本文的简化模型中,$g^{eq}$只规定三个矩。零阶为$E$,一阶(能量通量)为$(E+p)\mathbf{u}$,二阶为$c_s^2E\,\mathbf{I}+E\,\mathbf{u}\mathbf{u}$。一阶矩的非平衡部分$\mathbf{Q}$就是实际的热通量。

在Chapman–Enskog展开(以小Knudsen数求解非平衡的展开)的一阶,$\mathbf{Q}\approx-\tau\,\nabla\!\cdot\!\mathbf{R}^{eq}$。这里$\nabla(c_s^2E)$分成两部分。

第一项是Fourier热传导。第二项是粘性功$-\mathbf{u}\cdot\boldsymbol{\sigma}$的原料。它与$f$的非平衡应力$\boldsymbol{\Pi}^{neq}=\mathbf{P}-\mathbf{P}^{eq}\approx-\tau_\mu\rho c_s^2(\nabla\mathbf{u}+\nabla\mathbf{u}^T)$形状完全相同。问题出在前面的$\tau$。

如果只以$\omega_\kappa$把$g$拉向平衡,两部分都会带着$\tau_\kappa$出来。粘性功的系数正比于$\alpha$而不是$\nu$。热传导是对的,粘性加热却停留在Pr=1时的值。第一个viz中的红线就是这种情况。

准平衡g*:只把通量绑到应力的时钟上#

解决办法是$g^*$定义中的一行。$g^*$与$g^{eq}$的所有矩相同,只有能量通量不同。

$\mathbf{P}$是由$f$测得的实际动量通量,$\mathbf{P}^{eq}$是其平衡值。

对动理学方程取一阶矩,得到两个非平衡矩的单步更新。

$\mathbf{Q}$本身以$\omega_\kappa$弛豫,所以热传导服从$\alpha$。同时,若$\mathbf{Q}=\mathbf{u}\cdot\boldsymbol{\Pi}^{neq}$,下一步仍有$\mathbf{Q}'=\mathbf{u}\cdot\boldsymbol{\Pi}'^{neq}$。粘性功部分搭上了应力的时钟$\omega_\mu$。用Chapman–Enskog计算也一样:带$\tau_\kappa$的粘性功与$(\omega_\kappa-\omega_\mu)$项相消,只剩下$\tau_\mu$。

当$\omega_\kappa=\omega_\mu$时,$g^*$项消失。此时模型与单弛豫模型相同,Pr=1。$g^*$是只在Pr$\ne1$时才起作用的修正。只用一组分布的热格子,把 τ 加大,Pr 也还是贴在 1。这和那两个常数锁死的地方是同一堵墙。分布分成两组,如果不把粘性功的时钟挪过去,墙还在。问题只有一个:哪个矩以什么速率弛豫。

代码:同一个Couette通道,两种g#

解析解很简单。定常能量方程为$kT''+\mu(u')^2=0$。速度线性分布时,温度为抛物线。

中心最大值为$\theta_{max}=\mathrm{Pr}$。下面的代码把D2Q9在$x$方向设为一个周期单元,只沿$y$求解。壁面采用非平衡外推(Guo)边界。为简化起见,$f$使用等温平衡,能量只单向耦合。在低马赫($U=0.1$)Couette流中压力均匀,这一简化不改变结果。

import numpy as np
 
# D2Q9,x方向一个周期单元 → 沿y的一维数组
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w = np.array([4/9] + [1/9]*4 + [1/36]*4)
CV, T0 = 1.0, 1.0
 
def feq_lattice(rho, ux):
    cu = cx[:, None]*ux
    return w[:, None]*rho*(1 + 3*cu + 4.5*cu**2 - 1.5*ux**2)
 
def geq_energy(E, rho, ux):
    p = rho*(E/rho - 0.5*ux**2)/CV          # p = rho*T
    qx = (E + p)*ux
    cu = cx[:, None]*ux
    return w[:, None]*(E + 3*cx[:, None]*qx + 4.5*E*cu**2 - 1.5*E*ux**2)
 
def couette_heating(Pr, quasi, N=21, U=0.1, tau_mu=1.0, steps=8000):
    om_mu = 1/tau_mu
    om_ka = 1/(0.5 + (tau_mu - 0.5)/Pr)
    rho = np.ones(N); ux = np.zeros(N)
    E = rho*CV*T0*np.ones(N)
    f = feq_lattice(rho, ux); g = geq_energy(E, rho, ux)
    uw = np.zeros(N); uw[-1] = U
    for _ in range(steps):
        rho = f.sum(0); ux = (cx[:, None]*f).sum(0)/rho; E = g.sum(0)
        fe = feq_lattice(rho, ux); ge = geq_energy(E, rho, ux)
        Pxy = ((cx*cy)[:, None]*(f - fe)).sum(0)
        f = f + om_mu*(fe - f)
        if quasi:   # g*:在能量通量上加u·(P - P_eq)的准平衡
            gs = ge + w[:, None]*3*cy[:, None]*(ux*Pxy)
            g = g + om_mu*(ge - g) + (om_ka - om_mu)*(gs - g)
        else:       # 简单双τ:g只以ω_κ弛豫到平衡
            g = g + om_ka*(ge - g)
        for k in range(9):
            f[k] = np.roll(f[k], cy[k]); g[k] = np.roll(g[k], cy[k])
        for wall, nb in ((0, 1), (N-1, N-2)):   # 非平衡外推壁面
            r = f[:, nb].sum(); un = (cx*f[:, nb]).sum()/r; En = g[:, nb].sum()
            fn = feq_lattice(np.array([r]), np.array([un]))[:, 0]
            gn = geq_energy(np.array([En]), np.array([r]), np.array([un]))[:, 0]
            Ew = r*(CV*T0 + 0.5*uw[wall]**2)
            f[:, wall] = feq_lattice(np.array([r]), np.array([uw[wall]]))[:, 0] + f[:, nb] - fn
            g[:, wall] = geq_energy(np.array([Ew]), np.array([r]), np.array([uw[wall]]))[:, 0] + g[:, nb] - gn
    rho = f.sum(0); ux = (cx[:, None]*f).sum(0)/rho
    T = (g.sum(0)/rho - 0.5*ux**2)/CV
    theta = (T - T0)/(U**2/(8*CV))
    return theta.max(), ux
 
print(" Pr   exact   2-tau(g_eq)   g* quasi-eq")
for Pr in (0.5, 0.71, 1.0, 2.0):
    a, _ = couette_heating(Pr, False)
    b, u = couette_heating(Pr, True)
    print(f"{Pr:4.2f}  {Pr:6.3f}   {a:9.3f}    {b:9.3f}")
print(f"u(mid)/U = {u[10]/0.1:.4f}")

 Pr   exact   2-tau(g_eq)   g* quasi-eq
0.50   0.500       1.000        0.498
0.71   0.710       1.000        0.710
1.00   1.000       1.000        1.000
2.00   2.000       1.000        1.997
u(mid)/U = 0.5000

在浏览器里直接运行上面的代码。两个模型在同一个通道中并排推进。

拖动Pr滑块,绿色曲线的峰值跟随虚线解析值Pr,红色则贴在1.000不动。改变$U$或$\tau_\mu$,红色的$\theta_{max}$也不会离开1。不妨试试输出块中没有的Pr=3。

1.000这个数字说明了什么#

简单双τ一列无论Pr为多少都是1.000。这不是误差,而是结构。粘性功系数正比于$\tau_\kappa-\tfrac12$,导热系数也正比于$\tau_\kappa-\tfrac12$。在二者之比$\theta_{max}$中,$\tau_\kappa$恰好约掉。对空气(Pr=0.71),两壁之间的温升被高估41%。对Pr=7的水类流体,则被低估为七分之一。

$U=0.1$、$c_v=1$ 时,分母是 $U^2/(8c_v)=0.00125$。空气的真实 $\Delta T$ 为 $0.71\times0.00125\approx8.9\times10^{-4}$,简单双 τ 给出 $1.25\times10^{-3}$。格子温度是 1.0009 和 1.0013。只保留两位小数,两者都是 1.00。Pr=7 时真值是 $0.00875$,模型是 $0.00125$。1.009 对 1.001。要看 $\theta$,1 和 7 才会分开。

$g^*$一列为0.498、0.710、1.000、1.997。0.498 比 0.500 小 0.4%,1.997 比 2 小 0.15%。0.710 与 exact 相同。21 个格点上剩下的差别来自壁面和分辨率。中心速度$u/U=0.5000$完全相同,因为两个模型使用同一个$f$。温度上的错误不会出现在速度场里。

走向可压缩时,f也要加准平衡#

本文的$f$使用等温平衡。在可压缩双分布模型中,$f$一侧也要引入准平衡$f^*$。D2Q9的三阶矩有$c_{ix}^3=c_{ix}$这一缺陷。温度变化或马赫数增大时,应力中会出现误差项。$f^*$携带修正这一缺陷的压力张量。修正动量通量的是$f^*$,给能量通量加上粘性功的是$g^*$。

原理相同。守恒量(密度、动量、总能)由平衡和准平衡共享。唯一的不同是需要单独设定弛豫速率的那一个矩。非理想流体的压力项放在哪里,归根结底也是这个矩的位置。

热LB代码单独设定τκ时要检查的一行#

打开能量分布的碰撞式,看$\omega_\kappa$乘在哪里。如果整个$g$都以$\omega_\kappa$拉向平衡,粘性加热就停在Pr=1的值。检查只需跑一次Couette。若$\theta_{max}$不跟随Pr,说明缺了准平衡项。

只在$\omega_\kappa=\omega_\mu$下验证过,这个错误会安静地留下。热 LB 的检查里要放一个 Pr$\ne1$ 的粘性加热算例。静止流体里的温度扩散转不动这一项。把摩擦和传热绑在一个比值上的雷诺比拟,前提就是这个比值拴在 Pr 上。$\theta_{max}$ 不跟随 Pr,这个前提在格子里已经破了。