平滑粒子流体动力学的基本方程

平滑粒子流体动力学的基本方程

0sph方法的发展在研究中，我们发现了带有不同介质界面的三维流场三维系统的复杂变形过程的数值模拟问题。这不仅可以克服传统拉网法和欧洲网格法的缺点，而且可以比pic法更简单。经过调研,我们决定研究、开发SPH(SmoothedParticleHydrodynamics)方法。我们认为,作为PIC方法进一步演变而出现的若干新方向中的一个,SPH方法是比较具有新意的。它在1977年以后才发展起来,是一种纯拉格朗日方法,本质上不需要借助于网格。此方法面世以来一直在天体物理领域里成功地应用。1983年Monaghan等人提出了适用于SPH方法的人为粘性项,从而能模拟流场中的激波强间断现象。1990年Monaghan又提出了适用于SPH方法的人为热流项借以减轻由于引入人为粘性项后所产生的某些误差。近年来,SPH方法已扩展应用于高速碰撞等材料动载响应诸方面的数值模拟。SPH方法通过一个称为“核函数”的积分核进行“核函数估值”(KernalEstimation)近似,将流体力学基本方程组转换成计算用的公式。整个流场离散成一系列“粒子”,所有力学量由这些粒子负载。这些粒子可以按计算公式,亦即按流体力学流动的规律相当任意地流动,算法本身对其应用的限制大为减少。SPH方法一开始就是针对三维空间的。当然,也有低于三维的研究工作,例如,二维轴对称空间的SPH方法。所有这些特点,使SPH方法有可能满足惯性结束聚变中对界面不稳定性数值模拟的要求。本文内容包括方法的基本概念、基本方程的推导、计算方法、程序的运行环境和对程序的检验性试算等,并提供了一个球形弹丸高速斜碰撞靶板算例的结果。最后对尚存在的问题作了讨论。1sph的数值计算在SPH方法中,流体介质被分为N个小块,取各小块的质心{rj(t),j=1,\:,N}作为数值模拟用的N个插值基点。把小块j的质量mj赋与质心rj(t),j=1,\:,N,就可以视这些插值基点为质点而称之为粒子。粒子j的速度当然就是drj(t)dt≡uj(t),ddt是随体导数。这些粒子的数目分布密度是n(r,t)=Ν∑j=1δ(r-rj(t)),它们的质量分布密度是ρ(N)(r,t)=Ν∑j=1mjδ(r-rj(t))。当N增大时,ρ(N)趋近于流体介质的密度ρ(r,t)。可见,SPH用的是固着在流体介质上随之流动的活动插值基点。就这一点而言,与拉氏网格差分法中的网格点并没有什么不同。然而SPH不是通过含空间差分的差分方程,而是通过描述流动插值基点(“粒子”)位置及该处流场物理量对时间变化规律的一组常微分方程来实现数值模拟的。这些不出现任何流场物理量空间导数的常微分方程,是以流场物理量的积分估值关系式(“核函数估值”)为基础并利用所选积分核(“核函数”)的有利性质,由流体力学偏微分方程近似导出的。记流场函数f(r)在r处的核函数估值为〈f(r)〉,它的定义是〈f(r)〉=∫DW(r,r′,h)f(r′)dr′.(1)其中W(r,r′,h)是选用的核函数,D是积分域。显然,〈f(r)〉表示r附近各点处f值的加权平均。为了SPH实现数值模拟的需要,通常限定所选用的W满足下列要求:1°.W(r,r′,h)足够地光滑并满足归一化条件∫DW(r,r′,h)dr′=1.2°.量纲为长度的参量h用于指示使W显著不为零的r′(r固定)的取值域线度。当h→0时,W(r,r′,h)→δ(r-r′)。例如,当|r-r′|>2h,常选W是处处为零的。3°.W(r‚r′,h)=W(|r-r′|,h),因而有W(r,r′,h)=W(r′,r,h)和Δ´W=-ΔW=-(r-r′|r-r′|)dWdR(R≡|r-r′|,Δ≡∂∂r,Δ´≡∂∂r′)。为了简便,当r、r′取{rk(t),k=1,\:,N}时,还进一步约定fi≡〈f(ri)〉‚Wij≡W(ri,rj,h)‚ΔiWij≡∂Wij∂ri‚Wij,β≡∂Wij∂xβi‚xβi≡ri的直角坐标β分量以及uβi≡ui的β分量,而β=1,2,3。对(1)应用积分中值定理,由所选W的性质,不难看出〈f(r)〉=f(ˉr)→h→0f(r).(1′)SPH在数值模拟中一般总是用〈f(r)〉去逼近f(r),相当于用离r不远的ˉr处的流场函数值去近似r处的值。按约定的符号,在插值基点ri处,就是用fi去近似f(ri),或说为“用fi去近似粒子i携带的物理量f的值”。另外,积分域D(以r为中心)一般都取得足够大,使得W(r,r′)当r′位于D的边界面上时取值为零,从而保证在分部积分过程中出现的面积分项可以忽略。本文以下的讨论限于这种通常情形。在介质性质或物理量分布变化陡峭的地方,D的取法往往受到限制,无法保证面积分项为零。对于这类特殊情形,需要作更仔细的处理。对于前面提到的n(r,t)和ρ(N)(r,t),由(1)并采用约定符号,可立即得到它们在粒子i所在处的核函数估值为ni=Ν∑j=1Wij‚(2)ρ(Ν)i=Ν∑j=1mjWij.(3)显然,n-1/可用来估计粒子i附近粒子间距的大小,而ρ(Ν)i在h够小、N够大时会较好地逼近ri(t)处流体介质密度ρ(r,(t))的核函数估值ρi[因而也逼近ρ(r)本身]。由(3)看出,只要知道ri(t)、rj(t),就可确定ρ(N)i(mj和W为已知的)。所以,SPH用近似公式ρi=Ν∑j=1mjWij‚i=1,\:,Ν(4)或在通常情形与之等效的常微分方程来确定粒子i所携带的介质密度值(供插值使用)。粒子i所携带的介质压强pi随着ρi也将被状态方程及后面提到的能量方程所确定。为了导出用来确定粒子i的速度ui(t)和它携带的单位质量介质内能εi(=〈ε(ri)〉)的方程,需将(1)写成明显依赖于N个粒子(插值基点)的位置和它所携带的流场函数值(约定用对应的核函数估值)的形式〈f(r)〉=Ν∑j=1W(r,rj,h)fjΔˉVj.(1″)这可以通过将(1)中的积分写成N个对小块积分之和并使用积分中值定理而得到。ΔˉVj将偏离小块j的真实体积ΔVj,以补偿将W、f中本应采用的ˉrj改回rj和用fi代替f(rj(t))所带来的误差。同样道理,在允许作分部积分且略去面积分项从而把Δ´作用到W上的通常情形,就能将f=Δg的核函数估值写成〈Δg(r)〉=Ν∑j=1ΔW(r,rj,h)gjΔˉˉVj‚(1‴)其中的ΔˉˉVj一般略微不同于ΔˉV。进一步,我们需要采用下列近似等式ΔˉˉVj≈ΔˉVj≈mjρj(5)使(1″)、(1ue087)变成只需知道插值基点上的rj,fj即可算出结果的插值公式〈f(r)〉=Ν∑j=1W(r,rj,h)fjρjmj(6)和〈Δf(r)〉Ν∑j=1ΔW(r,rj,h)fjρjmj‚(7)其中fj,ρj都是因rj=rj(t)而随时间变化的,ρj由(4)或(4′)决定。当然(6)、(7)都是近似的,但当N取得够大,误差会足够地小。在通常情形下利用它们不难处理f为几个场函数的乘积并含有空间微分算子的情形。例如一种可能的方案是这样也能将Δ转到W身上去(不过,有的要先借助中值定理把乘积中不含Δ部分移到积分号外去并利用h够小的W会使ˉrj≈ri这一事实)。现在,我们可以由流体力学方程dudt=-1ρΔp想到,SPH中的粒子(插值基点)的运动方程应取为duidt=-〈1ρΔp〉i.不过,通常先作分解1ρΔp=Δ(pρ)+pρ2(Δρ)(9)再利用(8),得到决定粒子速度的常微分方程duidt=-Ν∑j=1mj(piρ2i+pjρ2j)ΔiWij,i=1,\:,Ν.(10)当两边乘以mi后,(10)右边第j项是粒子j对粒子i的作用力。由W的性质知,此力沿二粒子连线方向,交换i与j只改变力的方向不改变力的大小,是一种二体有心力。所以,(10)必能严格保证此N粒子系统的动量及角动量时时守恒。而这些优良性质却与分解(9)密切相关。不作这个分解而直接对1ρΔp使用(8),将得到不同于(10)的方程并且不保证动量、角动量守恒。这里的“毛病”就出在(6)、(7)、(8)都只是近似成立。仿上面,还可以得到决定单位质量内能εi的常微分方程不过,这两个方程都能与运动方程(10)相容,严格保持系统内能与动能之和守恒。鉴于在h小N大的情况下,ΔiWij的性质会使只有i附近的少数粒子参与求和,以致pρ2在求和号之内还是之外通常并不会引起太大的差别。总结以上所述,可以看到:SPH选用的活动插值基点{ri(t),i=1,\:,N}可以视之为携带着物理量mi,ρi,pi,εi并以速度ui(t)=dri(t)dt运动着的N个质点(粒子);决定这些质点状态(ui及各被携带物理量的值)随时间变化规律的是常微分方程(4′)、(10)、(11)并需辅以流体介质状态方程;定解的初始条件是令t=0时各质点的位置、速度及其携带的各物理量等于流体介质上对应小块质心的位置、速度及该处介质的有关物理量的实际值;通过求解如此定解的常微分方程组可以同时求得任意t时刻插值基点的位置和各该点上的各有关物理量(包括u)的值,从而完成数值模拟过程。补充一点。ni=Ν∑j=1Wij的倒数含有粒子i附近每个粒子平均占有之体积的意义。用它来逼近介质小块i的体积ΔVi是挺自然的。因此也可以用ΔˉˉVj≈ΔˉVj≈1ni,i=1,\:,Ν(12)去替代(5)。这将使插值公式(6)取如下形式〈f(r)〉=Ν∑j=1Wijfj(Ν∑k=1Wjk).(13)为了保持质量守恒,还需同时用去替代(4)、(4′)。其它方程[(10)、(11)等]的形式不必改动。显然,当mi=m(i=1,\:,N)时,这两种方案是完全一致的。当插值基点分布很接近介质密度分布,即mj≈mi时这两个方案也不会有明显的差别。核函数的具体形式不唯一,常用高斯分布型函数。我们采用下面的简化形式W(ν,h)=A{1-3ν2/2+3ν3/4,0≤ν≤1(2-ν)3/4,1≤ν≤20,ν≥2(15)式中ν=|ri-rj|/h,A=1/(πh3)是归一化常数。与差分方法一样,SPH方法模拟含有激波的流场时,也需要引入人为粘性项以平滑间断;同时引入人为热流以减轻人为粘性误差。Monaghan和Gingold导出适用于SPH的人为粘性项其中式中Ci,Cj是声速;α,β和ε是可调参数。Mognahan还导出了适用于SPH方法的人为热流项其中式中g1、g2是可调参数。计算时,当ρ˙i≤0时,置Hi=0。计及人为粘性πij和人为热流Hi后,(10)、(11)二式应改写成duiβdt=-∑jmj(Ρiρi2+Ρjρj2+πij)Wij,β‚(25)dεjdt=∑jmj(uiβ-ujβ)(Ρiρi2+12πij)Wij,β+Ηi.(26)不难看出,此两式相当于(10)、(11)二式在其每个含P项加上对于交换i\,j为对称的12πij,所以仍能保持动量、角动量与能量守恒。2计算方法和策略2.1步长和步长的确定为提高计算精度,时间进程采用蛙跳式方法。例如,求粒子的速度与坐标是按下式求解uβ,n+1/2=uβ,n-1/2+12(Δtn+Δtn-1)duβ,ndt(27)xβ,n+1=xβ,n+uβ,n+1/2Δtn(28)时间步长Δt的选取,由Courant条件确定Δt=ξhC+u|min‚(29)这是对每个粒子计算后取其中最小者,式中ξ是可调因子。2.2粒子界面密度粒子初始配置的任务是确定初始时刻粒子的空间坐标和质量。以后粒子坐标值随时间进程而变化,但质量不变。配置的原则主要有二:一是尽量保持空间分布的均匀性,二是要有一定的空间密度,各粒子之间的间隔一般要求小于或等于一个h值。对于不同构形的模型应采用不同的配置方法。对于具有良好对称性的简单模型,如直角六面体,通常可以直接均匀配置。对于较复杂的构形(如界面有初始扰动),则常需将模型划分成一系列多面体。多面体的中心坐标就是粒子坐标。多面体所含质量为粒子质量。在三维条件下,通常取八点多面体。就一般意义上说,八点多面体应该是十二面体,而且不是唯一的(八点六面体是唯一的,但那是特例)。求网格体积时,需要将它以统一的方式分成六个四面体来计算,以保证与邻网格之间能够密接在一起。见图1。2.3中欧网格划分文献中已说明了搜索近邻粒子的意义和作用。搜索工作借助于欧拉网格,以加快搜索速度。但这类网格的建立仅仅是为了加快搜索速度,它对于SPH方法来说并非是本质性的,所以与常规的拉氏网格法和欧氏网格法都有根本的区别。欧拉网格的划分要覆盖整个流场。为防止粒子流出网格范围,最简单办法是将最外层网格的边缘向外推移至“无穷远”。但这会出现一个网格中集中大量粒子的可能。进一步的处理是,随着流场的变化,随时重新划分网格。但欧拉网格每边长必须保持大于等于2h,以保证能搜索到所有近邻粒子。每个网格周围有26个邻网格,文献已说明只需搜索其中一半即13个网格。对于(i,j,k)网格而言,这13个网格可以这样来取:(i-1:i+1,j-1,k-1:k+1)9个网格,(i-1:i+1,j,k-1)3个网格以及(i-1,j,k)1个网格。2.4虚粒子力学模型的建立至今没有找到关于虚粒子具体使用技术的资料。我们是在选用质量方程为第2节(4′)式的条件下,只在固壁外侧布置了两排虚粒子,而对虚粒子的各力学量则以固壁为对称面作镜面反映式取值。数值试验表明效果良好。3密度和速度剖面内转SPH程序本身是三维的。为了检验其可靠性,曾用它计算过一维高速碰撞问题。这是两块铝板之间的碰撞,碰撞速度为500m/s。与解析解比较,压强、密度和速度剖面图都符合得相当好,但Hugoniot内能偏大。下面就人为粘性、人为热流及能量守恒(包括Hugoniot内能偏大)等问题作些讨论。3.1人为热流误差的影响SPH方法也可以引入人为粘性来处理冲击波间断。人为粘性项用两个可调因子α,β控制其大小。图2是没有人为粘性的结果,压力剖面出现很大振荡,但振荡的幅度没有无限发展。图上虚线是α,β=2.5的波形,作对比用。图3上绘出三组不同α,β值的压力波形。可见当α,β增大时,激波过渡层加宽,且有峰值下降趋势。当α,β减小时,过渡层变窄,但有参量振荡趋势。这些规律是与常规方法一致的。过渡层的一般宽度约占6个粒子,这与差分格式中只占约三个网格相比似乎要宽一倍。但如果考虑到SPH方法中有效作用距离为2h,相当于两个粒子;差分格式中则只有最近邻网格有作用,这样对比,也就相符了。SPH方法引入人为粘性后也存在人为加热误差,在碰撞面上表现尤为突出。因此需要引入人为热流来减少这类误差。人为热流项的作用由两个可调因子g1,g2来调整。图4列出用几组不同g1,g2值计算得出的比内能剖面图。g1,g2=0,碰撞面上内能明显地反常升高。随着g1,g2增大,尖峰不断下降,g1=1,g2=2这一组已基本上将尖峰削平。人为加热误差是因引入人为粘性产生的,所以g1,g2的选取要与α,β的选取相配合。图5显示同一组g1,g2值在不同α〈β值时计算所得的比内能剖面图,可见二者之间存在的依赖关系。3.2激波过渡层仿真两块铝板以500m/s速度相撞,Hugoniot内能SPH计算值与解析解相差达17%。但能量校核表明整个系统的总能量守恒情况良好,相对误差在千分之四以内,见图6。对于此模型,理论上Hugoniot内能应与动能相等。速度值误差很小,表明动能误差很小;总能守恒又好,则内能理应误差也小。进一步检查发现矛盾产生于人为粘性引起的激波过渡层上。图7是粒子能量随粒子分布图,三条实线是计算结果,分别表示总能、动能与内能。虚线是解析解。从图上可见,由于过渡层较宽,动能平台区虽与解析解重合,但系统的总动能计算值偏小许多,所以由增加内能来保证总能守恒良好。图上系统总内能明显偏大,促使平台区内能比解析解高出许多。为进一步说明此问题,做了x方向的粒子数加密4倍后的计算,见图8。由于过渡层变窄,内能误差明显减小了。图8上碰撞面处内能曲线有个突起,说明粒子数加密后,人为热流系数应加大一些,或者人为粘性系数减小一些。若对二者系数作精细调整,予期情况还可能得到改进。4粒子预撞击面本构模型模型是由一个铁质球形弹丸撞击一块铝质靶板。弹丸直径0.635cm,靶板为直角六面体,尺寸是2.5×2.5×0.65cm。取h=0.05cm。靶板中均匀配置粒子,间距为一个h,共布设粒子数为50×50×13=32500,弹丸中近似均匀地配置粒子,径向分7层,共布设粒子339个,粒子间距近似为一个h值。为了程序编制方便起见,总粒子数设置为50×50×14=35000,其中含虚粒子数为35000-32500-339=2161。图9为初始时刻粒子截面图。材料的本构关系简化为用Hugoniot曲线代替。撞击速度是720m/s,撞击方向与靶表面垂线成30°角。计算到15μs,共做了600个时间循环。图10是15μs时的粒子截面图,图11是15μs时靶板撞击面的变形图。程序是在486-DX2/66微机上运行的。运行速率为38min/(10cycle-particles),共化机时约13小时。5u3000eqp的运行速度计算流体力学的三维程序要求计算机有足够的内存和高速运行性能。我们研制SPH程序,开始是在YH-Ⅱ单CPU机上做的,后转移到SUN4工作站,目前在486-DX2/66微机上运行。微机内存是16MB(与YH-Ⅰ的内存相等),只能容纳约35000个粒子。运行速率与本构关系的复杂性、介质被压缩程度等因素有关。微机上运行的最高速度可达35min/(10cycle-particles)。Sun4上运行的最高速率约为75min/(10cycle-particles)。程序向量化后,在YH-Ⅰ上的运行速度为12min/(10cycle-particles)。预计程序并行化后在YH-Ⅱ4cpu机上最高速率将达到1.5min/(10cycle-particles)。用SPH程序模拟三维流体力学界面不稳定性,将需要配置较多的粒子,采用比较符合实际的本构关系,时间循环次数也要多一些。所以实用计算将需要在YH-Ⅱ多cpu机上进行。6存在的问题在实践中,发现还有些问题需要解决,主要是下述两个。6.1sph方法的缺陷激波通过不同介质之间的界面时,发现用SPH方法模拟出来的密度与压强有反常的跳动,见图12和图13。界面位于x=0.62附近。图12上密度反常跳动引起图13上压强的剧烈跳动。初步分析原因有二。一是当界面两侧粒子的质量与密度不同时,计算公式特别是计算密度的公式误差增大,例如公式ρi=mi∑jWij和ρi=∑jmiWij对界面附近粒子计算值会有明显的差异。二是SPH方法中粒子的坐标参量与热力学状态参量都设置在粒子所在空间点上,一旦出现力学量跳动,缺乏