（力学专业论文）基于二阶godunov格式的sph方法的研究.pdf

国防科学技术大学研究生院硕+ 学位论文 摘要 光滑粒子流体动力学方法作为一种无网格l a g r a n g e 粒子方法，特别适合不连 续、大变形和高应变率下的材料断裂等问题的数值模拟。然而传统的光滑粒子方 法存在一些固有缺陷。基于g o d u n o v 格式的光滑粒子法对传统的光滑粒子方法做 了改进，它引入了黎曼解，且计算格式本身具有数值粘性，无须引入人为粘性。 基于g o d u n o v 格式的光滑粒子法具有能提高接触界面计算精度、克服拉伸和压缩 不稳定性现象等优点，它的缺点是自由面计算精度、波阵面分辨率不高。 本文对基于g o d u n o v 格式的光滑粒子法进行了研究，针对它的缺陷提出了一 种新算法。采用二阶g o d u n o v 格式( 粒子内部物理量呈线性分布) 代替一阶g o d u n o v 格式( 粒子内部物理量呈常数分布) 来计算相互作用的粒子之间的黎曼问题初始 间断值，并用t a y l o r 展开对边界进行修正，从而建立了基于二阶g o d u n o v 格式的 s p h 方法，简称为二阶g s p h 方法。本文用新算法对气体的一维不定常流动、多 种冲击加载条件下产生的应力波在固体中的传播进行数值模拟，计算结果表明， 它具有以下优点： 1 ) 能有效提高接触边界的计算精度，克服拉伸和压缩不稳定现象，并提高了 激波的分辨率和计算精度，特别在模拟弹塑性应力波在固体中的传播时改进效果 较明显。 2 ) 在二阶g o d u n o v 格式的s p h 方法中，将t a y l o r 展开思想应用于核估计分 析，提高了该算法处理自由表面的能力，可以在边界上直接加载应力边界条件。 3 ) 由于该方法的计算格式本身具有数值粘性，因此计算过程中无须引入人为 粘性，不仅程序结构简单，并且免去了在其它光滑粒子方法中针对具体问题来选 择人为粘性中的可调参数的步骤，便于程序的调试和推广。 主题词：光滑粒子法；g o d u n o v 格式；t a y l o r 展开；黎曼解；线性函数 第i 页 国防科学技术大学研究生院硕+ 学位论文 a b s t r a c t s m o o t h e d p a r t i c l eh y d r o d y n a m i c s ( s p h ) i s a f u l l yl a g r a n g i a np a r t i c l e m e t h o d t h em e t h o ds e e m si d e a l l ys u i t e d t oh a n d l et h ep r o b l e m st h a tc a n n o tb e e f f e c t i v e l ys o l v e du s i n go t h e rn u m e r i c a lt e c h n i q u e sr e l i e d o nam e s h ，s u c h 觞t h e m o v i n gd i s c o n t i n u i t i e s ，l a r g ed e f o r m a t i o n s ，f r a c t u r e so f m a t e r i a lu n d e rh i g hl o a d i n gr a t e e t c h o w e v e rs t a n d a r ds p hm e t h o dh a ss o m ei n h e r e n ts h o r t c o m i n g s g o d t m o v 。s c h e m e s p hm e t h o dw h i c hu s e sr i e m a n ns o l v e r sh a si n h e r e n tn u m e r i c a lv i s c o s i t y ，t h e r e f o r en o a r t i f i c i a lv i s c o s i t yi sn e e d e da n dh a st h ea d v a n t a g e so fi m p r o v i n gt h es o l u t i o na c c u r a c y a r o u n dt h ec o n t a c ts u r f a c ee f f e c t i v e l ya n dr e d u c i n gt h ei n s t a b i l i t yb o t hi nt e n s i o na n d c o m p r e s s i o n i t sd i s a d v a n t a g ei sl o wa c c u r a c yn e a rb o u n d a r ya n d w a v et i o n t s i no r d e rt oa c h i e v eh i g h e ra c c u r a c yn e a rb o u n d a r ya n ds t e e p e rr e p r e s e n t a t i o no f w a v ef r o n t st h a ng o d u n o v s c h e m es p hm e t h o d ，an e wm e t h o dn a m e ds e c o n d 。o r d e r g o d u n o v s c h e m es p hm e t h o di sp r o p o s e di nt h i sp a p e r p h y s i c a lv a r i a b l e sa r el i n e a r f u n c t i o n si n s i d eap a r t i c l ei nn e wm e t h o dn o tc o n s t a n tf u n c t i o n sa si ng o d u n o v s c h e m e s p hm e t h o d t h e r e f o r et h ei n i t i a lv a l u e so fr i e m a n np r o b l e mb e t w e e np a r t i c l e sa r e c a l c u l a t e di nad i f f e r e n tw a y b o t ho n e d i m e n s i o n a lp r o b l e m so fn o n - s t e a d yf l o wo fg a s a n dt h ee l a s t i c p l a s t i cs t r e s sw a v e si n d u c e db yi m p a c tl o a d i n ga n di m p u l s i v es t r e s s l o a d i n ga r es i m u l a t e dn u m e r i c a l l yb ys e c o n d o r d e rg o d u n o v s c h e m es p h m e t h o d t h e m a i na d v a n t a g e so ft h en e wm e t h o da r es u m m a r i z e da sf o l l o w s ： l1t h er e s u l t ss h o w e dt h a tt h en e wm e t h o dh a dn o to n l ye f f e c t i v e l yi m p r o v e dt h e s o l u t i o na c e u r a c ya r o u n dt h ec o n t a c ts u r f a c eb u ta l s or e d u c e dt h ei n s t a b i l i t yb o t hi n t e n s i o na n dc o m p r e s s i o n a n dw a v ef r o n t sc a l c u l a t e db yn e wm e t h o di sm u c hs t e e p e r t h a nt h o s eb yg o d u n o v s c h e m es p hm e t h o d ，e s p e c i a l l yi ns t r e s sw a v ep r o b l e m s 2 1i n d u c i n gt a y l o rs e r i e se x p a n s i o ni n s e c o n d o r d e ro o d u n o v - s c h e m es p h ， b o u n d a r ya c c u r a c yi se f f e c t i v e l yi m p r o v e d i nt h i sc a s et h es t r e s sb o u n d a r yc o n d i t i o n c o u l db ei m p o s e dd i r e c t l y ，w h i l ei ns t a n d a r ds p h m e t h o dt h eg h o s tp a r t i c l em e t h o dw a s u s u a l l ya d o p t e d t oi m p o s eb o u n d a r yc o n d i t i o n 3 、) t 1 1 es c h e m eo ft h en e wm e t h o dh a si n h e r e n tn u m e r i c a lv i s c o s i t y t h e r e f o r en o a r t i f i c i a lv i s c o s i t yi sn e e d e d ，w h i l ei ns t a n d a r ds p hm e t h o dd i f f e r e n tp a r a m e t e r so f a r t i f i c i a lv i s c o s i t ya r ea d o p t e da c c o r d i n gt od i f f e r e n tp r o b l e m s k e yw o r d s ：s m o o t h e dp a r t i c l eh y d r o d y n a m i c s ，g o d u n o v s c h e m e 。t a y l o r s e r i e s 。r i e m a n ns o l v e r ，l i n e a rf u n c t i o n 一一一一 第i i 页 国防科学技术大学研究生院硕十学位论文 表目录 表3 1 激波管初始参数分布3 0 表3 2 钢板的性质参数3 l 表4 1 s h o c kt u b e 问题的初始参数分布3 6 表4 2 b l a s tw a v e 问题的初始参数分布3 7 表4 3w a l ls h o c k 问题的初始参数分布4 0 表4 4 s j 6 9 r e e n 问题的初始参数分布4 l 表4 5 流体模型材料性质参数4 3 表4 6 弹塑性模型材料性质参数4 5 第1 i i 页 国防科学技术大学研究生院硕十学能论文 图目录 图3 1 l a g r a n g e 网格中的g o d u n o v 间断分解思想。1 5 图3 2 光滑粒子法中粒子相互作用示意图1 9 图3 3r t s 坐标系以及在接触界面a b e 上的应力分布2 l 图3 4 一阶g s p h 方法中迎风格式的思想2 2 图3 5 粒子上物理量的线性分布2 3 图3 6 两种s p h 方法在自由边界的计算精度的比较。3 0 图3 75 l l s 时刻板中应力分布曲线3 1 图3 81 0 9 s 时刻板中压应力的分布曲线3 2 图3 92 2 9 s 时刻板中拉应力的分布曲线3 2 图3 1 0 不同光滑长度，5 9 s 时板中应力波的传播3 3 图3 1 1 不同初始光滑长度下板中应力波的传播3 3 图4 1 0 0 4 1 t s 时刻s h o c kt u b e 问题的参数分布曲线3 7 图4 2 b l a s tw a v e 问题的一阶g s p h 方法的计算结果3 8 图4 3 b l a s tw a v e 问题的二阶g s p h 方法的计算结果3 9 图4 4 加入人为热流后二阶g s p h 方法的计算结果4 0 图4 50 1 ”s 时刻w a l ls h o c k 问题的参数分布曲线4 l 图4 60 0 5 i - t s 时刻s j 5 9 r e e n 问题的参数分布曲线4 2 图4 70 3 1 a s 时刻铝飞片碰撞铝板的参数分布曲线4 4 图4 80 3 i t s 时刻钢飞片碰撞铝靶问题的参数分布曲线4 5 图4 9 粒子数为5 0 0 时一维弹性应力波的传播4 6 图4 1 0 粒子数为1 0 0 0 时一维弹性应力波的传播4 7 图4 1 1 理想弹塑性介质中一维弹塑性应力波的传播4 8 图4 1 2 不同时刻钢板对碰问题的应力分布曲线4 9 图4 1 3o 5 9 s 时刻钢板对碰问题的密度和速度分布曲线5 0 第1 v 页 独创性声明 本人声明所呈交的学位论文是我本人在导师指导下进行的研究工作及取得 的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含 其他人已经发表和撰写过的研究成果，也不包含为获得国防科学技术大学或其它 教育机构的学位或证书而使用过的材料。与我一同工作的同志对本研究所做的任 何贡献均已在论文中作了明确的说明并表示谢意。 学位论文题目：基王三险堡q 垦望垒q y 整盛鲍苎里旦虚洼鲍盟窥 学位论文作者签名：肄监墼址 日期： 加。7 年，月加日 学位论文版权使用授权书 本人完全了解国防科学技术大学有关保留、使用学位论文的规定。本人授权 国防科学技术大学可以保留并向国家有关部门或机构送交论文的复印件和电子 文档，允许论文被查阅和借阅；可以将学位论文的全部或部分内容编入有关数据 库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。 ( 保密学位论文在解密后适用本授权书。) 学位论文作者签名： 作者指导教师签名： ，月加日 ，月跏日 国防科学技术人学研究牛院硕十学位论文 第一章绪论 1 1s p h 方法的背景 随着计算机的飞速发展和广泛应用，以计算机为基础的数值模拟不再仅仅是 理论研究和实验研究的辅助手段，而成为独立于理论和实验的一种科学研究。它 可以用来模拟现实世界的各种过程，部分替代实验或作为实验的补充，检验理论 模型的j 下确性，尤其是还能呈现现实生活中不可重复或无法进行实验的现象，或 模拟耗资巨大的实验工程，从而极大地缩短了研究周期，节省了实验经费，还具 有一定的预测能力。 数值方法在近几十年有了突飞猛进的发展，目前基于网格的数值分析方法可 以分为有限差分法和有限元法两大类。有限差分法的一般求解过程是：首先将求 解域划分为网格，然后在网格的节点上用差分方程近似微分方程。有限元法的基 本思想是：将连续的求解区域离散为数目有限的单元，单元问只在数目有限的指 定点( 节点) 处相互联结，利用每个单元内假设的近似函数来分片地表示全求解域上 的待求的未知函数。将所有单元的这些特性关系按一定的条件集合起来引入边界 条件，构成了一组以节点变量为未知量的代数方程组，从而将一个连续的无限自 由度问题变成离散的有限自由度问题。有限差分法和有限元法已广泛地应用于计 算流体力学和计算固体力学的各个领域，成为数值模拟的主要方法之一。 然而，基于网格的数值方法在很多方面仍存在不足之处，使其在许多问题上 的应用受到了限制。在用拉格朗日法求解会属冲压成形、高速撞击、裂纹动态扩 展、流固耦合、局部化等涉及特大变形的问题时，有限元网格可能会产生严重扭 曲，不仅需要网格重构，而且严重地影响解的精度，导致前处理工作量增大，计 算结果振荡等问题【l 】。尤其对于高维问题，网格重构和细化将大量消耗计算资源， 使大型结构断裂和大变形计算分析等工作难以进行。鉴于网格方法的这些缺陷， 无网格法应运而生。 无网格方法的主要思想是：通过使用一系列任意分布的节点或粒子来求解具 有各种各样边界条件的积分方程或偏微分方程组从而得到精确稳定的数值解，这 些节点或粒子之间不需要网格进行连接。无网格方法基于点的近似，可以彻底或 部分地消除网格，不需要网格的初始划分和重构，保证了计算的精度。 近十几年来国际上许多著名的计算力学学者，对无网格方法进行了大量的研 究工作。典型的无网格法有光滑粒子流体动力学方法( s p h ) 、漫射元法( d e m ) 、无 网格g a l e r k i n 法( e f g ) 、再生核质点法( r k p m ) 、h p 云法( h p c l o u d s ) 、单位分解法 ( p 、局部p e t r o v g a l e r k i n 无网格法( m l p g ) 等掣2 1 。其中光滑粒子流体动力学方 第1 页 国防科学技术大学研究生院硕十学位论文 法( s m o o t h e dp a r t i c l eh y d r o d y n a m i c s ，简称s p h ) 是最早的尤剀格法之一，是由l u c y l 3 】 和m o n a g h a n 4 】等于1 9 7 7 年提出，主要应用于天体物理、流体力学以及两相流等研 究领域。 1 2s p h 方法的发展 光滑粒子流体动力学方法，是一种l a g r a n g e 形式的无网格粒子法，其特点是 在模拟物体大变形时，既可以克服e u l e r 方法难于跟踪物质变形和不能识别材料界 面的缺点，也解决了传统l a g r a n g e 方法在大变形下的网格扭曲问题。基于这些优 点，光滑粒子法现已被广泛地应用于冲击波模拟、流体动力学、水下爆炸模拟、 高速碰撞等材料动态响应的数值模拟等领域。 但是传统的光滑粒子法存在一些固有缺陷，主要表现为边界计算精度低，拉 应力不稳定性，一致性问题，零能模式，和难于处理本质边界条件等问题。国内 外学者对光滑粒子法作了大量的研究工作，主要包括三个方面：一是致力于改进 该方法本身的缺陷，提高精度和计算效率以及为该方法提供略实的理论基础和严 格的数学证明【5 】；二是与网格法相结合和其他数值方法如有限元方法相结合，提高 计算效率；三是不断拓展该方法的应用领域，光滑粒子法已经由最初的天体物理、 流体力学、两相流等相关研究领域向着水下爆炸、材料动态响应、多相流等更广 泛的领域发展。 由于传统的光滑粒子方法不具备数值粘性，所以为了使该方法适合于模拟冲 击波问题，必须引入物理耗散，以消除计算中的数值振荡。1 9 8 3 年，m o n a g h a n l 6 1 等人针对一维激波管问题，提出了适用于光滑粒子法的人为粘性项。1 9 9 7 年 m o n a g h a n t7 j 基于黎曼解思想，通过求解描述粒子间相互运动的特征速度来构造粘 性项。2 0 0 0 年，p a r s h i k o v 8 】等采用黎曼解来描述粒子间的相互作用，提出了粒子 间接触算法，该算法具有数值粘性，在一定程度上改善了数值计算精度和拉引力 不稳定现象。 在光滑粒子方法中，边界上质点处的近似值是区域内部和边界上质点值的线 性组合，因此本质边界条件的处理比较复杂。c a m p b e l l 9 l 重新考虑核估计的积分形 式，通过分部积分中的残余项来引入边界条件，该残余项只在边界的近邻区域内 非零。c h e n 等i lo 】把t a y l o r 展开用于核估计分析上，提出了一种改进的光滑粒子算 法( c o r r e c t i v es m o o t h e dp a r t i c l em e t h o d ，c s p m ) 。 光滑粒子方法中边界计算精度低的另一个主要原因是，该方法要求影响域中 的粒子分布具有完整性和对称性。为了提高边界的计算精度，克服边界粒子的影 响域内粒子分布不对称，l i b e r s k y 1 1 】采用“虚粒子 ( g h o s tp a r t i c l e ) 来描述真实 第2 页 国防科学技术大学研究生院硕十学位论文 粒r f 关于边界对称的镜像，虚粒子具有与真实粒子相i 司大小的密度和速度，但是 速度的方向与真实粒子相反。 s w e g l e 1 2 】采用f o u r i e r 稳定性分析方法，分析了光滑粒子法产生拉应力不稳定 的原因，并指出零能模式由节点积分产生。 b e n z 和w i n g a t e t ”】等将光滑粒子方法扩展应用于模拟脆性材料在受到碰撞后 的断裂。c h e n 1 4 】等采用光滑粒子法进行了复合材料的抗贯穿过程的数值模拟。 l e e 1 5 】等用光滑粒子法研究了陶瓷金属装甲的抗贯穿过程中的子弹消蚀、裂纹传 播、陶瓷锥的形成和背面金属板的破坏等。m e d i m 【m 】等在光滑粒子法中采用并行 算法，开展了高速碰撞和侵彻贯穿方面的数值计算。 为了克服光滑粒子法计算量大、效率低的缺点，a t t a w a y 1 。7 】将光滑粒子法与有 限元方法耦合起来，在大变形、高梯度、移动不连续的地方采用光滑粒子方法， 而其它区域采用有限元方法以提高整体计算效率和方便地处理边界条件。b e n m o u s s a 1 8 】通过引入g o d u n o v 型有限体积方法，使得归一化s p h 方法稳定，并可以 处理边界条件。j o l m s o n 【1 9 】等在初始时刻用有限元来离散整个区域，变形过程中只 有发生扭曲的单元才会转换为s p h 粒子，从而大大提高了计算效率。 2 0 0 6 年，s i g a l o t t i t 2 0 j 等提出了基于自适应核估计算法( a d k e ) 的激波捕捉法。 2 0 0 7 年，o g e r t 2 l j 提出了具有高阶收敛性的归一化的连续性方程和动量方程。2 0 0 9 年，h u 2 2 】等提出了不可压多相流s p h 的常密度法。 国内对光滑粒子法的研究起步较晚，且主要集中在光滑粒子方法的应用上。 张锁春【2 3 】对光滑粒子方法进行了综述，重点介绍了该方法的理论基础，流体动力 学方程组的推导，人为粘性和热传导，自引力、汇和磁场，光滑核的选取，以及 s p h 执行过程等有关问题；贝新源【2 4 】等对于s p h 方法在三维高速碰撞问题中的应 用进行了研究，指出了一些尚需解决的问题；宋顺成【2 5 】等将其与有限元方法结合 计算了侵彻贯穿方面的问题；王裂2 6 1 等对金属表面的微射流进行了模拟；丁桦【2 7 】 等将光滑粒子法进行了改进，用来模拟波在固体中的传播，并对边界条件的表述 和处理进行了探讨；汤文辉、毛益吲2 8 】等在高速碰撞、自由表面流动以及流体r t ( r a y l e i g h t a y l o r ) 不稳定性方面作了有益的尝试；张刚吲3 0 】等提出了一种新的光 滑粒子法的核函数，并将改进的光滑粒子法( c s p m ) 在冲击动力学领域进行了应 用。徐志宏【3 3 】等提出基于一阶g o d u n o v 格式的s p h 改进算法，并采用非迭代近似 黎曼解来描述粒子间的相互作用。韩旭【3 4 】等应用s p h 方法对两相流动问题进行了 研究。 第3 页 国防科学技术大学研究生院硕+ 学位论文 1 3 本文的主要内容 虽然国内外有许多学者在各个方面对光滑粒子方法进行了大量研究，提出了 许多改进的思路，例如：归一化光滑函数算法【3 7 i ( n s f ) ，守恒型光滑化算法，二元 粒子动力学方法【3 7 】，改进的光滑粒子算法【1 0 j ( c s p m ) ，粒子间接触算法等等。这些 改进的光滑粒子法都在一定程度上改善了数值计算精度和拉引力不稳定现象，但 是每种方法仍然存在各自的缺陷，例如当影响域内包含的粒子数目过多时，c s p m 方法存在压缩不稳定性现象；而粒子间接触算法在处理气体的自由边界问题时， 结果很不理想，有明显的反跳现象等等。因此，对s p h 方法有必要从理论上开展 进一步研究。 本文主要针对一阶g o d u n o v 格式的光滑粒子方法( 一阶g s p h 方法) 中自由表 面的计算精度，激波分辨率不高的问题，采用分段线性函数作为粒子内部的基本 插值函数，用来构造黎曼问题的初始间断值，并结合t a y l o r 展开修正边界，建立 了二阶g s p h 方法。本文的具体内容如下： 第一章讨论了基于网格的数值计算方法的特点和局限，简要叙述了无网格方 法的基本思想，特别是光滑粒子流体动力学方法的特性，并叙述了国内外光滑粒 子方法的发展和研究现状。 第二章详细介绍了光滑粒子方法的插值思想及数值计算方法、核函数的性质， 讨论了传统光滑粒子方法的精度；推导了函数导数、矢量的散度和旋度的核估计 形式；运用矢量运算法则推导了连续介质力学守恒方程的多种离散形式，介绍了 可变光滑长度；推导了本构方程离散形式。 第三章首先阐述了g o d u n o v 间断分解思想及黎曼解的两种近似声波近似 和非迭代近似；描述了在光滑粒子法中引入黎曼解的方法，还比较了s p h 方法与 含黎曼解的s p h 方法的动量方程的离散格式的差别。重点阐述了二阶g o d u n o v 格 式的基本思想，将其引入光滑粒子法中，并结合t a y l o r 展开思想在核估计中的应 用，建立了基于t a y l o r 展开的二阶g o d u n o v 格式的s p h 方法，简称为二阶g s p h 方法；推导出了二阶g s p h 方法的基本方程，分析二阶g s p h 方法的数值粘性。 本章最后对二阶g s p h 方法的自由边界计算精度、稳定性进行了验证。 第四章自编程序，用二阶g s p h 方法对气体的一维不定常流动、多种冲击加 载条件下产生的应力波在固体中的传播进行数值模拟。 第五章对全文进行了总结，并提出了需要进一步研究的若干问题。 第4 页 国防科学技术大学研究生院硕十学位论文 第二章光滑粒子流体动力学方法 2 1s p h 方法的概念和基本方程 2 1 1 核函数与核近似 光滑粒子流体动力学方法中，整个流场被离散成大量而有限的粒子。这些粒 子具有质量和体积，并且携带了流场的物理量，它们代表了流场的性质。粒子的 运动由流体动力学基本方程和初始条件及边界条件决定。粒子的运动和携带的物 理量的变化就代表流场的变化。光滑粒子法的核心是插值计算，根据当前时刻粒 子上的物理量，通过核函数插值计算，可得到下一时刻粒子上的物理量。 用贾表示位置矢量，任意函数f ( 2 ) 可利用6 函数表示为： f ( f f ) = l 厂( 元) 万( 孑一j ) 威 ( 2 1 ) r 1 ；一； 其中万 一牙) = ，、1二一二，q 为包含冤的积分域，若在q 内( j ) 是已定义的和连 【u z x 续的，贝j j ( 2 1 ) 式的积分式是精确的或严密的。 核近似也叫核估计，实际上就是用一种性质类比万函数的核函数来进行近似计 算的方法。对于任意的流场函数( 贾) ，若用核函数肜( 元一贾，h ) 取代万函数，则厂( 元) 可以通过核函数近似地表示为( 厂( 元) ) ，即： ( 厂( 元) ) = c ) 厂( 舅。) ( 元一贾，矗) 疵 ( 2 2 ) ( 2 2 ) 式即为函数核估计的表达式。其中核函数形( i i ，h ) 包括两个自变量：位置 矢量贾和光滑长度h 。 为保证插值计算的精度，核函数应具有如下性质： ( 1 ) 归一化条件【形( i - 2 ，办) 威= l ，在解域q 内核函数的积分值为l 。 ( 2 ) 强尖峰函数1 1 婵肜 一i 。，h ) = 万 - 2 ) ，由( 2 1 ) 式有l ，i m f ( 2 ) 。可 知光滑长度h 越小，核估计值越接近真实值。 ( 3 ) 局域性和非负性矿在i = i 处取最大值，且在其影响域内( 即k 一刘t c h ， 通常取誓_ 2 ) 有非负值形( i t h ) 2 0 ，在影响域外为零。 ( 4 ) 递减函数核函数随粒子间距的增加，相互作用应该减弱，所以核函数为距 离防一引的单调减函数。 ( 5 ) 对称性核函数w 为关于i 的偶函数，w ( y - i ，h ) = 形( 1 元一i i ，h ) 。 第5 页 国防科学技术大学研究生阮硕十掌何论文 常用的核函数有： ( 1 ) 样条( b s p l i n e ) 核函数 f l 一1 5 r 2 + o 7 5 r 3 ，o r 1 肚芳p 妁3 4 ，k 甏 其中d 是空间维数，r = k i 7 h ，q 为归一化常数，在一维、二维和三维情形下 分别为2 3 ，1 0 7 r r 和1 刀 。 ( 2 ) g a u s s i a n 核函数 形= ( 去几娟2 其中d 是空间维数，r = k 一7 i h 。 有关其他常用核函数的讨论可参考文献【2 】【3 8 】。 下面讨论核函数估计式( 2 2 ) 式的精度。将函数f ( 2 ) 在贾点进行t a y l o r 展开有： 们h + ( o x 冤吲+ 圭( 警姒冤。吲2 + ( 2 3 ) zd x 代入( 2 2 ) 式，得到： ( 体) ) - 旧+ ( 飘i 橱+ 互1 ( 羚口吲2 + 卜i 威 叫i ) 瞰旅+ ( 套肛司惭赤 q 固 + 互1g 、c 蕊3 z f ：- ) ，肛一牙) 2 肜( 2 - 2 , h ) 赤+ 由于核函数是关于计算点孑的对称函数，同时要求影响域内的粒子i 关于点孑 对称分布，所以在流场内部，( i 一贾) 在核函数影响域内是关于i 的奇函数，因此： 【( i 一牙 一面，办) 威= o ( 2 5 ) 根据核函数的归一化性质和局域性，( 2 4 ) 式变为： ( 饰) ) 叫卅2 以害”4 ) 因此，( 厂( i ) ) 作为函数厂( i ) 的估计值，具有二阶精度。对于边界点，由于( 2 5 ) 式不成立，所以函数f ( 2 ) 的核近似精度是一阶的。 21 2s p h 的基本方程 在s p h 方法中，整个流场由有限个具有质量的粒子表示的。厂( i ) 在粒子f 上 的值为厂( i ) ，记为z ，在粒子上的值为厂( 习) ，记为乃。用粒子的体积巧取 第6 页 国防科字技术人字研冗牛阮坝十字何论文 代( 2 3 ) 式中粒予的尤穷小体兀出7 ，m 粒子的质量m j = 杉乃，其中乃为粒予的 密度( = 1 ，2 ，n ) ，n 为在粒子f 的影响域中的粒子总数。 因此厂( i ) 的核估计表示式可写成离散化的粒子叠加求和形式： ( 厂( 元) ) = 厂( 习沙( 孑一i ，办) 巧 j = i = 芝乃形( 孑一i ，厅) ( 岛巧) =1,-j = 兰孚矿( i i ，乃) j = l 卜j 则粒子f 处函数厂( i ) 的核估计值为： ( ) ：nj m 形( 霉一t ，厅) ( 2 6 ) ，1 = 1 卜l ( 2 6 ) 式说明定义在粒子f 上的任一函数值可通过核函数对其影响域内的所有 粒子相对应的函数值进行加权平均近似。 利用( 2 6 ) 式，可将任意流场参量f ( f f ) 近似地表示出来。例如密度可表示为： p = m ，( 墨一t ，办) ( 2 7 ) ，= i 下面求函数空间导数耵( 贾) 的核估计( 可( 牙) ) 。 记尹= j i ，、一。分别为位置矢量舅、i 的分量，则有嘉= 一嘉，所以： a ( 尹。h )a ( 尹，办) 莎a ( 尹，办) 莎a ( f ，办) 一= 一一= 一= 一一 叙口 莎苏。莎苏h叙台 即： v j ( i i ) = 一v ；( 孟一i ) ( 2 8 ) 利用( 2 2 ) 式、( 2 8 ) 式和g a u s s 公式，( 可( 贾) ) 表示成： ( 耵( i ) ) = v ，t 厂( 孟) 形( 夏一曼，办) 威 = v ； 厂( 曼) 形( i i ，办) 彬一厂( 贾) v ；( i 一贾，办) 废 = 厂( 牙) 形( 冤一是。，办) 元击+ 上厂( 元) v ，形( i i ，而) 威。 其中为q 的边界，亓是面元击的外法向单位矢量。记v = v ，因为核函数在边 界上的值为零，所以上式可进一步化为： ( 夥( 元) ) = 厂( 贾) v 。形( i i ，矗) 威 = v l 厂( 贾) ( i 一贾，乃) 击= v ( 厂( j ) ) ( 2 9 ) 再利用( 2 6 ) 式，可以得到： 第7 页 国防科学技术大学研究生院硕十学位论文 ( 耵( 元) ) 蛳= 知( 暑一t ，办) ( 2 1 0 ) 由( 2 9 ) 和( 2 1 0 ) 式可知，s p h 方法将函数的空间求导转化为对核函数的求导。 由此可见，在s p h 方法中，通过核估计可以将连续函数及其各阶空间导数用离散 点上的函数值的加权平均得到，权因子便是核函数或其相应的导数。 若在( 2 8 ) 式中，取函数厂为常数1 ，则有： o = 兰形( 墨一训 ( 2 1 1 ) i - 、卜j 若核函数满足( 2 1 1 ) 式，则核函数具有一阶一致性。 采用与( 夥( i ) ) 同样的方法，可以将矢量力的散度( v 力) 和旋度( v 彳) 的核估 计的求和形式表示成： ( v 五) ，= 兰等互v 嘶一己瑚 ( 2 1 2 ) j 2 lr j ( v 力) 。= _ 兰等互x v ( 置一i ，办) ( 2 1 3 ) 夥( 墨) 的核估计的其他形式： ( 耵) ，= v 子m ( 乃一，形( 暑一己，办) ( 夥) ，2 吉丢( 乃一，) v 形( 墨一墨，功 c v 厂，= 口兰j = l ，一( 考+ 吾 v c 墨一己，办， ( 耵) ，= 口一l 告+ 告i v ( 墨一己，办) ( v 力) ，和( v j ) 还可表示成： ( v 彳) ，2 善万m jq - 锄v 嘶一办) ( v 彳) ，= 万1 蔷n 一( 互一互) v ( 墨一墨，功 ( v j ) ，2 蔫詈( 刁匈删( ( v 力) ，= 一百1 善n ( 刁一互) v ( i 一弓，厅) ( 2 1 4 ) ( 2 1 5 ) ( 2 1 6 ) ( 2 1 7 ) ( 2 1 8 ) ( 2 1 9 ) ( 2 2 0 ) 第8 页 国防科学技术大学研究生院硕十学位论文 2 2 粒子间的相互作用及其对称性 在s p h 方法中，粒子间的相互作用表现在每个粒子的运动状态由其影响域内 的与其相互作用的近邻粒子的物理量决定，粒子的运动遵循三大守恒定律。 连续介质的质量、动量和能量守恒方程分别为： - - d p ：一叩 ( 2 2 1 ) 一= 一，删 l 二1j f， 口l 丝：一1 盯筇， ( 2 2 2 ) 国p 生：丝u a ，芦 ( 2 2 3 ) d t p 其中口，= x , y ，z ，表示空间坐标轴方向，z ，。是速度矢量历的分量，仃筇是c a u c h y 应力张量方的分量，e 是比内能。 下面对每个方程的离散形式分别进行推导。 2 2 1 质量守恒方程 在s p h 方法中，( 2 2 1 ) 式有两种离散形式。第一种是直接对密度进行离散，即： 上 层= m ，形( 霉一墨，办) ( 2 7 ) j = l ( 2 7 ) 式为密度求和法，它表示粒子密度可以通过其影响域内所有粒子的密度的 加权平均近似求得。 若对( 2 2 1 ) 式右边的速度散度进行离散，直接利用( 2 1 0 ) 式可得至l j ( 2 2 1 ) 式的另 一种离散形式： 鲁一碡考( 硼 ( 2 2 4 ) 通常称之为连续性密度法。它表现了粒子密度变化率和粒子与其影响域内所有粒 子之间的相对速度有关。 密度求和法在整个求解域内积分，严格遵守质量守恒定律，但由于边界点影 响域内粒子分布不对称，边界点计算精度低。对于广义流体问题的模拟，应用经 修正的密度求和法可以得到较好的结果。对于具有强间断问题的模拟( 如爆炸、高 速冲击等1 ，应优先选取连续性密度法【3 9 1 。 2 2 2 动量守恒方程 与质量守恒方程的离散化类似，动量守恒方程的离散也有多种形式。 第9 页 国防科学技术人学研冤生院帧十学何论文 ( 1 ) 直接利用( 2 1 5 ) 式，方程( 2 2 2 ) 右边u 丁以离散为： 丝a c t = 万1 蔷n 鸭( 一塌 ( 2 2 5 ) ( 2 ) 若利用关系式： 一v a ：v ( 丝) + 二v ( p ) pp p 则( 2 2 2 ) 式的离散形式为： 等= 善c 鲁+ 等巩 c 2 2 6 ， 由于( 2 2 6 ) 式具有对称性，能保证系统总动量守恒，而( 2 2 5 ) 式则不能。 通常只要动量方程的离散形式具有对称性就能保证系统动量守恒。动量方程 更一般的离散形式可以写为： 等= 姜c 禹+ 磊瞩 当y = 2 时，上式退化为( 2 2 6 ) 式。 2 2 3 能量守恒方程 能量守恒方程也有多种禺敌形式。 ( 1 ) 直接利用( 2 1 5 ) 式，方程( 2 2 3 ) 右边可以离散为： 生d：譬兰(“歹一甜?嵫tdt p 葺j i ( 2 ) 若利用关系式： - y v j ：v ( 笪) + 彳v ( ) 则( 2 2 3 ) 式的离散形式为： 生d t = 争m 聊i 霉p i ( “；一甜 ( 3 ) ( 2 2 7 ) 式与( 2 2 8 ) 式相加得到更普遍的离散形式： 鱼a t = 争哆睁种m 峨 离散形式( 2 2 9 1 式具有对称性，因此能够确保系统的总能量守恒。 2 2 4 可变光滑长度 ( 2 2 7 ) ( 2 2 8 ) ( 2 2 9 ) 在光滑粒子法中，光滑长度h 非常重要，直接影响到求解的精度与效率。光滑 长度的选取应满足充分覆盖求解域的要求。简单地说，光滑长度选取同粒子局部 第1 0 页 国防科学技术大学研究生院硕十学何论文 密度有关。通过改变光滑长度的人小，以保证影响域内的粒子数量相对不变。对 于均匀变化的变量场，选用稍大的光滑长度也能满足对变量场的近似精度，但对 于变化剧烈的变量场，选取较大的光滑长度将造成局部失真。若选取较小的光滑 长度近似程度明显提高。另一个需要考虑的因素是计算效率，光滑长度较大，影 响域内粒子数多，计算量就大，计算效率就低。因此选取合适的光滑长度十分重 要。 当流体变形很小时，h 可近似地取为常数。对于局部膨胀或局部压缩等大变形 问题，为保持全局精度一致性，必须对光滑长度进行修正，使每个粒子影响域内 的粒子数相对保持不变。 最简单的方法为通过影响域内粒子质量保持不变，对光滑长度进行动态变换： i 时 ( 2 3 0 ) 、岛分别为初始光滑长度和密度，d 为空间维数。 b e n z ( 1 9 8 9 ) 提出另一种光滑长度的动态变换法，在每个时间步，求光滑长度随 时间的导数，用于计算下个时问步的光滑长度： 塑：一! 鱼塑( 2 3 1 ) 一= 一一 i 二jl ， d tdpd t 利用连续性方程，可以将上式写成离散形式： 一d h , ：一t 7 , 厶u , g 了m j ( 甜歹一u ；) e j ，(232)dt d鼍pj 、 诤 由( 2 3 0 ) 式和( 2 3 1 ) 式可知，光滑长度是时间和空间的函数，所以每个粒子都具 有独立的光滑长度。假设粒子f 位于粒子的影响域内，若瑰h j ，则粒子不一 定位于粒子f 的影响域内，即粒子f 对粒子，施加作用力，而不会受到粒子的反作 用力，这与牛顿第三定律矛盾。此外，在验证总动能和总能量守恒时，用到了 = ，也即缈( 1 i 一霉l ，曩) = 形( e i i ，一) 。因此必须保证光滑长度的对称性。 实现粒子间相互作用的对称性，通常有两种方式： ( 1 ) 采用对称形式的光滑长度来代替曩和办，其中有多种计算方法，例 如：( h i + h i ) j 2 ，m