基于增量本构关系弹塑性分析的无网格伽辽金法.pdf

第3 4 卷第4 期 2 0 0 5 年7 月 中国矿业大学学报 J o u r n a lo fC h i n aU n i v e r s i t yo fM i n i n g ＆T e c h n o l o g y 文章编号1 0 0 0 1 9 6 4 2 0 0 5 0 4 0 5 0 9 0 5 基于增量本构关系弹塑性分析 的无网格伽辽金法 赵光明1 ’2 ，宋顺成2 ，杨显杰2 1 ．安徽理工大学资源开发与管理工程系，安徽淮南2 3 2 0 0 1 } 2 ．西南交通大学应用力学与工程系，四川成都6 1 0 0 3 1 V 0 1 ．3 4N O ．4 J u l ．2 0 0 5 摘要在无网格伽辽金法的基础上，采用增量形式的无网格伽辽金法的插值方法，并利用应力应 变增量形式表征材料的弹塑性本构关系，在小变形假设的前提下，提出了基于增量本构关系的弹 塑性分析的无网格伽辽金法；采用罚参数修正了能量变分方程式，方便地实现了无网格伽辽金法 的本质边界条件．对二维平面问题及轴对称问题的计算结果表明运用增量形式的无网格伽辽金 法进行弹塑性问题的数值分析，具有较高的计算精度，且适应性强，从而为材料的弹塑性分析计 算提供了一种实用方法． 关键词无网格伽辽金法；弹塑性；本构关系；移动最小二乘 中图分类号02 4 1 ．8 2文献标识码A E l e m e n t ．．F r e eG a l e r k i nM e t h o df o rE l a s t o ．．P l a s t i c A n a l y s i sB a s e do nI n c r e m e n t a lC o n s t i t u t i v eE q u a t i o n s Z H A OG u a n g m i n 9 1 ”，S O N GS h u n c h e n 9 2 ，Y A N GX i a n j i e 2 1 ．D e p a r t m e n to fR e s o u r c eE x p l o r a t i o na n dM a n a g e m e n tE n g i n e e r i n g 。A n h u iU n i v e r s i t yo f S c i e n c ea n dT e c h n o l o g y ，H u a i n a n ，A n h u i2 3 2 0 0 1 ，C h i n a ；2 ．D e p a r t m e n to fA p p l i e dM e c h a n i c s a n dE n g i n e e r i n g ，S o u t h w e s tJ i a o t o n gU n i v e r s i t y ，C h e n g d u ，S i c h u a n6 1 0 0 3 1 ，C h i n a A b s t r a c t O nt h eb a s i so ft h ee l e m e n t f r e eG a l e r k i nm e t h o d E F G M ，t h ee l e m e n t f r e eG a l e r k i n m e t h o df o re l a s t o p l a s t i ca n a l y s e sw a sp r o p o s e db a s e do nt h ei n c r e m e n t a lc o n s t i t u t i v ee q u a t i o n sa n d t h es m a l ld e f o r m a t i o nh y p o t h e s i s ．I nt h i sp r o p o s e dm e t h o d ，t h ei n c r e m e n t a li n t e r p o l a t i o nm e t h o di n t h eE F G Mi sa d o p t e d ，t h ei n c r e m e n t so fs t r e s sa n ds t r a i na r eu s e dt oc h a r a c t e r i z et h ee l a s t o p l a s t i c c o n s t i t u t i v er e l a t i o n s h i p ，a n dt h ep e n a l t yp a r a m e t e r sa r ea d o p t e dt or e v i s et h ee n e r g yv a r i a t i o n e q u a t i o n st oe a s i l yr e a l i z et h ee s s e n t i a lb o u n d a r yc o n d i t i o n si nt h eE F G M ．T h ec a l c u l a t i n gr e s u l t so f 2 Dp l a n ea n da x i s y m m e t r i c a lp r o b l e m ss h o wt h a tn u m e r i c a la n a l y s i sf o re l a s t o p l a s t i cp r o b l e m sb yo f t h ei n c r e m e n t a lE F G Mh a sa d v a n t a g e ss u c ha sh i g hc a l c u l a t i o na c c u r a c ya n ds t r o n ga d a p t a b i l i t y ． T h e r e f o r e ，t h i sm e t h o dc a ns o l v et h ec a l c u l a t i o np r o b l e m so fm a t e r i a le l a s t o p l a s t i ca n a l y s i s ． K e yw o r d s e l e m e n t - f r e eG a l e r k i nm e t h o d ；e l a s t o p l a s t i c i t y ；c o n s t i t u t i v ee q u a t i o n ；m o v i n g l e a s t s q u a r e s 近几年，无网格法已成为计算力学领域里十分 热门的研究课题，该方法的基本要点是不需要网格 单元，只需要具体的结点信息，直接利用结点信息 对计算域内任一点的位移进行拟合，这样可以消除 收稿日期2 0 0 4 1 0 一2 0 基金项目国家自然科学基金项目 1 0 3 7 2 0 8 6 作者简介赵光明 1 9 7 6 一 ，男，安徽省桐城市人，讲师，博士研究生，从事高速冲击动力学和数值模拟方面的研究 万方数据 5 1 0中国矿业大学学报 第3 4 卷 或部分消除 网格存在带来的困难，因此，在处理 非线性问题、不连续问题，尤其滑移面不连续问 题1 1 - 5 3 时，具有明显的优点． 无网格法起源于L ．B ．L u c y 在解决无边界天 体物理时提出的光滑质点流体动力学方法 S P H ， 而后，N a y r o l e s [ 6 3 等人提出了扩散单元法 D E M ， B e l y t s c h k o c T J 在其基础上，提出了无网格伽辽金法 E l e m e n t f r e eG a l e r k i nm e t h o d 。E F G ． 本文在文献[ 8 9 ] 的基础上，通过引入新的能 量泛函式，并将增量形式的应力应变关系代入泛函 式中，推导出基于增量本构关系的无网格伽辽金法 计算时的离散控制方程，这种计算方法能方便地适 应弹塑性行为与加载以及变形过程有关的特性． 1无网格伽辽金法的基本原理 在域力内，任意一点z 的位移Ⅳ z 近似函数 可由移动最小二乘原理构造为 H 工 ≈U h x 一∑P f x 口f x 三 i f f i l P T 工 口 z ， 1 式中P ， z 为任意阶的基函数；m 为基函数的次 数；a x 为系数矩阵．采用移动最d x - “ 乘技术，构 造加权范数J 工 ，并使得其取最小值． J x P 口一H 。 T W 工 P 口一Ⅳ。 ． 2 由最小二乘原理必须有 i d d A 工 口 工 一H x u 。一O ， 3 d 口 其中A P T W x P ， 日 P T W 工 。 P W I 户1 z 1 户2 z 1 ⋯ I P l z z P 2 z z ⋯ L 。0 。，户。0 。，⋯ 呦 x x 1 O l 以z ■ loo 户。 z 1 ] 户。 z 。 I i P 。 z 。 J ⋯ O ⋯ O ⋯ W x - x 。 式中w x z i 为点X i 在工处的权函数；砧。为名义 结点值．由此，可以得到 a x A 一1 工 H x Ⅱ。． 4 将式 4 代入式 1 ，有 U “ x P T 工 A 一1 x 日 x H 。 ∞ x H 。， 5 其中 称为形函数． 显而易见，形函数受到权函数和基函数的影 响，只要选定适当的权函数和基函数便可以对问题 进行求解．基函数一般选择一阶或二阶，权函数的 选择对无网格伽辽金法的贡献非常大，它直接影响 到问题的求解，文献[ 1 7 ] 中提出了很多的权函数， 如幂权函数、．指数权函数、样条权函数。其中本文采 用了较为常用的幂权函数 w d i fexp[--d,／c2“]-expE-dm．x／cz-]di≤厶J，1-ex d 。。 ， 式中d 。 l | x 一麓l | 为两点间距离 d m a 。为权函数支 持域的最大半径；f 为控制常量，决定权函数的形 状． 2 弹塑性增量无网格伽辽金法 弹塑性问题的分析通常与加载及变形过程有 关，因此在计算中通常将载荷分成若干增量，在每 一增量步中按线性问题来处理．假设任意两连续的 增量步时刻为t 和f ，在小变形的弹塑性分析 中，除应力应变关系以外，其它方程和边界条件都 是线性的．假设在t 时刻的条件已知，即体载荷 ‘F ，，边界载荷‘_ l ，位移条件‘面已知．当过渡到t 时，载荷和位移条件有一增量，即 件血F i ‘F f △F f ， 7 a 件m T f ‘T i △7 ’f ， 7 b 件血阮 ‘U i A u f ． 7 c 则位移增量血和应变增量△E 需满足能量泛函式 ／／ I , a T A e 百1m T D e P A e d a Q ～ 乞 厶 I A u T ‘F A F d , O I △矿 t 亍 △亍 d P ， 8 8 r d 式中叮，E 分别为应力向量和应变向量；左上标t 表 示时刻；D e 。表示弹塑性矩阵；△表示增量．将无网 格伽辽金法的插值方法表示成增量形式 ， 纽一I z I 血。， 9 a A e B A u 。，。 9 b 式中形函数矩阵 ．r j 5 ，0 声z 0 ⋯丸0 ] 归【．o 庐，o 声。⋯o 丸J ’ p 1 ．。0 声2 ．。0 ⋯ 丸．。0 ] 曰 1 0 声l 。，0声2 ．，⋯0 丸．，1 ． b ，。， 奴，。 丸，， 九。。 ⋯ 丸，， 丸．。j t D x p T x A 一1 x 日 x 6 将式 9 a ， 9 b 代入式 8 得 万方数据 第4 期赵光明等基于增量本构关系弹塑性分析的无网格伽辽金法 5 1 1 Ⅱ一l A u 。T B T t 叮 丢△l I 。T B T D e p B A u 。 d O J厶 l A u 订∞T ‘F A F d O 一 尝 I △Ⅱ’T ①T ‘亍 △亍 d P ． 1 0 P 4 其变分式为 泗 I 艿△l l 订B T ‘叮 艿衄订B T D 。p B L X u d O 一 矗 1 8 血订∞T ‘F A F d O 一 尝 I8 △髓。T ∞T ‘亍 A T d P ． 1 1 r d 对于上式的变分式还必须满足本质边界条件， 我们知道，通过移动最小二乘拟合出的形函数，不 具备K r o n e c h e rd e l t a 的性质，即办 工』 ≠如，本质边 界n 上H 一五或A u △五的条件不易实现，本文通 过罚因子口来实现本质边界条件，在变分式中引入 能量增量函数 A ／- ／ 口IA u T A u Z 正 d P ． 1 2 ≠ 则式 9 a 引入式 1 2 的变分式可得 艿△1 I 2 al6 △口。T O r S O A u 。d r L 口I6 △Ⅳ订∞T 瓜d r ， 1 3 L 式中s [ 詈兰] ，如果本质边界点上给定z 或 y 方向位移，则S 。 1 或S ， 1 ，否则S 。 0 或 S ， O ． 则修正后的能量变分式为 令8 Ⅱ一0 可得 K o p 置。 z a u ’ Q F 十Q T Q 。一Q 。， 1 5 简记为K A u ’ a Q ， 1 6 其中K e p IB T D e 。B d J 2 ；K 。 2 口I ∞r S O d F ； 乞 屯 Q F l ∞T ‘F A F d O ； 乞 Q T I 矿 ‘亍 △亍 d 力； 。 Q 。一a l ①T J S △_ d r ；Q 。 l B r ‘叮 d O ． 屯 南 3求解方案 3 ．1迭代公式 对于弹塑性增量问题的分析，本文采用牛顿一 拉斐逊法，可将式 1 6 改写为 ‘ 趾眉‘”’△脚。‘”’ a Q ‘⋯， 1 7 其中件“K “’一I B T 件血D 窜’ B d O K 。， 当 a Q I 矿 件m F ‘n ’ d n I 矿 ‘ 血亍 “ ’ d r 乞 ≠。 o j 矿I S 氙d r j B T 件吲” d O r “ 仃 3 ．2 弹塑性矩阵D e 。的计算方法 假设弹性和塑性矩阵分别记为鼠和D p ，则弹 塑性矩阵D 。，一D 。一D ，．以平面应力为例，理想塑性 材料的屈服函数F i 1S i i s o - - 号商，s 一[ 趴 娩 轧] T 为应力偏量向量，仉。为屈服应力，由塑性流动 法则可以推导塑性矩阵 r s 1 1 芦s 2 2 2 s 1 1 p s 2 2 s 2 2 芦s 1 1 D p - 南I 也z 怕1 2 l 对称 式中E ，卢分别为弹性模量和泊松比； G s } 1 s i 2 2 即1 1 s 2 2 2 1 一卢 s ；2 ． 3 ．3 本构关系积分 在对任一高斯积分点的计算时，容易求得该点 得血和△E ，对于应力增量计算一般采用积分式 △叮 ID e ，d E 计算，可通过弹性因子法E 们来计算． 为了提高计算精度，需要对高斯积分点的结果通过 切向预测径向返回子增量法进行修正． 4算例 4 ．1受均布载荷的悬臂梁 图1 为悬臂梁受到均布载荷g 的作用，梁的基 本参数L 8m ，H 3m ，泊松比卢 0 ．3 ，弹性模 量E 1 9 6G P a ，屈服应力以。一2 1 6M P a ，材料为理 想塑性．循环加载过程如表1 ．无网格法计算中，划 分为8 3 的积分网格，每一网格中采用3 3 高斯 点，d 。。一41 T I ．通过不等步长的载荷增量步来实现 外加载荷．图2 绘制了载荷q 1 7M N ／m 梁的中 心线挠度图；图3 是由g 0N ／m 到q - - 1 7M N ／m 8l，L 1，●●●●●●●j n n 、，＼， 舰胁也 户 n 毖 ∥ 0 0 一 1 A 段 一 一 1 1，＼，＼ 万方数据 5 1 2 中国矿业大学学报 第3 4 卷 A 点位移与外加载荷关系；图4 是循环加载下的B 点位移图．本文也将A N S Y S 软件计算的结果绘制 在对应的图表中，无网格伽辽金法与A N S Y S 软件 结果比较，图2 ，3 中的最大的相差分别是4 ．2 ％， 4 ．3 6 ％，说明无网格伽辽金法解决弹塑性问题是可 行的． 表1 循环加载路径 T a b l e1 P a s s w a yo fc i r c u l a rl o a d s 时间／s 载荷q ／ M N m 1 Y q 图1 受均布载荷的悬臂梁 图2 沿梁中心线的挠度 F i g ．2D i s p l a c e m e n ta l o n gt h em i d d l el i n e l ＼ 稔 迥 图3A 点Y 方向位移与载荷g 的关系 F i g ．3 R e l a t i o nb e t w e e nd i s p l a c e m e n ta n dp r e s s u r e o np o i n tAi ny - d i r e c t i o n 2 0 一1 5 昌1 0 主5 皇0 弼 ．2 0 位移／m ” 图4 循环加载下B 点的位移 F i g ．4 D i s p l a c e m e n to np o i n tBu n d e rc i r c u l a rp r e s s u r e 4 ．2 受内压圆筒 图5 是一轴向受约束的承受内压的厚壁圆筒， 材料为理想塑性，服从y o n ．M i s e s 屈服条件．圆筒 的基本参数口 1c m ，b 2a m ，E 6 5G P a ，∥ 0 ．4 ，O “ s o 一1 7 ．3 2M P a ，户 1 2 ．5M P a ．弹塑性理论解 采用文献[ 1 0 ] 的公式计算，图6 表示径向位移的计 算结果，无网格伽辽金法的计算结果与理论解极为 接近．图7 是应力分布的计算结果，实际结果表明 在同等网格 结点 布置情况下计算出的弹性与塑 性分界面的半径，无网格方法比A N S Y S 软件与理 论解结果更为接近． 图5 受内压的圆筒 F i g ．5At h i c kt u b eu n d e ri n t e r n a lp r e s s u r e x ／m m 图6 径向位移分布 F i g ．6 D i s t r i b u t i o no fd i s p l a c e m e n ti nr a d i a ld i r e c t i o n 理论解 芒 ◆无网格解 三 A N S Y S 翩t b 图7 沿壁厚方向的应力分布 F i g ．7 S t r e s sd i s t r i b u t i o ni nt h et h i c k n e s sd i r e c t i o n O 8 6 4 2 O 8 6 万方数据 第4 期赵光明等基于增量本构关系弹塑性分析的无网格伽辽金法 5 1 3 5结论 本文探索如何用一种新的数值方法一无网格 伽辽金法来求解弹塑性问题，通过引入新的能量泛 函式和罚因子能量增量式，推导计算弹塑性问题的 无网格伽辽金法控制方程．从数值算例的计算结果 来看，无网格方法具有较高的计算精度，与 A N S Y S 软件计算出的结果相差甚少，说明在弹塑 性数值计算中是可行的．如果将这种方法推广应用 到弹塑性的不连续问题，将具有有限元法不可比拟 的优越性． 参考文献 I - 1 - 1 B e l y t s c h k oT ，K r o n g a u zY ，O r g a nD ．M e s h l e s sm e t h o d A no v e r v i e wa n dr e c e n td e v e l o p m e n t s [ J ] ． C o m p u tM e t h o d sA p p lM e c hE n g r g ，1 9 9 6 1 3 9 3 4 7 ． 1 - 2 3 F l e m i n gM ，C h uYA ，O r g a nB ．E n r i c h e de l e m e n t f r e eG a l e r k i nm e t h o d sf o rc r a c k t i pf i e l d [ J ] ．I n tJ N u m e rM e t hE n g n g ，1 9 9 7 4 0 1 4 8 3 - 1 5 0 4 ． [ 3 ] B e l y t s c h k oT ，O r g a nD ，G e r l a n c hC ．E l e m e n t f r e e G a l e r k i nm e t h o d sf o rd y n a m i cf r a c t u r ei nc o n c r e t e [ J ] ．C o m p uA p p lM e c hE n g r g ，2 0 0 0 1 8 7 3 8 5 3 9 9 ． [ 4 ] L i e wKM ，N gTY ，W uYC ．M e s hf r e em e t h o df o r l a r g e d e f o r m a t i o n a n a l y s i s - ar e p r o d u c i n g k e r n e l p a r t i c l ea p p r o a c h [ J ] ．E n g i n e e r i n gS t r u c t u r e s ，2 0 0 2 2 4 5 4 3 5 5 1 ． [ 5 ] H a g i h a r aS ，T s u n o r iM ，I k e d aT ．E l e m e n t f r e eG a l e r k i nm e t h du s i n gd i r e c t e dg r a p ha n di t sa p p l i c a t i o nt O c r e e pp r o b l e m [ J ] ．C o m p u t a t i o n a lM e c h a n i c s ，2 0 0 3 3 1 4 8 9 4 9 5 ． [ 6 ] N a y r o l e sB ，T o u z o tG ，V i l l o nP ．G e n e r a l i n gt h ef i n i t e e l e m e n tm e t h o d ．d i f f u s ea p p r o x i m a t i o na n dt h ed i f f u s e e l e m e n t [ J ] ．C o m p u t a t i o n a lM e c h a n i c s ，1 9 9 2 1 0 3 0 7 3 1 8 ． [ 7 1B e l y t s c h k oT ，L uYY ，G uL ．E l e m e n t f r e eG a l e r k i n m e t h o d [ J ] ．I n tJf o rN u m e rM e t hE n g n g ，1 9 9 4 3 7 2 2 9 2 5 6 ． [ 8 ] 丁皓江，何福保．弹性和塑性力学中的有限单元法 [ M ] ．北京机械工业出版社，1 9 8 9 ． [ 9 ] 王勖成．有限单元法[ M ] ．北京清华大学出版社， 2 0 0 3 ． [ 1 0 ] 严宗达．塑性力学[ M ] ．天津天津大学出版社， 】9 8 8 ． 责任编辑邓群 万方数据