Interface Compression Technique in PPM
-
摘要: 高精度多组分分段抛物线法(Piecewise Parabolic Method,PPM)在对可压缩多相流问题进行模拟计算时,在不同组分交界面上存在界面扩散。为此,通过引入包含界面压缩和密度修正的人工界面压缩方法,抑制界面扩散现象。采用一个界面函数表示运动的物质界面,在多组分质量守恒方程和输运方程中添加考虑人工压缩和人工黏性的压缩源项,并在伪时间内采用二阶中心差分法和两步Runge-Kutta方法进行离散求解,采用Strang型分裂格式实现了整体算法的时间二阶精度。一维与二维数值模拟试验表明,结合人工界面压缩之后的PPM能有效抑制界面上数值扩散问题,在长时间的数值模拟中,人工界面压缩能够将扩散界面厚度维持在一定网格之内且保持界面形状不改变,尤其对于涉及稀疏波的问题,如激波引起的水中气泡坍塌,界面压缩效果更为显著。Abstract: This paper describes an artificial interface compression technique for the multi-fluid piecewise parabolic method (PPM). The proposed approach enables the simulation of interfaces between compressible multi-fluid flows with high density ratios and strong shock waves. A compression source term incorporated both interface compression and density correction is added to the mass conservation equation. The compression source term is solved in pseudo-time steps using the interface compression technique and the advection part is solved by multi-fluid PPM. The Strang splitting algorithm achieves second-order accuracy by combining the solutions of the advection operator and the interface compression operator. Numerical tests on the interaction of shock waves with interfaces in compressible multi-fluid flows reveal that multi-fluid PPM combined with the artificial interface compression technique can effectively prevent the smearing phenomenon, which is often observed at the contact interface. For long-time simulations, artificial interface compression with interface sharpening can constrain the thickness of the diffused interface to a few cells and maintain the interface profile. This artificial interface compression technique works well with multi-fluid PPM and the effect is obvious. It is a significant step in the accurate simulation of the collapse of air cavities in water, which involves strong rarefaction waves.
-
磁驱动固体套筒内爆是指电流通过金属套筒表面时,在洛仑兹力的作用下金属套筒径向向内箍缩内爆的物理过程。1973年,Turchi等[1]首次提出磁驱动固体套筒内爆的概念。自20世纪90年代以来,磁驱动固体套筒实验被广泛应用于高压状态方程[2]、材料本构[3]、层裂损伤[4]、磁瑞利-泰勒(Magneto-Rayleigh-Taylor,MRT)不稳定性发展[5–6]、Richtmyer-Meshkov(RM)不稳定性发展[7]等研究。
磁驱动固体套筒实验涉及热扩散、磁扩散、焦耳加热、弹塑性、断裂、层裂等物理过程,并伴有大变形、界面不稳定性等现象。磁驱动固体套筒理论有薄壳模型[8–10]、不可压缩模型[11–13]、电作用量-速度模型[14–15]、全电路模型[15]和磁流体力学模型[16–17]等。这些理论模型已被用于脉冲功率装置、磁驱动固体套筒实验的模拟、设计和研究[7–17]。阚明先等[17]采用二维磁流体力学程序MDSC2模拟回流罩结构磁驱动固体套筒实验时发现,根据回流罩结构磁驱动固体套筒实验测量的电流或回路电流不能直接模拟磁驱动固体套筒,模拟的套筒速度总是比测量速度大,即回路电流并不完全从固体套筒表面流过。回路电流与固体套筒上通过的电流之间存在一个电流系数。由于MDSC2程序[17]以外的理论计算或数值模拟都未提到电流系数,因此,本研究采用其他理论模型对磁驱动固体套筒实验进行模拟,分析回路电流与通过固体套筒的电流之间的关系,通过模拟分析不同回流罩结构固体套筒实验,进一步探讨磁驱动固体套筒实验中电流系数的影响因素和变化规律。
1. 负载结构
大电流脉冲装置上的固体套筒实验通常采用回流罩结构[15, 17–18]。回流罩结构固体套筒实验的初始结构的rz剖面如图1所示,其中,虚线为对称轴。回流罩结构固体套筒实验装置从外到内依次为金属回流罩、绝缘材料和金属套筒,套筒两端为金属电极,上端为阳极,下端为阴极。回路电流从回流罩金属流入,绕过绝缘材料,经过套筒的外表面从阴极流出。电流加载后,电极外面的固体套筒被切割成与阴阳极之间的间隙等高的套筒,在洛仑兹力作用下沿径向向内箍缩。表1为FP-2装置[19]中回流罩结构磁驱动固体套筒实验的套筒参数。图2显示了FP-2装置上不同实验测得的电流变化曲线,电流的上升时间约为
5500 ns,电流峰值为9~11 MA。表 1 磁驱动固体套筒实验的套筒参数Table 1. Liner parameters of the magnetically driven solid liner experimentsExp. No. Liner material Liner’s inner radius/mm Liner’s thickness/mm 1 Al 45 0.6 2 Al 30 0.6 3 Al 45 1.6 4 Al 30 1.9 2. 电流系数的不可压缩模型验证
在薄壳模型、不可压缩模型、电作用量-速度模型、全电路模型、磁流体力学模型等[8–16]适用于磁驱动固体套筒的理论模型中,固体套筒边界的磁感应强度(B)为
B(t)=μ0Iexp(t)2πro (1) 式中:μ0为真空磁导率,Iexp(t)为磁驱动实验测量电流,ro为固体套筒的外半径。
二维磁驱动数值模拟程序MDSC2是由中国工程物理研究院流体物理研究所开发的二维磁流体力学程序[20–21]。该程序已被广泛应用于磁驱动飞片发射、超薄飞片、磁驱动准等熵压缩、磁驱动样品等实验的模拟研究[22–25]。最近,研究人员发现,采用MDSC2程序模拟FP-2装置上的磁驱动固体套筒实验时,基于实验测量的电流或回路电流并不能正确模拟套筒的动力学过程,模拟的套筒速度总是比实验测量值大。为正确模拟FP-2装置上的磁驱动固体套筒实验,需将边界磁感应强度公式[17]修正为
B(t)=μ0fcIexp(t)2πro (2) 式中:fc为回流罩结构rz柱面套筒的电流系数,fc<1。由于文献[17]之外的理论计算或数值模拟中均未提到电流系数fc,因此,需要确定fc是回流罩固体套筒实验固有的,还是MDSC2程序造成的。下面采用固体套筒的不可压缩模型理论确认电流系数是否存在。
在磁驱动固体套筒的不可压缩模型[11–13]中,不考虑套筒的磁扩散,假设磁压只作用于套筒的外表面,且磁压做功全部转化为套筒动能,套筒不可压缩,只作径向运动。设ρ为套筒密度,h为套筒高度,vo为套筒外界面速度,ri、vi分别为套筒内半径和内界面速度,r、v为套筒内某点的径向位置(ri≤r≤ro)和速度,由不可压缩假设,有
rivi=rovo (3) rv=rovo (4) 则套筒总动能Ek为
Ek=∫roriρπrhv2dr=πρhr2ov2olnrori (5) 由于磁压只作用于套筒的外表面,且磁压做功全部转化为套筒动能,则
dEkdt=2πμ0rohvoB2 (6) 将式(5)代入式(6)并积分,可得
dvodt=−v2oro−1ln(ro/ri)[B22μ0ρro+v2o2ro(1−r2or2i)] (7) dvidt=−v2iri−1ln(ro/ri)[B22μ0ρri+v2o2ri(1−r2ir2o)] (8) 采用上述不可压缩模型,对固体套筒实验4进行不可压缩模型模拟验证。图3给出了采用不可压缩模型模拟得到的套筒内界面速度。显然,采用回路电流或测量电流直接模拟的套筒速度明显比实验测量速度大,后者是前者的0.82倍,即计算不可压缩模型的边界磁感应强度时不能用式(1),而是用式(2)。不可压缩模型的模拟结果表明,对于回流罩固体套筒实验,回路电流或测量电流与固体套筒上通过的电流之间的电流系数不是MDSC2程序造成的,而是回流罩固体套筒实验固有的。
3. 电流系数规律
从第2节的模拟可知,磁驱动固体套筒理论的边界磁感应强度公式中包含电流系数,它反映了有多少回路电流从套筒实际流过。在磁驱动实验中,实验测量的电流是流入回流罩之前的电流,即回路电流,而不是从套筒直接流过的电流。从套筒流过的电流很难被直接测量,因此,电流系数难以预知。回流罩的结构比较复杂,阴阳电极之间连有金属套筒、绝缘材料,金属套筒与绝缘材料之间是真空,回流罩结构的分流机制包括阴阳极间的并联电路分流、漏磁、真空击穿等。事实上,电流系数是通过数值模拟发现的,由磁流体力学程序模拟速度与磁驱动套筒实验测量速度的对比确定。当前的固体套筒实验的模拟都是后验的,无法直接正确预测,因此,研究电流系数的变化规律非常重要,是正确设计和预测固体套筒实验的基础。
由于磁流体力学模型[21, 26]是包含固体弹塑性、热扩散、磁扩散等物理过程的可压缩模型,能够比不可压缩模型更加准确地描述磁驱动固体套筒实验,因此,下面将采用MDSC2程序对FP-2装置上开展的磁驱动固体套筒实验的电流系数变化规律进行研究。
图4给出了实验1~实验4的套筒内界面模拟速度。可以看出,应用式(2)的磁流体力学模型能正确描述磁驱动固体套筒实验。然而,不同的磁驱动固体套筒实验对应的电流系数是不同的。回流罩结构磁驱动固体套筒实验的电流系数和套筒的初始尺寸列于表2。
表 2 磁驱动固体套筒实验的电流系数Table 2. Current coefficients of the magnetically driven solid liner experimentsExp. No. Liner’s inner radius/mm Liner’s thickness/mm fc 1 45 0.6 0.87 2 30 0.6 0.90 3 45 1.6 0.85 4 30 1.9 0.88 由表2可知:电流系数是常数,不随时间的发展而变化,即电流系数与实验过程无关;对于不同的套筒,电流系数有所不同,说明电流系数与套筒的初始结构有关。由实验1和实验2可知,当套筒厚度相同时,若套筒内半径不同,则电流系数不同,且内半径越大,电流系数越小。对比实验1和实验3,或者实验2和实验4可知,当套筒内半径相同时,若套筒厚度不同,则电流系数不同,且套筒厚度越大,电流系数越小。
4. 结 论
采用不可压缩模型验证了回流罩结构磁驱动固体套筒实验中电流系数的存在,即回流罩结构磁驱动固体套筒实验的实验电流/回路电流并不完全从负载套筒的表面通过,实验电流/回路电流与套筒表面流过的电流之间存在一个电流系数。采用包含固体弹塑性、热扩散、磁扩散的磁流体力学模型,对回流罩结构磁驱动固体套筒实验的电流系数进行了确定和分析,结果显示,磁流体力学模型和有电流系数的边界磁感应强度公式能正确模拟回流罩结构磁驱动固体套筒实验。电流系数与套筒结构的关系为:
(1) 不同套筒对应的电流系数不同;
(2) 电流系数与实验过程无关,由套筒初始结构决定;
(3) 套筒厚度相同时,电流系数由套筒内半径决定,套筒内半径越大,电流系数越小;
(4) 套筒内半径相同时,电流系数由套筒厚度决定,套筒厚度越大,电流系数越小。
正确认识磁驱动固体套筒实验的电流系数变化规律,使磁驱动固体套筒实验的磁流体模拟从后验模拟发展成先验的准确设计和预测,有助于降低实验成本,加快柱面相关的实验研究。
-
表 1 水中气泡塌陷问题中状态方程参数及初始参数
Table 1. Equation of state parameters and the initial time data of an air cavity collapse in water
Material ρ/(kg·m–3) p/Pa u/(m·s–1) v/(m·s–1) γ Water (Post-shock) 1.325 1.915×104 68.52 0 4.4 Water (Pre-shock) 1 1 0 0 4.4 Air 0.001 1 0 0 1.4 表 2 空气-R22气柱相互作用问题中的状态方程参数及初始参数
Table 2. Equation of state parameters and the initial time data of air-R22 shock-cylinder interaction
Material ρ/(kg·m–3) p/MPa u/(m·s–1) v/(m·s–1) γ Air (Post-shock) 1.686 0.159 –113.5 0 1.400 Air (Pre-shock) 1.225 0.101 0 0 1.400 R22 3.863 0.101 0 0 1.249 -
[1] MESHKOV E E. Instability of the interface of two gases accelerated by a shock wave [J]. Fluid Dynamics, 1969, 4(5): 101–104. [2] BROUILLETTE M. The Richtmyer-Meshkov Instability [J]. Annual Review of Fluid Mechanics, 2002, 34(34): 445–468. [3] LOMBARDINI M, PULLIN D I, MEIRON D I. Turbulent mixing driven by spherical implosions (Part 1): flow description and mixing-layer growth [J]. Journal of Fluid Mechanics, 2014, 748(2): 85–112. [4] LOMBARDINI M, PULLIN D I, MEIRON D I. Turbulent mixing driven by spherical implosions (Part 2): turbulence statistics [J]. Journal of Fluid Mechanics, 2014, 748(10): 113–142. [5] CLEMENS N T, MUNGAL M G. Large-scale structure and entrainment in the supersonic mixing layer [J]. Journal of Fluid Mechanics, 1995, 284(284): 171–216. [6] KAWAI S, LELE S K. Large-eddy simulation of jet mixing in supersonic crossflows [J]. American Institute of Aeronautics and Astronautics, 2010, 48(9): 2063–2083. doi: 10.2514/1.J050282 [7] JOHNSEN E, COLONIUS T. Shock-induced collapse of a gas bubble in shockwave lithotripsy [J]. The Journal of the Acoustical Society of America, 2008, 124(4): 2011–2020. doi: 10.1121/1.2973229 [8] KLASEBOER E, HUNG K C, WANG C, et al. Experimental and numerical investigation of the dynamics of an underwater explosion bubble near a resilient/rigid structure [J]. Journal of Fluid Mechanics, 2005, 537(537): 387–413. [9] RANJAN D, OAKLEY J, BONAZZA R. 3D shock-bubble interactions [J]. Physics of Fluids, 2013, 25(9): 117–140. [10] THEOFANOUS T G. Aerobreakup of newtonian and viscoelastic liquids [J]. Annual Review of Fluid Mechanics, 2011, 43(1): 661–690. doi: 10.1146/annurev-fluid-122109-160638 [11] COLELLA P, WOODWARD P R. The piecewise parabolic method (PPM) for gas-dynamical simulations [J]. Journal of Computational Physics, 1984, 54(1): 174–201. doi: 10.1016/0021-9991(84)90143-8 [12] 马东军, 孙德军, 尹协远. 高密度比多介质可压缩流动的PPM方法 [J]. 计算物理, 2001, 18(6): 517–522. doi: 10.3969/j.issn.1001-246X.2001.06.008MA D J, SUN D J, YIN X Y. Piecewise parabolic method for compressible flows of multifluids with high density ratios [J]. Chinese Journal of Computational Physics, 2001, 18(6): 517–522. doi: 10.3969/j.issn.1001-246X.2001.06.008 [13] BAI J S, WANG B, WANG T, et al. Numerical simulation of the Richtmyer-Meshkov instability in initially nonuniform flows and mixing with reshock [J]. Physical Review E, 2012, 86(6): 066319. doi: 10.1103/PhysRevE.86.066319 [14] BAI J S, ZOU L Y, WANG T, et al. Experimental and numerical study of shock-accelerated elliptic heavy gas cylinders [J]. Physical Review E, 2011, 82(2): 056318. [15] XIAO J X, BAI J S, WANG T. Numerical study of initial perturbation effects on Richtmyer-Meshkov instability in nonuniform flows [J]. Physical Review E, 2016, 94(1): 013112. doi: 10.1103/PhysRevE.94.013112 [16] SHYUE K M. A fluid-mixture type algorithm for compressible multicomponent flow with van der Waals equation of state [J]. Journal of Computational Physics, 1999, 156(1): 43–88. doi: 10.1006/jcph.1999.6349 [17] SHYUE K M. An efficient shock-capturing algorithm for compressible multicomponent problems [J]. Journal of Computational Physics, 1998, 142(1): 208–242. doi: 10.1006/jcph.1998.5930 [18] ALLAIRE G, CLERC S, KOKH S. A five-equation model for the simulation of interfaces between compressible fluids [J]. Journal of Computational Physics, 2002, 181(2): 577–616. doi: 10.1006/jcph.2002.7143 [19] JOHNSEN E, COLONIUS T. Implementation of WENO schemes in compressible multicomponent flow problems [J]. Journal of Computational Physics, 2006, 219(2): 715–732. doi: 10.1016/j.jcp.2006.04.018 [20] KOKH S, ALLAIRE G. Numerical simulation of 2-D two-phase flows with interface [C]//TORO E F. Godunov Methods.Boston, MA: Springer, 2001: 513–518. [21] II S, XIE B, XIAO F. An interface capturing method with a continuous function: the THINC method on unstructured triangular and tetrahedral meshes [J]. Journal of Computational Physics, 2014, 259: 260–269. [22] SHYUE K M, XIAO F. An Eulerian interface sharpening algorithm for compressible two-phase flow: the algebraic THINC approach [J]. Journal of Computational Physics, 2014, 268(2): 326–354. [23] XIAO F, HONMA Y, KONO T. A simple algebraic interface capturing scheme using hyperbolic tangent function [J]. International Journal for Numerical Methods in Fluids, 2005, 48(9): 1023–1040. doi: 10.1002/fld.975 [24] XIAO F, II S, CHEN C. Revisit to the THINC scheme: a simple algebraic VOF algorithm [J]. Journal of Computational Physics, 2011, 230(19): 7086–7092. doi: 10.1016/j.jcp.2011.06.012 [25] KOKH S, LAGOUTIÈRE F. An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model [J]. Journal of Computational Physics, 2010, 229(8): 2773–2809. doi: 10.1016/j.jcp.2009.12.003 [26] FRIESS M B, KOKH S. Simulation of sharp interface multi-material flows involving an arbitrary number of components through an extended five-equation model [J]. Journal of Computational Physics, 2014, 273(273): 488–519. [27] DELIGANT M, SPECKLIN M, KHELLADI S. A naturally anti-diffusive compressible two phases Kapila model with boundedness preservation coupled to a high order finite volume solver [J]. Computers and Fluids, 2015, 114(1): 265–273. [28] OLSSON E, KREISS G, ZAHEDI S. A conservative level set method for two phase flow II [J]. Journal of Computational Physics, 2007, 225(1): 785–807. doi: 10.1016/j.jcp.2006.12.027 [29] SHUKLA R K, PANTANO C, FREUND J B. An interface capturing method for the simulation of multi-phase compressible flows [J]. Journal of Computational Physics, 2010, 229(19): 7411–7439. doi: 10.1016/j.jcp.2010.06.025 [30] SHUKLA R K. Nonlinear preconditioning for efficient and accurate interface capturing in simulation of multicomponent compressible flows [J]. Journal of Computational Physics, 2014, 276(1): 508–540. [31] FREUND J B, SHUKLA R K, EVAN A P. Shock-induced bubble jetting into a viscous fluid with application to tissue injury in shock-wave lithotripsy [J]. The Journal of the Acoustical Society of America, 2009, 126(5): 2746–2756. doi: 10.1121/1.3224830 [32] SHU C W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws [C]//QUARTERONI A. Advanced Numerical Approximation of Nonlinear Hyperbolic Equations: Lecture Notes in Mathematics (Vol.1697). Berlin, Heidelberg: Springer, 1998: 325–432. [33] NOURGALIEV R R, DINH T N, THEOFANOUS T G. Adaptive characteristics-based matching for compressible multifluid dynamics [J]. Journal of Computational Physics, 2006, 213(2): 500–529. doi: 10.1016/j.jcp.2005.08.028 [34] HU X Y, KHOO B C. An interface interaction method for compressible multifluids [J]. Journal of Computational Physics, 2004, 198(1): 35–64. doi: 10.1016/j.jcp.2003.12.018 [35] SHYUE K M. A wave-propagation based volume tracking method for compressible multicomponent flow in two space dimensions [J]. Journal of Computational Physics, 2006, 215(1): 219–244. doi: 10.1016/j.jcp.2005.10.030 [36] HAAS J F, STURTEVANT B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities [J]. Journal of Fluid Mechanics, 1987, 181(1): 41–76. -