基于MOF界面重构的多物质ALE方法

［ 8 － 11 ］ . 该模型用声波 Riemann 解估计界面速度， 分数保持不变. ⑥一维子网格力学模型 并由此计算每种物质 的体积分数的改变量. 利用总能量守恒原理计算混合网格中各种物质的共同的压力 ， 该压力用在动量方程
中. 每种物质的内能用它自己的内能方程进行更新 . 此外， 还考虑了混合网格中各种物质之间内能的交换 ， 即 相互做功. 界面重构的多物质 ALE 方法. 空间离散采用显式相容有限元方法 . 对混合网 ［ 13 － 14 ］ 格提出了一种新的二维子网格力学模型 ， 用近似 Riemann 解 HLLC 计算界面的法向速度， 界面切向速 度用混合网格的节点速度插值计算 . 提出了一种精确积分守恒重映方法 ， 该方法用计算机图形学中的多边形 ［ 15 ］ 剪裁算法 计算老网格和新网格的交点. 重映后采用 MOF 方法重构界面. 给出了若干数值算例.
收稿日期： 2008 － 03 － 23 ； 修回日期： 2009 － 09 － 15 10671025 ） 、 基金项目： 国 家 973 项 目 （ 2005CB321700 ） 、 国 家 自 然 科 学 基 金 （ 10771019 ， 中国工程物理研究院科学技术基金 （ 2008A0201009 ） 和计算物理重点实验室基金（ 9140C690101070C69 ） 资助项目 作者简介： 贾祖朋（ 1967 － ） ， 湖南慈利， 副研究员， 从事计算数学和计算流体力学的研究， 北京市海淀区丰豪东路 2 号 5 室 100094.
0
引言
［ 1］ 很多实际问题的数值模拟会遇到界面大变形和网格大变形问题 。 对界面大变形问题， 有 level set 、 ［ 3］ VOF［2］等界面追踪方法；对网格大变形问题， 有 ALE 方法 . 多物质大变形问题的数值模拟需要同时处理界
面大变形和网格大变形. 传统的拉氏方法和 ALE 方法无法解决这个难题. 多物质 ALE 方法是解决该类问题 的比较有效的方法. 该方法有以下几个特点：①允许一个网格内含有多种物质， 物质界面可以穿过网格；②用 VOF 方法等界面追踪方法处理界面的大变形 . 重分重映处理网格的大变形；③用 level set 方法、 多物质 ALE 方法的关键之一是对混合网格提出一个封闭模型 ， 用以计算经过一个拉氏步后混合网格中 ［ 4 － 5］ . 在这个模 每种物质的体积分数和物理量的变化 . 文献中的方法大概分为六种： ① 平均应变率混合模型 型中， 混合网格中的每种物质的应变率取为网格的平均应变率 . 当混合网格中界面两侧密度比或压力比很大 ［ 4 － 6］ . ②压力平衡混合模型. 如果计算网格的尺寸足够小， 应用这个模型会得到非物理的结果 使得混合网 时， 格中的压力达到平衡的时间尺度小于计算所用的时间步长 ， 在单元水平上假设压力达到瞬间平衡是物理过 ［ 7］ 程的一个合理的近似. 压力平衡是非线性问题， 需要设计复杂的迭代格式进行计算 . ③ 压力松弛模型. 这
（ 12 ） （ 13 ）
更新节点的位置以后， 就可以计算混合网格 z 在 t ρz 与 2. 3 节中类似， 可以分 5 步计算 t 时刻的坐标. 然后用几何工具箱计算 t 物质的体积 V 混合网格 t 令 Qp z = －
+1 / 2 f zr ( mp )·v n ， ∑ r r n +1 n +1 z， k n +1 n +1 / 2 n +1
时刻的量
n 2
首先计算 t 时刻混合网格的压力， 它定义为混合网格中两种物质的压力的体积加权平均 ， 即 pn z =
n r n n α z， ∑ k p z， k， k =1
（8）
其中 α z， k 是体积分数 . 然后计算角力 f z 和节点力 f r . t 关于节点 r， 其加速度、 v
n +1
* } S L = min { v nL － a L ， ， v* n － a * } S R = max { v nR + a R ， . v* n + a * * n
图1 Fig. 1
混合网格示意图
Illustration of mixed cells
（5）
中间波的波速为
第3 期
贾祖朋：基于 MOF 界面重构的多物质 ALE 方法
［ 8］ 种模型引进一种类似于粘性的松弛机制使混合网格中的压力趋于平衡 ， 有的采用等熵假设 . ④ 接触混合 ［ 5］ 模型 . 通过求解应力和应变率在界面上的跳跃条件 ， 允许界面两侧滑移和分离. ⑤ 体积模量加权平均混合
模型
［ 6］
. 在这个模型中， 当混合网格被压缩时， 采用体积模量加权平均， 当混合网格膨胀时， 每种物质的体积
355
SM =
ρ R v nR ( S R － v nR ) － ρ L v nL ( S L － v nL ) + p L － p R ， ρ R ( S R － v nR ) － ρ L ( S L － v nL )
（6）
交点 ( xcs ( 1 ) ， ycs ( 1 ) ) 的法向速度就取为 S M ； 5 ） 用前面两步中求出的两个交点 t n 时刻的切向速度和法向速度计算这两个交点 t n 时刻的速度. 求出两个交点 t 时刻的速度以后， 就可以计算这两个点在 t t
r z
（1）
用单元形函数定义角质量和单元 z 在节点 r 处的角力 f . 利用角质量定义单元质量 M z 和节点质量 M r . 与文 ［ 26 ］ 一样， 用边人工粘性消除激波振荡， 用子网格扰动压力抑制网格的非物理变形 . 单元 z 的密度用 （2） ρz = Mz / Vz ( t ) 计算， 其中 V z ( t ) 表示 t 时刻单元 z 的面积. 节点 r 的全离散动量方程为 M r Δv r = ΔtF r = Δt∑ f rz ，
z
（3）
单元 z 的内能的改变量由下式计算： Δe z =
+1 / 2 f zr ·v n Δt / M z . ∑ r r
（4）
2
2. 1
二维子网格力学模型
混合网格的描述
假设混合网格内有两种物质， 界面是穿过网格的一条线段. 考虑如图 1 所示的一个混合网格， 其 4 个节 点按逆时针方向 排 序. 对 4 个 节 点 定 义 物 质 状 态 变 量 ns，如 果 节 点 k ( k = 1 ， 则 2， 3， 4 ) 在 物 质 1 里 面， ns ( k ) = 1 ， ics ( 2 ) = 4 . 相应的两个 否则 ns ( k ) = 2 . 变量 ics 表示与直线界面相交的边， 在图 1 中 ics ( 1 ) = 2 ， 交点为 ( xcs ( 1 ) ， 在与之对应的边上属于物质 1 和物质 2 的节点 ycs ( 1 ) ) 和 ( xcs ( 2 ) ， ycs ( 2 ) ) . 对每个交点， ls2 ( 1 ) = 3 ， ls1 ( 2 ) = 1 ， ls2 ( 2 ) = 4 . 由此可用 2. 2 节中的几何 分别用 ls1 和 ls2 表示， 在图 1 中， 有 ls1 ( 1 ) = 2 ， 工具箱确定由直线界面截取的两个多边形， 计算这两个 ( ( ) ， 多边形的重心 即两个物质中心点 xce 1 ， yce ( 1 ) ) 和
n +1 / 2 n n +1 / 2
时刻的坐标. 同样可以计算 4 个节点在
n +1 / 2
时刻的坐标. 然后用几何工具箱计算 t p
n +1 / 2 i n i n i
n +1 / 2
时刻界面与混合网格的边的交点 ， 再计算 t
n +1 / 2 n i
时刻混合网格 （7）
中两种物质的体积分数和密度. 利用等熵假设， 可以计算 t = p r .
n r
n +1 / 2
t 时刻的速度、
n +1
时刻的坐标分别为 （9） （ 10 ） （ 11 ）
ar = fr / Mr .
n +1 r +1 / 2 n +1 ) vn = ( vn /2. r r + vr +1 n +1 / 2 ( 1)， xn = xn r r + Δ tv r +1 n +1 / 2 ( 2). yn = yn r r + Δ tv r n +1
（ 北京应用物理与计算数学研究所，北京 100088 ）
摘
ofFluid） 界面重构的多物质 ALE （ Arbitrary LagrangianEulerian） 方法. 流 要： 提出一种基于 MOF （ Moment-
体力学方程组采用相容有限元方法进行空间离散 . 提出一种新的二维子网格力学模型，用来计算混合网格中的 物理量经过一个拉氏步后发生的变化，混合网格内的界面重构采用 MOF 方法. 提出一种精确积分守恒重映方 法. 给出数值算例，如空气和水的 Riemann 问题， Dukowicz 问题，水中强激波与空气泡相互作用问题等 . 结果 表明，方法具有较高的精度，能够处理物质界面和网格的大变形问题 . 关键词： 多物质 ALE 方法； 相容有限元方法； MOF 方法 中图分类号： O241. 82 文献标识码： A
本文提出一种基于 MOF
［ 12 ］
354
计
算
物
理
第 27 卷
1
显式相容有限元方法
26 ］ 文［ 提出了三维笛卡尔坐标系中流体力学的显式相容有限元拉氏方法 . 在二维笛卡尔坐标系中可以 建立类似的方法. 采用交错网格， 即压力、 密度和内能定义在单元中心， 速度和加速度定义在网格节点上. 采 用四边形单元对计算区域进行离散 . 参考坐标系 ( ξ， η ) 中的标准单元及相应的形函数分别为 { ( ξ， k = 1， …， 4， ψ k ( ξ， η ) － 1 ≤ ξ， η ≤ 1 }， η ) = ( 1 ± ξ ) ( 1 ± η ) /4，
n ( i
时刻混合网格中两种物质的压力 i = 1， 2.
c
n )2 i
δV
n +1 / 2 i
/V ，
其中 c 是声速， δV 权平均. 2. 4 计算 t
n +1
n i
n +1 / 2 i
=V
n +1 / 2 i
－ V . 混合网格 t
n +1 / 2
时刻的压力定义为混合网格中两种物质的压力的体积加
考虑如图个节点定义物质状态变量ns如果节点里面则ns1否则ns变量ics表示与直线界面相交的边在图相应的两个交点为的节点分别用ls1ls2表示在图节中的几何工具箱确定由直线界面截取的两个多边形计算这两个多边形的重心即两个物质中心点几何工具箱给定混合网格中直线界面的方程和个节点的坐标可以用几何工具箱确定节中所定义的变量nsicsxcsycsls1ls2xceyce确定节点位置
( xce ( 2 ) ， 以及两种物质的体积分数. yce ( 2 ) ) ，
2. 2
几何工具箱
给定混合网格中直线界面的方程和 4 个节点的坐 标， 可以用几何工具箱确定 2. 1 节中所定义的变量 ns， ics， xcs， ycs， ls1 ， ls2 ， xce， yce 的值. 几何工 具 箱 包 含 4 个 几何工具： ① 确定节点位置； ② 求直线与直线的交点； ③截取多边形；④求多边形的面积和重心. 2. 3 计算 t n + 1 / 2 时刻的量 首先计算 t 时刻的界面速度. 分以下 5 步： 1 ） 计算 t n 时刻界面的法向单位矢量 n. 假设其方向是由物质 1 指向物质 2 ； 2 ） 计算 t n 时刻 4 个节点的法向速度和切向速度； 3 ） 用 t n 时刻两个交点各自对应的两个网格节点的切向速度通过线性插值计算这两个交点的切向速度 ； 4 ） 通过用近似 Riemann 解 HLLC ， 求解两个法向 Riemann 问题， 计算两个交点的法向速度. v nR 分别表示节点 ls1 ( 1 ) 和 ls2 ( 1 ) 的法向速度， pL ， a L 分别 例如， 考虑交点 ( xcs ( 1 ) ， ρL ， ycs ( 1 ) ) . 令 v nL ， pR ， a R 分别表示混合网格中物质 2 的密度、 vn 表示混合网格中物质 1 的密度、 压力和声速， ρR ， 压力和声速， 和 a 分别表示左右法向速度和声速的 Roe 平均值. 左波和右波的波速估计为
n +1 n +1 时刻的体积 V z . 混合网格 z 在 t 时刻的密度为
第 27 卷 第 3 期 2010 年 5 月
计
算
物
理
CHINESE JOURNAL OF COMPUTATIONAL PHYSICS
Vol. 27 ， No. 3 May， 2010
246X（ 2010 ） 03-0353-08 文章编号：1001-
基于 MOF 界面重构的多物质 ALE 方法
贾祖朋