【COMSOL】修改内置物理场方程 模拟水沸腾的蒸汽泡运动

前言

在COMSOL中通过弱贡献（weak contribution）节点，可以向任何内置模型方程中添加约束。这里，以模拟水沸腾的蒸汽泡运动为例，来讲一讲当我们研究的方程与模块方程不一样时，应该怎么增加新的贡献项。

物理场

涉及到水的非等温流动，液相和气相的转化。添加三个模块，层流、相场和传热。主要的控制方程为：

流体：

ρ ∂ u ∂ t + ρ ( u ⋅ ▽ ) u = ▽ ⋅ [ − p I   + K ] + F + ρ g \rho \frac{\partial \textbf{u}}{\partial t}+\rho \left ( \textbf{u} \cdot \bigtriangledown \right )\textbf{u}=\triangledown \cdot \left [ -p\textbf{I }+\textbf{K}\right ]+\textbf{F}+\rho \textbf{g} ρ∂t∂u​+ρ(u⋅▽)u=▽⋅[−pI +K]+F+ρg (1)

∂ ∂ t ρ + ▽ ⋅ ( ρ u ) = 0 \frac{\partial }{\partial t}\rho+\bigtriangledown \cdot (\rho \textbf{u})=0 ∂t∂​ρ+▽⋅(ρu)=0 (2)

相场：

∂ ϕ ∂ t + u ⋅ ▽ ϕ = ▽ ⋅ γ λ ε 2 ▽ ψ \frac{\partial \phi }{\partial t}+\textbf{u}\cdot \bigtriangledown \phi =\bigtriangledown \cdot \frac{\gamma \lambda }{\varepsilon ^2}\bigtriangledown \psi ∂t∂ϕ​+u⋅▽ϕ=▽⋅ε2γλ​▽ψ (3)

ψ = − ▽ ⋅ ε 2 ▽ ϕ + ( ϕ 2 − 1 ) ϕ \psi =-\bigtriangledown \cdot \varepsilon ^2\bigtriangledown \phi +(\phi ^2-1)\phi ψ=−▽⋅ε2▽ϕ+(ϕ2−1)ϕ (4)

传热
ρ C p ∂ T ∂ t + ρ C p u ⋅ ▽ T + ▽ ⋅ q = 0 \rho C_p\frac{\partial T}{\partial t}+\rho C_p\textbf{u}\cdot \bigtriangledown T+\bigtriangledown\cdot \textbf{q}=0 ρCp​∂t∂T​+ρCp​u⋅▽T+▽⋅q=0 (5)

以相场方程为例，Cahn - Hilliard 方程中不含源项，若考虑液气相变，需要增加源项：

S f = m ˙ δ ( V f , v ρ v + V f , L ρ L ) S_f=\dot{m}\delta \left ( \frac{V_{f,v}}{\rho_v}+\frac{V_{f,L}}{\rho_L} \right ) Sf​=m˙δ(ρv​Vf,v​​+ρL​Vf,L​​)

Cahn - Hilliard 方程改为：

∂ ϕ ∂ t + u ⋅ ▽ ϕ = ▽ ⋅ γ λ ε 2 ▽ ψ + S f \frac{\partial \phi }{\partial t}+\textbf{u}\cdot \bigtriangledown \phi =\bigtriangledown \cdot \frac{\gamma \lambda }{\varepsilon ^2}\bigtriangledown \psi+S_f ∂t∂ϕ​+u⋅▽ϕ=▽⋅ε2γλ​▽ψ+Sf​ （6）

对其乘以一个试函数 test(psi) 再对域内积分：

0 = ∫ Ω t e s t ( p s i ) ∗ ( ∂ ϕ ∂ t + u ⋅ ▽ ϕ − ▽ ⋅ γ λ ε 2 ▽ ψ − S f ) d Ω 0=\int_{\Omega}test(psi)*(\frac{\partial \phi }{\partial t}+\textbf{u}\cdot \bigtriangledown \phi -\bigtriangledown \cdot \frac{\gamma \lambda }{\varepsilon ^2}\bigtriangledown \psi-S_f)d\Omega 0=∫Ω​test(psi)∗(∂t∂ϕ​+u⋅▽ϕ−▽⋅ε2γλ​▽ψ−Sf​)dΩ

运用散度定理，再整理一下：

0 = ∫ Ω [ t e s t ( p s i ) ∗ ∂ ϕ ∂ t + ▽ ( t e s t ( p s i ) ) ⋅ γ λ ϵ 2 ▽ ψ ] d Ω + ∫ Ω [ u ⋅ ▽ ϕ − S f ] ∗ t e s t ( p s i ) d Ω − ∫ ∂ Ω t e s t ( p s i ) γ λ ε 2 ▽ ψ ⋅ n d s 0=\int_{\Omega}[test(psi)*\frac{\partial \phi }{\partial t}+\bigtriangledown(test(psi))\cdot\frac{\gamma \lambda}{\epsilon^2}\bigtriangledown\psi]d\Omega+\int_{\Omega}[\textbf{u}\cdot\bigtriangledown\phi-S_f]*test(psi)d\Omega-\int_{\partial\Omega}test(psi)\frac{\gamma \lambda }{\varepsilon ^2}\bigtriangledown \psi\cdot\textbf{n}ds 0=∫Ω​[test(psi)∗∂t∂ϕ​+▽(test(psi))⋅ϵ2γλ​▽ψ]dΩ+∫Ω​[u⋅▽ϕ−Sf​]∗test(psi)dΩ−∫∂Ω​test(psi)ε2γλ​▽ψ⋅nds (7)

(7)式内第一项，对应COMSOL中弱形式：

-pf.lam*pf.mobility*psix*test(psix)/pf.epsilon_pf^2-pf.lam*pf.mobility*psiy*test(psiy)/pf.epsilon_pf^2-phipft*test(psi)

(7)式内第三项，对应于边界条件。

(7)式内第二项，对应于COMSOL中弱形式：

(-pf.ux*phipfx-pf.uy*phipfy)*test(psi)

然而，该弱形式不含源项，需要手动添加。可以在相场节点下，添加弱贡献（weak contribution）， t e s t ( p s i ) ∗ S f test(psi)*S_f test(psi)∗Sf​。

至于内置模块方程用什么试函数，可以在组件里查看。

对于流动的连续性方程，原方程不包含质量源项，应该增加，用上述方法比较容易知道应该增加弱贡献 t e s t ( p ) ∗ S m test(p)*S_m test(p)∗Sm​.

至于相变热量的吸收和释放，因为传热方程中本身带有热源项，可以直接处理为热源，不需要对方程进行改动。

以上就介绍了如何使用弱贡献修改内置方程，具体案例提供一个视频教程。案例文件可在我的主页下载。