（机械设计及理论专业论文）一种改进的mlsph方法理论研究及其应用.pdf

硕f 学位论文 摘要 光滑粒子流体动力学( s p h ) 方法是一种典型的无网格拉格朗同方法， 该方法只需要节点信息，不需要背景网格，可以避免基于网格的拉格朗 日方法中的网格畸变，适宜处理高速冲击、大变形等问题。但该方法引 入到固体中时，遇到两个重要的问题：一是边界处理问题；二是拉伸不 稳定性问题，即在拉伸状态下会产生非物理的粒子聚集现象。为了改善 拉伸不稳定性，将移动最小二乘法引入到s p h 方法中，进而演化出移动 最小二乘粒子流体动力学( m l s p h ) 方法。 本文对m l s p h 方法进行了详细研究，主要工作如下： 首先，基于m l s p h 方法的基本理论，自主开发了一维、二维m l s p h 程序，并通过分片试验验证了程序的有效性；讨论分析了m l s p h 方法的 稳定性，并结合s w e g l e 扰动数值分析得出，m l s p h 方法比s p h 方法具 有更好稳定性。同时，还讨论了ml sp h 方法中一些问题，如人工粘性、 光滑长度和相邻粒子搜索法等。 其次，分析了m l s p h 方法存在的问题，并有针对性地提出了改进措 施；通过引入边界向量，给出了三种不同形式的离散方程组，并讨论了 它们的守恒性；提出了隔步搜索法以提高m l s p h 方法的计算效率，并通 过数值分析验证了该方法的可行性以及所开发的基于隔离搜索m l s p h 程序的有效性。 最后，利用m l s p h 方法对激波管问题进行研究，得到了较好的结果。 并利用改进的m l s p h 方法对一些由拉伸起主导作用的问题进行应用研 究，并与s p h 方法结果对比得出，m l s p h 方法可以改善s p h 方法的拉 伸不稳定性，而且采用隔步搜索法后能有效地提高原m l s p h 方法的计算 效率。 关键词：光滑粒子流体动力学( s p h ) ；移动最小二乘粒子流体动力学 ( m l s p h ) ；隔步搜索法；拉伸不稳定性；数值模拟 a b s t r a c t 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 ) m e t h o d i sap u r em e s h l e s s l a g r a n g i a nt e c h n i q u ew h i c hi sa p p e a l i n g a san e wn u m e r i c a im e t h o du s e dt o s o l v et h ep r o b l e m s w i t h h i g h v e l o c i t y a n dl a r g e d e f o r m a t i o n t h e a d v a n t a g e so fs p ha r e ：( 1 ) o n l yn o d a ld a t a a r en e c e s s a r y ；( 2 ) b a c k - g r o u n d g r i dd o e sn o tb en e e d e da n ds ot h es e v e r ep r o b l e m sa l w a y s a s s o c l a t e dw l t n m e s hd i s t o r t i o n i nt h el a g r a n g i a ng r i d b a s e d m e t h o d sc a nb e a v o i d e d h o w e v e r t h el o wc o m p u t a t i o n a la c c u r a c y a tt h eb o u n d a r ya n dt h et e n s i l e i n s t a b i l i t ya r ee x i t e d ，w h e n s p hm e t h o d i su s e dt os o l v e t h es o l i d m e c h a n i c s ，p r o b l e m s n o n p h y s i c a lp h e n o m e n a o ft h ep a r t i c l e sa g g r e g a t l o n a p p e a r si nt e n s i l e s t a t e i no r d e rt oi m p r o v et h e t e n s i l ei n s t a b i l i t y ， t h e m o v i n g l e a s t s q u a r e si n t e r p o l a n t m e t h o di sa p p l i e dt ot h es p hm e t h o d , w h i c he v 0 1 v e sm o v i n gl e a s ts q u a r e sp a r t i c l e h y d r o d y n a m i c s ( m l s p h ) m e t h o d t h i sp a p e ri sf o c u s e do nt h es t u d yo fm l s p h m e t h o d s ，t h em a l nw o r k s i nt h ep a p e ra r e ： f i r s t l y ，1dm l s p hp r o g r a ma n d2 dm l s p hp r o g r a m a r ed e v e l o p e do n t h eb a s i co ft h et h e o r yo ft h e m l s p hm e t h o d t h ee f f e c t i v i t y o ft h e p r o g r a m si s v e r i f i e db yp a t c ht e s t t h e n t h es t a b i l i t yo ft h em e t h o d 1 s i n v e s t i g a t e d a n dn u m e r i c a ls i m u l a t i o no fs w e g l ed i s t u r b i n g t e s t1 sa p p l l e d ， w h i c hs h o w s t h a tm l s p hm e t h o d i sm o r es t a b l e t h a ns p hm e t h o d m e a n w h i l e ，s o m en u m e r i c a la s p e c t s s u c ha sa r t i f i c i a lv i s c o s i t y ，v a r i a b l e s m o o t h i n gl e n g t h ， n e a r e s tn e i g h b o rp a r t i c l e s e a r c h i n ga l g o r l t h m s a r e d i s c u s s e d s e c o n d l y ，t h es h o r t c o m i n go fm l s p h m e t h o da r ei n v e s t i g a t e di nd e t a i l , a n ds o m ee f f e c t i v em e a s u r e sa r ep r o p o s e d t h r e ed if f e r e n tf o r m so fd i s c r e t e e q u a t i o n sa r eg i v e nt h r o u g ht h ei n t r o d u c t i o no f t h eb o u n d a r yf l u x e s ，a n dt h e c o n s e r v a t i o no ft h e ma r ed i s c u s s e d t h el e a p i n gs e a r c hm e t h o d i sp r e s e n t e d t oe n h a n c et h ee f f i c i e n c yo fm l s p hm e t h o d ，t h ee f f e c t i v i t yo f t h i sm e t h o d a n dt h ep r o g r a ma r ev e r i f i e db ys o i b _ en u m e r i c a lt e s t s f i n a l l v ，t h e n n u m e r i c a le x a m p l e w i t hs h o c kt u b ep r o b l e m i s s a t i s f a c t o r yb y m l s p hm e t h o d s o m ep r o b l e m s w i t h s t r e t c h i n g a r e 硕l j 学位论文 i n v e s t i g a t e db ym o d i f i e dm l s p hm e t h o d ，w h i c hg e tm u c hb e t t e rr e s u l t c o m p a r e dw i t hs p hm e t h o d ，t h es i m u l a t e dr e s u l t s s h o wt h a tm l s p h m e t h o dc a ni m p r o v et h et e n s i o ni n s t a b i l i t yo fs p hm e t h o de f f e c t i v e l y ，a n d t h ee f f i c i e n c yo fm l s p hm e t h o di se n h a n c e db yl e a p i n gs e a r c hm e t h o d k e y w o r d s ：s m o o t h e dp a r t i c l e h y d r o d y n a m i c s ( s p h ) ；m o v i n g l e a s t s q u a r e sp a r t i c l eh y d r o d y n a m i c s ( m l s p h ) ；l e a p i n gs e a r c h m e t h o d ；t e n s i l ei n s t a b i l i t y ；n u m e r i c a ls i m u l a t i o n i v 硕i j 学位论义 第1 章绪论 1 1 引言 近几十年来，有限元方法。1 已成为解决科学和工程问题的主要数值 手段。目前，人们已经成功地开发了大量的有限元商业软件，并在工程 分析和计算中得到了广泛的应用。 有限元法求解问题的基本思想是：先将一个表示结构或连续体的求 解区域离散为若干个子域( 单元) ，并通过它们边界上的结点相互联结成 为组合体；利用每个单元内所假设的近似函数来分片地表示整个求解区 域上待求的未知场变量，而每个单元内的近似函数由未知场函数及其导 数在单元各个结点的数值和与其对应的插值函数来表达；由于在联结相 邻单元的结点上，场函数应具有相同的数值，因而将它们作为数值求解 的基本未知量，这使得求解原来待求场函数的无穷多自由度问题转换为 求解场函数结点值的有限自由度问题；然后通过和原问题数学模型( 基本 方程、边界条件) 等效的变分原理或加权残值法，建立求解基本未知量( 场 函数的结点值) 的代数方程组或常微分方程组；最后通过一定的算法求解 此方程组，得到问题的解答。 虽然有限元法通用、灵活，并广泛应用于工程分析和计算中，但是 随着计算对象复杂程度的增加和应用工作的深入，也逐渐显现了其本身 难以克服的一些不足，如所求解函数的导数精度低、高阶不收敛、自锁 问题等，而且在有限元法中，单元和网格既是分析解决问题的载体，同 时也是对其应用的限制。 一般地，有限元法有两种类型：l a g r a n g e 型和e u l e r 型。l a g r a n g e 型有 限元法，在整个计算过程中网格是附着或者固定在物质上的，网格会随 着物质的运动而运动。这使得在求解金属冲压成形、高速冲击、裂纹动 态扩展、流固耦合、局部化等涉及特大变形的问题时，网格可能会产生 严重扭曲，不仅需要网格重构，而且会严重地影响解的精度；当网格变 形太大的时候，计算精度和求解都会受到很大的影响；另外，由于时间 步长是由最小单元尺寸控制，若网格太小，则会影响时间历程的效率， 甚至会导致计算失败。对裂纹的动态扩展问题，由于裂纹的扩展方向不 能事先确定，因而在计算过程中需要不断地重划分网格以模拟裂纹的动 态扩展过程，这大幅度地增加了计算工作量，降低了计算效率；有限元 近似基于网格，因此必然难于处理与原始网格线不一致的不连续性和大 一种改进的m i 。s p h 办i ：理论研究及j e 心用 变形；单元网格划分等前处理数据准备工作量大，尤其是对复杂三维问 题，其工作量往往比有限元分析本身还大。 e u l e r 型有限元法，整个计算过程中网格是固定在模拟对象所处的空 间上，模拟对象在固定网格单元上运动，虽然可以很好地处理物体大变 形问题，但是很难分析物体上固定点的场变量的时间历程；对于材料或 者介质几何图形不规则或复杂的问题，很难处理，往往需要将不规则几 何形状问题转化为规则的计算域；在初始布置网格时，需要足够大以至 于能覆盖到物质流动的整个区域，使得有时为了提高计算效率而使用较 为粗糙的网格，因而降低了离散化区域的分辨率和计算精度o 鉴于有限元在这些方面存在缺陷，近几年来国际上许多著名的计算 力学学者都对无网格方法进行了大量的研究工作口。83 。无网格方法采用基 于点的近似，可以彻底或部分地消除网格，不仅可以保证计算的精度， 而且可以大幅度减小网格生成工作量。无网格方法的主要思想是：采用 一系列任意分布的节点( 或粒子) 离散所求解问题区域，通过一种与权函 数( 或核函数) 有关的近似，使某个域上的节点可以影响研究对象上任何 一点的力学特性，这些节点或粒子之间不需要网格进行连接。无网格方 法在涉及网格畸变、网格移动和不确定边界等问题中显示出明显的优势， 被认为是一种很有发展前途的数值分析方法1 。 1 2 无网格方法的发展与研究现状 无网格法的研究可以追溯到2 0 世纪7 0 年代初对非规则网格有限差分 法的研究h 5 1 ，但由于当时有限元法的成功，这类方法没有得到重视。19 7 7 年l u c y 6 1 ，g i n g o l d 和m o n a h a n t 7 】分别提出了光滑粒子流体动力学法 ( 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 ) ，并成功应用于天体物理方面的 研究中。不过，到2 0 世纪8 0 年代末，无网格法都没有突破性地进展。直 到19 9 4 年b e l y t s c h k o 9 】在修正了d e m ( d i f f u s e re l e m e n t sm e t h o d s ) 8 】的基 础上提出了e f g ( e l e m e n t sf r e eg a l e r k i n ) 法后，无网格法才得到迅速的发 展。 19 81 年，l a n c a s t e r 等【10 1 为了从散乱数据点中拟合曲线和曲面，提出 了移动最小二乘法( m o v i n gl e a s ts q u a r e s ，ml s ) 。m l s 是无网格方法中构 造近似函数的一种重要方法，对无网格方法的发展起了积极有效的推动 作用。 l9 9 2 年，n a y r o l e s 等i s 】人将m l s 方法引入到g a l e r k i n 方法中，提出了 离散元法( d i f f u s e e l e m e n t sm e t h o d ，d e m ) 。l 9 9 4 年，b e l y t s c h k o 等对 d e m 进行了有效的改进，即在计算形函数导数时保留了被n a y r o l e s 等忽略 形! i j 学位论文 掉的项，并利用l a g r a n g e 乘子法对边界条件处理，采用高阶高斯积分完成 区域积分等方法，提出了无单元g a l e r k i n 法( e l e m e n t s f r e eg a l e r k i n m e t h o d s ，e f g m ) ，并成功地解决一些有限元法不能很好解决的问题【l 1 ，l 2 1 。 l i u 等【1 3 】将e f g 和边界元法相耦合，用于固体的应力分析：b e l y t s c h k o 和 h e g e m 等i 1 4 ，15 1 将e f g 法和有限元法耦合，以发挥各自的优势；张雄等【16 1 将e f g 法的思想应用于节理岩体的分析中；为了避免使用背景网格， b e i s s e l 等【l7 1 提出了节点积分方案，但计算稳定性较差；周维垣等【18 】对e f g 法进行了详细地介绍，并应用于裂纹扩展分析中。这些研究结果表 明，e f g m 法精度和收敛速度都高于有限元法，且没有体积锁死现象，但 e f g 法计算量大，并且需要借助背景网格进行数值积分。 19 9 5 年，o d e n 和他的学生d u a r t e 等【”】利用m l s 法建立单位分解函数， 并构造权函数和试函数，再通过g a l e r k i n 法建立离散模型，进而提出了h p 云团无网格法( h pc l o u d sm e t h o d s ，h p c m ) ，并对此进行了严格的数学证 明。随后，o d e n 等 2 0 】又将有限元形函数作为单位分解函数，提出了基于 云团法的新型h p 有限元( n e wc l o u d s b a s e dh pf e m l 。l i s z k a 等【2l 】采用配 点格式，避免了g a l e r k i n 格式中用于数值积分计算的背景网格，提出h p 无网格云团法( h pm e s h l e s sc l o u d sm e t h o d ，h p m c m ) ，它是一种纯无网 格方法。 19 9 6 年，o n a t e 和z i e n k i e w i c z 等 2 2 - 2 4 】合作，利用基于高斯权函数带权 正交m l s 来构造近似函数，并采用配点型加权残值格式把支配方程离散 成非积分的形式，再结合广义有限差分法，提出了有限点法( f i n i t ep o i n t m e t h o d ，f p m ) ，并成功地应用于计算流体动力学问题。有限点法( f p m ) 由于采用配点格式对求解域进行离散，不需要背景网格和高斯积分，计 算效率较高，主要应用于流体动力学领域。宋康祖等【25 】将其应用于弹塑 性分析中。 19 9 8 年，a t l u r i 和z h u 等 2 6 , 2 7 】基于m l s 近似函数和局部 p e t r o v g a l e r k i n 法离散，提出了无网格局部p e t r o v g a l e r k i n 法，并用于求解 带调和算子的l a p l a c e 方程和p o i s s o n 方程、稳态不可压流的n a v i e r s t o k e s 方程、应变梯度材料等问题，并对m l p g 方法进行了误差分析。2 0 0 2 年， a t l u r i 币t ls h e ns p 2 8 , 2 9 1 通过采用不同的近似函数和加权函数构造了六种不 同的无网格局部p e t r o v g a l e r k i n 方法，并对其计算效率和精度进行了深入 研究，进一步发展了m l p g 法的理论和应用范围。l i u 等【3 0 】将ml p g 和有 限元及边界元相耦合，充分发挥它们各自的优势；张见明等【31 ，3 2 1 将用于 杂交边界元的修正变分原理与移动最小二乘近似结合，提出了杂交边界 点法，其输入数据只是求解域边界上的离散点。 一种改进的m l s p h 方法胖沦研究肢j cj 世用 2 0 01 年，陆明万、张雄等 3 3 - 35 】在f p m 法格式的基础上引入辅助点和 采用加权残值法、最小二乘法进行求解，提出了最小二乘配点无网格法 ( l e a s t s q u a r e sc o l l o c a t i o nm e s h l e s sm e t h o d ，l s c m ) $ 1 1 加权最小二乘配点 无网格法( w e i g h t e dl e a s ts q u a r e sc o l l o c a t i o nm e s h l e s sm e t h o d ，w l s c m ) ， 同时还在子域插值和配点法的基础上提出分阶拟合直接配点无网格法。 这三种配点型无网格法可以很好地解决f p m 的不足。 2 0 0 6 年，王新宇 3 6 1 和马上【3 7 】等分别对物质点法( m a t e r i a lp o i n t m e t h o d ，m p m ) 进行了研究，并应用于冲击领域。m p m 法是在质点网格法 ( p i c ) 基础上发展而来的，它利用了欧拉法和拉格朗日法两者的优点，仅 简单地划分背景网格，物质点信息通过形函数映射到背景网格节点，然 后根据相应的控制方程和本构关系，计算物质点在冲击载荷下的应力和 应变，通过物质点来跟踪材料体的变形和破损，但在整个计算过程中背 景网格始终固定不变，避免了重新划分网格。 1 3 光滑粒子法的发展与研究现状 目前无网格方法已有3o 余种，并且还将有更多的新的无网格方法出 现，其中s p h 方法是最早的无网格方法之一，其发展已有三十多年。s p h 方法在求解过程中不需要任何网格，其通过核函数近似进行插值计算， 将流体力学方程或固体方程组转化为s p h 方程组，适宜于处理流体、侵彻 和爆炸等各种大变形问题。 2 0 世纪9 0 年代，m o n a g h a n 38 ，3 9 】在对s p h 方法深入研究后，提出了适 用于s p h 方法的人工粘性项，可以很好地模拟流场中的强激波间断现象； 随后又提出了适用于s p h 方法的人工热量项，可以减小由于引入人工粘性 项后所产生的某些误差。j o h n s o n l 4 0 】在计算中用六种不同形式的人工粘 性，研究光滑粒子法中人工粘性的效果，结果发现，在某些情况下人工 粘性对计算结果有十分明显的影响。人工粘性系数通常是常数，但m o r r i s 等【4 1 】提出了一种系数随时间变化的人工粘性，在波阵面引入人工粘性， 而衰减项使得波后人工粘性迅速减小。王肖钧等 4 2 】采用通量修正输运法 ( f c t ) 代替人工粘性法，在处理激波振荡方面取得了更好的效果。 19 9 0 年，l i b e r s k y 等【4 3 】在s p h 法中引入了偏应力项，并考虑了简单的 弹塑性本构关系，从而率先将材料强度效应引入到s p h 方法，成功地丌展 了高速碰撞数值模拟计算。j o h n s o n 等人 4 4 , 45 】提出了归一化光滑函数算 法，提高了s p h 法的精度，使其能够通过分片试验，可以正确模拟常应 变状态；并将s ph 方法和有限元方法耦合起来，发挥它们各自的优势，即 在变形大的地方采用s p h 法，开展侵彻方面的数值计算，得到了一些有意 - 4 一 硕l j 学位论文 义的结果。贝新源”6 】开发了三维s p h 程序，应用于高速斜侵彻碰撞问题； 韩旭、杨刚等 4 7 , 4 8 】利用s p h 方法模拟两相流动和近水面爆炸等问题，取 得了较好的结果。 然而，s p h 方法在处理高速碰撞冲击等固体领域内问题时存在严重的 缺陷，这主要有两个方面：一是边界处理问题，例如边界施加困难而导 致自由面的计算精度难以保证；二是拉伸不稳定问题。l9 9 5 年，s w e g l e t 4 9 】 运用v o nn e u m a n n 稳定性分析方法对传统s p h 方法的稳定性进行了比较 细致的研究，指出这种不稳定性仅与核函数的二阶导数与应力乘积的符 号有关。d y k a 5 0 ，5 1 】等人认为s p h 方法的拉伸不稳定性与有限元法所遇到 的相类似，因而采用与有限元处理该问题相类似的方法，在s p h 的粒子间 额外的插入积分点，即应力点法( s t r e s sp o i n t ) ，解决了标准s p h 的不稳定 性问题。随后，c h e n 等【5 2 】把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 e m e t h o d ，c s p m ) ，有效地解决了计算精度、相容性、拉应力不稳定等一系 列实际问题。随后，b o n e t 等 5 3 , 5 4 采用核修正和积分修正得出了完备的改 进的光光滑粒子法( c s p h ) ，改善了计算精度和拉伸不稳定现象。m o n a g h a n 和g r a y 等【5 5 】则认为，如果固体处于拉伸状态时，内部原子之间就会产生 相互的吸引力以抵抗拉伸，于是提出了人工应力法( a r t i f i c i a ls t r e s s ) 来解 决s p h 方法的拉伸不稳定。s w e g l e 和h i c k s 等人【56 】采用守恒光滑法 ( c o n s e r v a t i v es m o o t h i n g ) 来控制s p h 拉伸不稳定性。l i u 等人【57 】提出再生 核粒子法( r e p r o d u c i n gk e r n e lp a r t i c l em e t h o d s ) 以改善s p h 法的拉伸不稳定 性。19 9 9 年，d i l t s 58 ，5 9 】利用移动最小二乘法( m l s ) 进行场函数插值，进 而提出移动最小二乘粒子流体动力学法( m l s p h ) 来改善s p h 方法的拉伸 不稳定性等。 1 4 本文的选题依据 随着科技的发展，数值计算越来越被人们重视，要求也愈来愈高。 尽管有限元方法非常的成熟，并广泛应用于各个领域，然而在模拟爆炸、 冲击、大变形、大梯度、移动材料界面和自由表面等问题时存在一定的 困难和局限。 无网格方法的提出，使得网格可以彻底或部分地消除，其采用基于点 的近似，只需节点信息，不需将节点连成单元。光滑质点流体动力学是 最早出现的无网格方法，其突出特点是计算空间导数时不需要使用任何 网格，而是通过一个称为“核函数”的积分进行“核函数估值”近似，将动力 学基本方程转换成离散形式的近似方程，从而避免了高维l a g r a n g e 网格方 种改进的m l s p h 方法理论研究成j cj ! 用 法中的网格缠绕和扭曲等令人头痛的问题。近三十年来，国内外的研究 表明：作为一种真正的无网格方法，s p h 方法在解决结构高度大变形和界 面不连续等问题时具有独特的优势，尤其适合于求解高速碰撞等大变形问 题。 然而，在s p h 方法求解固体问题时，它存在了两个重要的缺陷：一是 边界处理困难；二是拉伸不稳定性。s p h 方法求解固体问题时，在拉伸状 态会下产生非物理的粒子聚集现象，即所谓的拉伸不稳定性；尤其是在 求解由瞬态拉伸应力区产生的损伤、断裂、破片等问题时，拉伸不稳定 性问题显得尤为突出。这两个缺陷严重地阻碍s p h 方法应用于固体领域， 因而是必须要克服的。为了解决拉伸不稳定性问题，d i l t s 5 8 , 5 9 1 将移动最 小二乘方法应用到粒子流体动力学法中，进而提出了移动最小二乘粒子 流体动力学法。该方法采用移动最小二乘法进行插值计算，不仅可以提 高了计算精度，而且可以很好地改善s p h 方法的拉伸不稳定性。但该方法 由于采用m l s 插值计算，则在计算过程中需要求解逆矩阵，大大增加了 计算量，降低了计算效率；而且关于该方法在固体问题中的应用研究较 少。由此可见，改善s p h 方法的拉伸不稳定性，提高m l s p h 方法的计算 效率以及推广m l s p h 方法在固体问题中的应用仍需要深入研究。 总结以上论述可以看出，s p h 方法在解决大变形问题时所具有独特 性的优势和m l s p h 方法所具有的优点，使得m l s p h 方法可以更广泛地应 用于冲击等固体领域，这也是选本题的意义所在。 1 5 本文的主要研究内容 基于上述，本文对m l s p h 无网格方法及其应用进行了详细研究，主 要内容有： 在第1 章中，调研了国内外大量关于无网格方法的相关文献，对该方 法的研究进展和发展趋势有了深入的了解；回顾了s p h 方法的发展历史及 其研究进展；阐述了本文的选题依据及主要内容。 在第2 章中，阐述了滑粒子流体动力学方法和移动最小二乘粒子流体 动力学方法的基本理论。结合m l s 法，详细推导了ml s p h 方法的离散表 达式。讨论分析了m l s p h 方法的稳定性、人工粘度、光滑长度的更新及 相邻粒子对的搜索方法等相关问题。利用f o r t r a n 语言开发了一维、二维 ml s p h 程序，并进行了分片试验、s w e g l e 数值稳定性试验。 在第3 章中，分析了m l s ph 方法存在的问题，并给出了相应的改进措 施；针对计算效率较低问题，提出了隔步搜索方法以提高效率。阐述了 隔步搜索法的基本思想及其具体实施过程。丌发了基于隔步搜索的 硕i j 学位论文 m l sp h 程序，并进行了验证。 在第4 章中，利用改进的m l s p h 方法对激波管问题和一些由拉伸起主 导作用问题进行了应用研究，并与s p h 方法的计算结果进行对比。 最后总结了全文的工作及展望了今后的研究工作。 一种改进的m l s p h 方法理论研究及j c 膨用 2 1 引言 第2 章m l sp h 方法及其程序实现 光滑粒子流体动力学( s p h ) 方法是l u c y t6 】和g i n g o l d t7 】等在l9 7 7 年分 别提出的一种纯l a g r a n g e 无网格方法。该方法适宜于处理高速冲击、碰 撞、爆炸以及各种大变形问题。然而s p h 方法在计算固体方面时，遇到 了两个重要的问题：一是边界处理问题；二是拉伸不稳定性问题【4 引；尤其 是在处理由瞬态拉伸应力区产生的损伤、断裂、破片等问题，拉伸不稳 定性问题显得尤为突出。为了解决该问题，d i l t s 5 8 , 5 9 】将移动最小二乘方 法f m l s ) 应用到s p h 方法中，进而提出移动最小二乘粒子流体动力学法 ( m l s p h ) 。该方法采用移动最小二乘法进行插值计算，不仅可以提高s p h 方法的精度，而且很好地改善s p h 方法的拉伸不稳定性。 2 2s p h 方法的基本理论 s p h 方法的主要思想是将连续介质离散成任意分布的质点，所有物 理量都集中于质点上，未知函数的数值通过核积分的形式来表示。其中， 核函数的物理意义为粒子物理量分布的概率密度，而数学意义上为未知 函数的光滑函数。 2 2 1s p h 方法的基本方程 s p h 方法有两个关键步骤：第一步积分表示法，即场函数核近似法； 第二步粒子近似法。 s p h 方法的核心就是插值理论。这种方法在计算空间导数时不需要 任何网格，而是通过核函数的积分核进行核函数估值近似的，将连续性 介质的偏微分方程转化为积分形式。整个物体被离散成一系列粒子，所 有物理量( 密度、压力、速度、体积、能量等) 由这些粒子携带。由于计算 时只涉及各个离散粒子的信息，所以积分必须要通过对邻近粒子的求和 得到。 在数学上，对于任意连续函数g ( 工) ，它可以通过狄拉克万函数表示为 g ( x ) = i g ( x ) 6 ( x - x ) d x ( 2 1 ) 硕i j 学位论文 式中，狄拉克万函数为万。一砷= 0 1z x i x x ，l 在s p h 方法中，任意函数( 工) 在空间某一点x 处的核估计值都可以通 过函数厂( x ) 在域q 中的积分获得 ( 工) i 厂( 工7 ) 形( 工一z ， ) 出7 ( 2 2 ) 矗 式中，w ( x 一工，j i ) 为核函数，h 为定义核函数影响域的光滑长度，出为小体 元。该式对矢量和张量均适用。 空间导数v ( 石) 的近似式可利用( 2 2 ) 式，用v 厂( j ) 代替厂( 功得到 v 厂( 功i ( v f ( x ) ) 矽( x x ，j 1 1 ) 出 ( 2 3 ) 矗 式中，v 为哈密顿算子。 根据散度定律：v ( 州= u v v + v v u ，式( 2 3 ) 变为， v ( 功v ( 厂( x ) 形。一x ，1 1 1 ) ) 出一p ( x ) v w ( x x ，| 1 1 ) 出 ( 2 4 ) 再根据高斯定理： v 力f = 扣西= 扣再嬲，将( 2 4 ) 式的等号右边第一 项在体积域q 上的积分转化为在面积域s 上的积分，即 ，v ( f ( x ) w ( x - x ，j i i ) ) 出= f ( x ) w ( x - x ，1 1 ) 亓据 ( 2 5 ) 如果积分域位于计算域内，且与材料边界不相交，则( 2 5 ) 式的面积 积分项为零。如果积分域与材料边界相交，则需要另对边界条件进行处 理。 于是，函数的最终微分形式为： v 厂( 力= 一f f ( x ) v w ( x x ，j l i ) 出( 2 6 ) 矗 因此，函数导数的核估计值可以通过函数的值和核函数的导数来确 定，而不是由函数本身的导数来确定。 在s p h 方法中，整个系统是由具有独立的质量、占有独立空间的有 限个粒子表示的。若使用粒子的体积矿，取代积分中粒子处的无穷小体 元出，则粒子的质量m ，可表示为： m = p j ( 2 7 ) 式中，p ，为粒子，的密度( 歹= 1 , 2 ，忉，其中n 为粒子的支持域内的粒子总 量。 于是，( 2 2 ) 和( 2 6 ) 式的连续方程可以分别写为粒子的求和形式： 种改进的m l s p t t 方法理论研究及j cj 避用 ! 一! ! - ! ! ! ! ! ! ! ! ! ! ! = = = ! = = = = = = = = = = = = = = = = = = = = = = ! ! ! = = = ! ! ! = = 竺! = = = ! = = = ! = = ! = = = = m ) = 姜考似 ( 2 8 ) 口m ) = 姜考似v ( 2 9 ) 式中w u = w ( x - x , h 坍= 警。 上述中，核函数一般满足以下基本条件： ( 1 ) 正则化条件，即在积分域q 内核函数的积分为1 ，j r v ( x x ，h ) d x = l ； ( 2 ) 当光滑长度h 趋于零时，具有狄克拉函数性质： i i m w ( x 一工，五) = j ( x 一工7 ) ； - - t o ( 3 ) 紧支性条件，在核函数的作用域内， w ( i - x ， ) 为非负值，当 l x - - x 7 l k h ( 一般k 取2 ) 时，w ( x x ，五) = o 。 一般常用的核函数有： ( 1 ) g a u s s i a n 核函数 w ( r ， ) = c j p 一足2 ( 2 1 0 ) 式中d 是空间维数；r 为在点工和x 7 处两粒子之间的相对距离， r ：三：堕，其中，为两点之间的距离；白为归一化常数，在一维、二维 和三维空间中嬲崂土，嘉，玄。 ( 2 ) 样条( s p l i n e ) 核函数 w ( r ，h ) = c d 三c 三r 一，尺2 ，o r ， 丢( 2 一，1 r 2 ( 2 11 ) 0 ，其他 式中，白的值在一维、二维和三维情况下分别为三3 h 乏，1 砌，； ( 3 ) b 样条( b s p i n e ) 函数 硕f j 学位论文 w ( r ，h ) = c d 2 r 2 + ! r 3 o r l 3 2 。 丢( 2 - r ) 3 ， l r 2 0 ，其他 ( 2 1 2 ) 式中，c d 的值在一维、二维和三维情况下分别为_ i 1 1 ，7 砌1 5 ，荔万3 ； ( 4 ) 二次光滑函数 。 矽( r ，| i i ) = c 。( 三一号r + 3 1 6 r 2 ) ，o r 2 ( 2 1 3 ) 式中，白的值在一维、二维和三维情况下分别为i 1 ，孑2 ，i 矛5 。 2 2 2 守恒：h - 程组的离散 s p h 的控制方程基于l a g r a n g e 形式的流体动力学方程 m - - n a v i e r s t o k e s 方程，其中包括了质量守恒方程、动量守恒；b - 程和能 艮施引 亿 式中，v 为速度，e 为能量，仃为应力张量。 在计算过程中要不断地更新应力张量项 仃妒= s 筇一p ( 2 1 5 ) 式中，口和表示张量坐标，p 为压力项，为万函数，为偏应力项 = 2 g ( 锄一三v v ) + 嘞+ ( 2 1 6 ) 式中g 为剪切模量，毒叩为应变率项，为旋转率项，分别为 = 圭c 象+ 薏，= 圭c 薏一薏， c 2 7 ， 利用2 2 1 中的基本方程对控制方程整理，可推导出s p h 方法的离散 墨主岛 = 手l ：- ：( v 支j o - v 一, ) h ，m ：o w o 一 c 2 ，8 ， 一种改进的m l s p h7 j 法理论”究及j e 成用 2 3m l sph 方法的基本理论 m l s p h 方法用移动最小二乘法代替s p h 方法中的核近似法进行插值 计算。这不仅可以提高方法的计算精度，而且能提高方法的稳定性。 2 3 1 移动最小二乘法 假设待求函数“( x ) 在求解域q 中的n 个结点x ，( ，= 1 , 2 ，i v ) 处的函数值 “，= “( x ，) 是已知的，在域q 内构造待求函数“( x ) 的全局近似函数“( x ) ，待求 函数”( x ) 在计算点x 的邻域q 。内可以局部近似为 “6 ( x ，i ) = p l ( x ) a ，( x ) = p r ( i ) 口( x ) ( 2 19 ) 式中i = i xyz 】r 是计算点x 的邻域q ，内任一点的空间坐标， p r ( i ) = b 。( i ) n ( i ) p m ( 2 ) 】，a ( i ) 是基函数， m 是基函数的个数， 口( x ) = 【口。( x ) 口：( x ) a m ( x ) r ，a i ( x ) 是待定系数。基函数a ( i ) 应满足以下条件： p l 【叫2 1 ( 2 2 0 ) 鼽( i ) c k ( q ) ，i = 1 , 2 ，m 、 式中c x ( q ) 表示域q 内具有直到k 阶连续导数的函数空间。通常使用单项式 作为基函数。如二维空间中单项式基函数为 线性基：p r ( i ) = 【l ，x ，y 】，m = 3 二次基：p r ( i ) = 【1 ，茗，y ，x 2 , x y ，y 2 】，m = 6 ( 2 2 1 ) 在移动最小二乘近似( m l s ) 中，系数a i ( x ) 的选取使得近似函数“6 ( x ，i ) 在计算点x 的邻域q 。内是待求函数“( x ) 在某种最小二乘意义下的最佳近 似。计算点x 的邻域q 。称为m l s 近似函数在该计算点处的定义域，简称 为计算点x 的定义域。 将求解域q 用n 个结点离散，在每个结点x = 1 , 2 ，i v ) 处定义一个权函 数( 工) = w ( x 一而) 。权函数彤( 卫) 只在结点而周围的一个有限区域q ，中大于 零，而在该邻域外为零，即该函数是紧支的。区域q ，称为权函数( x ) 的 支撑域，也称为结点x ，的影响域。 设计算点x 的邻域q 。包括个结点，近似函数“( x ，i ) 在这些结点i = x ， 处的误差的加权平方和为： ，= w l ( x ) u 6 ( x ，x ，) - - u ( x ，) 】2 = 形( 工) 易( z ，) 口，( x ) - u ( x ，) 】2 ( 2 2 2 ) 令j 取最小值，即 丽o j= 2 荟( x ) 【蔷b ( x ，h ( 工卜h ，】p ( x ，) = o ，j = l ，2 ，m ( 2 2 3 ) 硕f j 学位论文 由此可得 mn n 【( x ) 易( 西) p ，( x ，) 】口，( x ) = ( x ) 所( _ ， ( 2 2 4 ) f = ii = 1 ，= l 即 彳( 工( 功= b ( x ) u 式中 a - - z , w l ( x