《隐式Newmark法分析波动问题精度的探讨.pdf》由会员分享，可在线阅读，更多相关《隐式Newmark法分析波动问题精度的探讨.pdf（9页珍藏版）》请在三一文库上搜索。
1、第 21 卷第 1 期 1 9 29 年 1 月 爆炸与冲击 PEXL OSIONADNSCH屹KV铲A VSE V o l . 1 2 , J a n . N.ol 1992 隐式 N e wma rk 法分析波动问题精度的探讨 方秦 (南京工程兵工程学院 , 南京2 10007 ) 摘要本文研究 了用平均加速度的 Ne wma rk积分格 式求解一维平面弹性波动间题的精 度 。 由于平均加速度的 N ew ma rk 积分格式没有振幅衰减 、 不产生数值阻尼 , 因此采用此法分析 波动问题会产生较强的 “ 高频振荡 ” , 尤其是加速度的计算 值误差甚大 。 本文对防护工程中常见 的平 台
2、型荷 载和升压三角形荷载作用下的波 动问题 , 应用 有限元分析的实例表明 : 采用有数值 阻尼的 N ew ma r k 积分格式是改 进有限元求解波动问题的一个有效方法 , 可以大大地提高应力 和加速度求解精度 。 文中推荐 的有数值阻尼N e a r k 积分格式的积分参数( 。 . 9 , 声 0 . 49 ) 可供工 程中土与结构动力分析时参考使用 . 关键词有限元法应力波 N ew ma l k法 数值阻尼 在诸如抗核爆炸 、 常规武器爆炸以及抗偶然性爆炸的地下结构等的动力分析中 , 常遇 到瞬态波传播问题 。 求解瞬态波动间题的一种有效方法是特征线法 。 但特征线法只适用于 一些
3、简单的波动问题 , 因 此 , 对于复杂的瞬态波动间题 的求解常常采用有限元法 。 众所周 知 , 有限元法求解波动间题是采用时空域分离的半离散化方法 , 即将波动间题的空域用有 限元离散 , 把波动方程转化为二阶常微分方程 , 然后利用各种方法求得常微分方程的解 。 由此可见 , 由于有限元求解波动间题是采 用半离散方法 , 在空域的单元划分 , 时域的步 长 选取时 , 常常不满足原波动问题的特征线上关系 , 因此必然导致有限元解的误 差 。 文献l 5讨论 了时空域离散产生的各种误差间题 。 在这些文章中 , 一般都以稳态波动 间题为例 进行分析 。 关于瞬态波动 问题的有限元解的精度
4、, 尤其是有限元求解加速度误差问题的讨 论较少 。 有限元求解波动问题存在 “ 低通 滤波 ” 现象 。 即一个单元存在一个载止频率 。豁 , 瞬态 波中的成份高于载止频率将被滤掉 , 只有低于 功思的瞬态波才可通过 。 这是用有限元法求 解波动间题时误差的主要来源之一 。 波动问题与振动问题一个主要区别在于运动的高频 分量对振动问题影响不大 , 即结构可取前几个低 振型反应的叠加 ; 而运动的高频分量对波 动间题影响较大 一般 不能轻易地舍弃阁 。 因此 , 采用有限元法求解波动问题时 , 单元尺寸 和时间步长应足够小 , 以充分反映瞬态波的高频分量的响应 。 为了在运动方程的数值积分 中反
5、映高频分量的影响 , 常常建议采用无振幅衰减的平均加速度的 N e wma rk 积分格式 。 本文以一维平面弹性波动问题为例 , 比较系统地研究了平均加速度的 N ew ma rk 积分 格式求解 波动间题的精 度 。 算例表明 : 平均加速度的 N o wma rk 积分格式没有振 幅衰减 , 不 19 90年1 0 月4日收到原搞 , 199 1 年5月 1 2 日收到修改稿 。 爆炸与冲击第 12卷 能消除由于时空域离散和荷载的间断等产生的应力和加速度的 “ 高频振荡 ” 。 为了消除有 限元解的高频振荡 , 本文利用 Newma rk 积分格式产生的数值阻尼技巧 , 大大地改善了有
6、限元法求解波动问题中应力和加速度的精度 。 一 、 N ew ma rk积分格式的滤波特性分析 波动方程经有限元的空域离散化后 , 成为 M岔+戊十 Kd一F(才) (1) 式中M为质量矩阵 , c 为阻尼矩 阵 , K 为刚度矩阵 , F( t)为荷载矢量 , d 、 泣和注分别为 位移 , 速度和加速度矢量 。 方程(l)常用直接积分法进行求解 , 在直接积分法中使用较多的 是N e wma rk 法 。 ( 1)式的N e wma rk 积分格式可表述为 , 、.、尹 ? 一OJ 了叮、Z毛 M a一+ , + C巧+ , + Kd , +【 = F 一+ , 卫 2 ,. _ _、 .
7、 _ _ 、 + , 己. 一邹 u . 十 二犷至又l 一 艺万)a . 十 Za l l ,+ 1 ) 仇+ 、 一 t t . 十邹(l一 J) a , + 如 .+ : (4) 式中 、 刀为N e wma rk 法积分参数 。 为了研究 N ew ma rk积分格式的稳定性和收验性 , 将(2 ) (4 ) 式改写为 夕 ,+ l 一A夕 , + 丫 。 (5) 式中 、.了 、. 了 内卜甘了Z . f 、 , d :、 价一 、百 A一A护A Z l+ , 2加2 2公,刀 亡 。门 A , 一I _ L戈占。 乙 1+2八td雪 司 (8) 、 产 、 Q UnU 护 了、,.
8、 了 !lw e. ,.e s J 一IJ ;一 半 (:一2刀)。 2 A Z l 一戈(l一的扩 , z一 ,(l一2刀)右 。 1一2戈(l一d)雪。 一!tl.e eJs e j . J . 玉 + 念 ,。 一A 厂 J Z尸 , . 下丁L又I一 2刀)F . 十 2朋 . ,(l一占)F . +研 .+ , 其中 A 为放大矩阵 , 穿 。 为荷载矢量 , 。 为有限元系统的频率 , 右为材料的内阻尼比 。 对方程(5 )进行谱半径分析 , 结果为 , 。 当 。、 、 , , ” 、 ,/2 , ,)(。+ 告 ) 2 /4时 , 积分格式为无条件稳定 ; o 2当o ( 唁蕊
9、 1 , d ) 1/2 , 口 0 . 5时 , 风) A随,/T 增加而减小 , 且 越大 , 冰)A衰减越快 。 这表明 0 . 5 的 N ew ma r k 积分格式对运动的高频 分量起了 “ 阻尼 ” 作用消值越大 , 频率越高 ,“ 阻尼 ” 就越大 。 当 ,/T较小时(比如习/T二 . 0 0 5) , 风)A、 1 . 0 , 即不存在 “ 阻尼 ” 作用 。 这种阻尼并不是材料本身的内阻尼 , 而是数值方 法产生 的 “ 数值阻尼 ” 或称 “ 算法阻尼 ” 。 可以证明 : N e wma rk 积分格式的数值阻尼考可表达 为 1 、,。 二 , 叭 o 一万 ,/乙(
10、11) 由此可见 , 在保证数值积分稳定性条件下 , 要使 N ew ma r k 积分格式产生数值阻尼 , 则必须 使 0 . 5 。 由于 风A)随,/T增大而减小 , , /T太小 , 不产生数值阻尼 , ,/T太大时会产生过大 的数值阻尼 。 根据笔者经验 , 选用 0 . 1( ,/T续1 . 0比较合适 。 考虑到 , / T镇, /几 l。, 其中 T m,。 为有限元系统 中的最小的周期 , 即 T ml。 2汀 m a xo儡 (口) (12) 其中ma x o福为 有限元网格中所有单元 中的单元最高频率 。 对某些形式单元 , 可计算得 (一) 。招 。 例如 , 四结点矩
11、形单元 , 谧 一 粤“ aZ + 。2 aO (13) 式中 c, 为介质压缩波波速 , 。 和 b分别为矩形单元的两个边长 。 根据以上分析 , 参照图 可给出 ,/T的最高值 ,/T ( ,/几 ,。 成 丛 2二/ma x o 4 . x ) _ c, 如 v石厂不丁 汀。b (14) 若取 a 一b , e, ,= a , 则 ,/T成0 . 45 。 以上从理论上分析了 N e wma r k积分格式的滤波特性 。 下面通过典型算例说明有数值 阻尼的 N o wma rk 积分格式改善有限元法求解波动 问题的精度 。 爆炸与 , 冲击第1 2卷 二 、 算例 有限元分析波动问题时常
12、常遇 到时间步 长的选取 , 目前还未给出统一的选取准则 。 文 献6根据荷载函数的升压时间给出步长选取的经验公式如下 、 t/ 升 、 壳平 ) 5 (15) 其中协为荷载的升压时间 , l 为波传播方向有限元网格中的最大单元尺寸 。 (15 ) 式的 经验准则似乎是从有限元法计算应力的精度分析得 出的 , 而对于求解加速度的精度是否 合适 , 值得进一步研究 。 一般认为 , 有 限元法计算加 速度的精度最 差 , 应力和速度次之 , 位 移最好 。 本文采 用文献7推荐的时间步长与单元尺寸的关系式 卫 l/e , (16) 、 、 6 、 、 3 _ 生 甲甲甲 甲甲下 下 9 7 98
13、 . . . LL 9 9 1。奄 奄 0 0 0三 三 1 弹性模量 刀二62 00MPa 2 . 泊桑比 一 0 . 1 5 3 密度 l ,一1 . 637入10 一 ,N S, / em . 4 单元边长 。 一I ,一, ,f t 5 . 波速 e, 20 000em/s 图 2 有限元计算网格及计算参数 Fi g . 2 M e sha nd m aterial Pr oe P rtiesem Ploy ed inF E anal ysis 1 . Elastiemo d u lu s , 2 , P oi s s o n ,5 r atio , 3 . De nsity , 4
14、. Elem e ntsie z , 5 . W a vev elo c lty 图 2 为有限元法的计算网格及其计算参数 。 例 1 平台型荷载作用下情况 : 抓) t一加H( t ) 利用一维波动方程的特征线关系 , 可以求得平台型荷载作用下 一维弹性波动问题的 理论解 位移 u ( x , t) 一尹 。 (t一t。)/烬 , (一7) 速度 。 (: , t) =p o H(z一 t。) /户 e, (一8) 应力 a ( x , t) 一p oH ( ,一z。) (19) 加速度 a ( x , )p。(t一t。) /那 , (20) 式中 H(t)为Ha e v isjd e 函数
15、 , d()为Dsa r e 函数 , t。一x / ( ,。 图 3 为单元L的位移时程曲线 。 从 图 3 可知 平均 加速 度 法和有数值阻尼的 N ew ma rk 积分格式(乃一 0 . 9 , 刀一0 . 4 9或 石一0 . 7 , 刀 0 . 36 ) 均得 到准 确 的 位移计算值 。 由此可以推断 : 位移主要 由运动的低 频分量所组成 第 l 期方秦 : 隐式 New m a k 法分析 波动问题精度的探讨 的 , 荷载函数的间断以及时空域离散等产生的 “ 高频振荡 ” 对其影响不大 , 因此有数值阻尼 和无数值阻尼的 N e wma rk 法均得到准确的位移计算值 。
16、图 4为突加平台型荷载作用下单 元L的应力时程曲线 。 由图 4 可见 , 由于荷载强间断特性 , 所以计算得应力时程曲线出现 了 “ 初始病态 ” 和 “ 超越现象 ” , 其中 , 平均加速度法计算得应力误差为 5 0 纬 , 而有数值阻尼 的算法(取 J=0 . 9 , 刀二0 . 4 9或 d0 . 7 , 声 0 . 3 6 )得到了与理论值基本一致的结果 。 图 5为 突加平台型荷载作用下单元L的加速度时程曲线 , 由该图可见 , 有数值阻尼的算法计算得 到的加速度值与理论值相比较比不 采 用数值阻尼的算法接近 。 这个算例说明了 : 波动有限 元法采用的半离散化法以及荷载本身间断
17、产生解的 “ 高频振荡 ” 可采用数值阻尼进 行消 一 1 . 理论 解 J二 0 . 9 , 声二0 . 4 9 , . 占, 0 . 5 , 声= 0 . 25 一一 J ”0 . 7 , 声二 0 . 36 精 确值 J二0 . 9 d二 0 . 7 10 20 叮“、 ,助 1 . D“” _ I 分冀 护尸户户 。. 1 汗 入入 犬村 汤产一 一 “ 1 1 畔V了 犷一 一 卜 -11一。 。 。5 屏 3 0住0 50 / n w e L s e-盛s e- - 10152 02530 /血 肛介卜 图 3 突加载作用下位移时程 曲线 F琢 3Di sP扭e Cm e nr一r
18、ime h i sto ris e f o rsteP f o ading 1 . Ex actvalue 图 4 突加载作用下应力时程曲线 Fig . 4Str es s一tim e h j so tris e fo r se tP IOai dng 1 . Th eo re t i侧目阳luUo n 1500 10 00 500 0 一 500 1理 论解 乃二 0 . 5 , 声 0 . 2 5 一6一0 . 9 一6=0 . 7 尹二 0 . 49 君0 . 3 6 一1500 例 2 图 5 突加载作用下加速度时程曲线 F jg . 5Ac cele r ai o t n 一tim
19、e hi s 奴i e s f o r s e t P I。di ng1 . Th e份etic als o加tlo n 有升压的三角形荷载作用情况 : 夕(t) t/ t l t , 一 t 亡2 一 之 1 0 ( t 镇 八 忿 : t 镇 tZ (2 1) 爆炸 与 冲击第1 2 卷 由于三角形荷载导数不连续 , 所以其加速度出现间断 。 不难求得理论解为 位移 匹 . (忿一 / c, ) 2 那 , 2亡 , 0( ( tl 夕0 2那 ,止- ( : 一 x / e, ) + P0 C P , (艺 : 一 t:) t Z (z一三一 l) c, (2 2) 一 会 ( , 一
20、: / e, ) , 一 , t, 镇 镇 t: ! ! Y. l I , . e e e e s e l 、 一一 、 .声 念 Z 气 扭 速度 0 P t 卿 ,亡l 0镇 t 续 t l (23) 尹0 脚 , 九一t 才 2 一 t l t l t ( tZ r . . . .悦. 召 性|, 、 一一 、 2 口 刃 ,劣 Z、 沙 应力 a (劣 , t) = 尹。t/t、 tZ 一 忿 p“ 而二万 ; 0镇 名蕊tl t : t ( tZ (24) 加速度 !业 l那声 、 a Lx , t) =凡 一 . 一一卫一一 - t 那 , ( t: 一 乙 1 ) 0镇 t 攫
21、t , (2 5) , 成 2 计算时取 t;= l。 , tZ=l0 00ms 。 图 6为三角形荷载作用下单元L的应力时程曲线 , 由图 6可知 , 当 , =0 . lm s 时(满 足( 15)式 , 无论N e wma rk积分参数取d=0 . 5 , 声= 0 . 25 或 J二0 . 7 , 刀 一 0 . 36 或 J=0 . 9 , 夕二 0 . 4 9 均得到与理论值一致的应力计算结果 。 图 7 为三角型荷载作用下单元L加速度时程 曲线 。 由此可见 , 即使平均加速度法计算得应力值是正确的 , 但计算得加速度却增大了 33% , 而有数值阻尼的算法得到了与理论值基本一致
22、的结果 。 取积分参数二 0 . 7 , 声二 0 . 3 6 , 算得加速度峰值增大了 9 . 3写 , 而 d二0 . 9 , 刀=0 . 4 9算得加速度峰值与理论值一致 。 图 8为取, =0 . Zms(即, t/ 开二1/5 ) 计算得单元L的应力时程曲线 , 当 ,/妞一l5 / 时 , 用平均加速度法计算得应力峰值误差为2 0% , 而有数值阻尼的算法(取积分参数 占 0 . 9 , 声= 0 . 49或d=0 . 7 , 刀二0 . 36均得到与理论值一致的结果 。 图 9 为 ,二Zm s (即,/之 升 二1/5 ) 时的单元L加速度时程曲线 , 由图 9可见 , 平均加
23、速度法计算得加速度值出现 了高 频振荡 , 峰值误差高达 12 5 % , 0 . 7 , 声 0 . 3 6计算得加速度峰值误差为 3 3 % , 而 d一. 0 9 , 刀 0 . 4 9得到的加速度峰值与理论值基本一致的结果 。 图 6 、 图 7 与图 8 、 图 9相比 , 当 时间步长减小一倍时 , 用平均加速度法得到的应力计算值虽有一定的改善 , 但其加速度计 算值误差仍较大 。 算例表明8 : 通过减小时间步长来提高加速度计算精度是不可取的 , 也 不能得到与理论值一致的加速度计算值 。 第 1 期方秦 : 隐式 N ew ma rk 法分析波动问题精度的探讨15 1 00 _
24、 _占二0 . 5 , 声一 0 . 2 5; J一 0 、 9 , 刀二 0 . Jg ; J二0 . 7 , 声 0 . 3 6 a /Ml场 1理论值 一100 0 . 1 0 一 2 0 0 0 . 05 一 3 00 1 0 2 0 3040 忿/jr 图 6 三角形荷载作用下应力时程 曲线 (A亡0 . lm s ) F返 . 6St r e s s 一 timeh 妇to r i韶fo r t r认ngu扭r Ia Odingw irh r如e tin r e (配=0 . l m s ) 1 . Th e ot改i c81valu必 一J o o 十 、 图 7 三角形荷载作用
25、下加速度 时程曲线(配= 0 . l ms) F论 . 7Ac份拓ratio n一itm e histo r油 l o r t ria ngula:la oi d n g 州h t r 娜 time(配=0 . l m s ) 1 . The o retic a l value 0 0 0 0 0 00 nj, .1 口/ 、!下场) 任 一10 0 一200 理论值 一30 0 一 5 00 沙 0 . 5 , 尸一0 . 25 6二0 . 9 , 声一0 . 49; d0 . 7 ,尹一 0 . 36 一 . t0 0 z 搜矛 0 . 0 5 0 . 0 J 0 203Ot/亡 图 8
26、三 角形 荷载作用下应力时程曲线(从= 0 . Zm s) Fig . 8Stres s一tim e h i sto rie s fo r triangu坛rlo ai d ng 侧th r俪time (却= 0 . Zm s ) 1 . T he o re tic al s olution 三 角形荷载作用下加速度时程 曲线 (超 0 . Zms) Aeeele了atio n 一tim e h妇to rie sfo r tria ngu加r 10 a山 ng wih t rl加 i tm e (川=0 . Zms) 1 . T he o rei t c al v alue 爆炸与冲击第2 1
27、 卷 三 、 结论 通过以上理论分析及算例 , 可得到以下几点看法 : 1 . 有限元法分析波动间题时 , 无论采用有无数值阻尼的协w ma rk 法均能得到准确的 位移计算值 。 2 . 当满足一定条件时 , 如式( 15 ) 平均加速度法能得到可靠的应力计算值 , 但加速度计 算值误差较大 。 误差随荷载本身间断性的强弱 , 时间步长的大小等而变化 。 荷载本身间断 性越强 , 时间步长越大 , 误差越大 ; 反之亦然 。 固此 , 平均加速度法不宜用于波动间题有限 元分析 . 3 , 采用有数值阻尼的N e wma rk 积分 格式时 , 无论是在突加平台型荷载 , 还是有升压 三角形荷
28、载作用下均得到与理论值基本一致的应力和加速度计算结果 。 建议在波动问题 有限元分析时采用积分参数为 d 0 . 9 , 声二0 . 钧 的决w ma rk 积分格式 。 4 . 有限元法分析波动问题时 , 由于荷载间断 , 时空域离散等所产生解的 “ 高频振荡 ” 可 田数值阻尼进行消除 。 固此较其它办法s 简单方便 。 5 . 应用数值阻尼的 N ew m, k法分析波动问题时 , 对于有升压的荷载 , 选取时 间步 长 可放宽为 ,/ ,升 音 。 本文的结论虽在一维平面弹性波情况下得出的 , 但其结论可供二 、 三维波动问题有限 元法分析时参考使用 。 由于篇幅所限 , 关于中心差分
29、法 、 w ”s o n 一。法 , 显式 N ew ma rk 法以 及二维相互作用分析另文讨论 。 参考文献 决Iy拓 hkoT , H.助.TJR . C帅pua t ti叨目 M .th司, f倪T t a耐c ntA刀aly由: Els e吹rSd e们。e u P b l台h翻活 , 1983 廖振腆 , 刘晶波 . 地震工程与工程振动 , 1 98 6 , 6(2) : l一6 刘晶波 , 序振肠 . 地震工程与工程振动 , 1 9 8 9 , 9(2) , 1一 1 1 刘晶波 , 序振鹅 . 地展工程与工程振动 , 1 99 0 , 1 0 2 ) : 19 ValliaP伸
30、ns ,人ng K K . Dy nal n 沁人刀.1邓肠人pP】 i司o tRO CkM e C汕”i。巧曲全 m s sth功tCo nf m Nu砒Me th , 0加n 斌尤h, 1 985 : 1 1 91 3 4 . 孙钧 , 吴逸群 . 岩土工程学报 , i o ax , 5(4) : 2 8一 2 伪t址 KJ . l F nite Ec Ime ntr PO Cd e ure s inEngine七nns 人naly雍 . 阶w 介 n 把y : Pr七nt晚一枷1ll n c E ng jew心 a f枯 , 19 82 方秦 , 显式N e Ml ln a k t 法求
31、解波动问题精 度探讨(待发表) ,.f J.J, J,.J工.J ,J,J.,J ,占,臼 九J怪翻台启O , .且 甘r.lL .J L .L p胜 Jr .L J.L J.J L.J S TUDIESONTHE ANAL Y S ISO F FOR W AVE ACC U R ACYO F I F NI TE EL EMENT IMP L ICI TNE专VMARK ME THOD PRO PA G AT I ONPR OBL EMS a F ng Qin (孙咖户,瓜扣田喃夕爪对秘 曰他,Ar邓咨 已帅产 袱肠户曲 习, N如乒, 2 10 0 0 7) AB S T RACTThis
32、a Pp e rstudis e the a c e u r a eyo f usjng the av era g e ae“ la e r tjo n m ethd oo fN e w - m a rk,5jntegr a l n er in s o lvjng o n e一 dim e n sio nale las t iew av e Pro a P助ti on prob lem s . e Bc a u s e the 第 1 期 方秦 : 隐式 N e w m ok法 分析 波动问题精度的探讨 5 3 avera 解 a c c e ler a t i o nem tho d h a
33、s no ca阶b iliti eos f a m P l j tu d e d ec a y an d nu m e ric a lda ”nPin g , h t e st rongoc川a t o r y re sPo ns emay b e r e s ulte d whe nusing this m etho d to a n a】 ys ew a鸭 t P o伸醉tio n r P ob】e m s , e se p eja l lya r即er r o r m a y b e intro d uee din e omPu a t tion alac c elerai t o n
34、. n I a d 山tl助 , h tis 伴p e r s t u di e s h t eela s t iew ave Pr o 件梦ti on Problems u nde rste Plo ading a n d tri助即la r 】 创川加g , bothwith ris e ti们。 e s , u singnu m erieal da m i P n g te ehnjque s . N u m e riealex a mPlo h s ow t抽t h ti s te Chniq讹15 a nef介etive w a J , to im讲 o v e th e a c c
35、 uraeyof F E a n aly幻5 fo r a w v e PrOa P 醉do nPO r卜 lem s , 夕e atly in e r创招ing th e ae eur aeyofeom 四 tatio n alstr侧沼 a n d a c沈lea r ti叻 . N ew m肚k nletho d wih t th e i f le t ring 讲 o e p rtyof 阵 ra jr neter v a lue sJ 0 . 9 , 声” 0 . 49jss ug ge s te dto 加 u创月 i n com - P】ieate d 5011一stru ctu re inter a etiona n a ly sis K E Y WO R D S finite ele m ent m etho d , ste r 貂w a ve , N ew m ark m etho d , n u m ei r a C l da t n Ping