第2章 显式分析理论

7 第二章第二章第二章第二章 显式时间积分显式时间积分显式时间积分显式时间积分 如图所示，先考虑简单的单自由度线性弹簧阻尼系统，根据达朗贝尔动力学原理可得： )(tpkuucum =++ &&& （2.1） u&& 为加速度，u&为速度，u为位移， )(tp 为外力。 大家知道，对于该线性问题，可以用解析方法来求解该常微分方程。若为非线性问题， 比如 k为位移 u的函数，方程为： )()( tpuukucum =++ &&& （2.2） 此时用解析方法一般很难求解，所以应用数值解法来求解，常用的为有限差分法和有限 元法。上述公式具有普遍意义，对于有限元法而言，上述运动方程的矩阵形式为： )(tPKUUCUM =++ &&& （2.3） 式中U&&为节点加速度列阵，U&为节点速度列阵，U 为位移列阵， )(tP 为外力向量列阵，M 为质量矩阵，C为阻尼矩阵，K为刚度矩阵。 求解该运动方程目前有两种方法用的较多，一种是振型叠加法，一种是逐步积分法，对 于复杂问题，一般采用逐步积分法，大体可分为增量法，迭代法和混合法。 隐式的求解方法一般是采用增量迭代法，需要转置换刚度矩阵，通过一系列线性逼近 (Newton-Raphson) 来获得解，对于存在内部接触这样的高度非线性动力学问题往往无法保 证收敛。LS-DYNA采用显式中心差分法来进行时间积分，在已知 0，……， nt 时间步解的 情况下，求解 1+nt 时间步的解，运动方程为： )()()()()( int nnnnn tUCtHtFtPtUM &&& −+−= （2.4） 8 式中 )(tP 为外力向量列阵， )(int ntF 为内力矢量，为单元内力和接触力之和，表达式 为 ∑∫ +Ω= Ω contact n T n FdBtF σ)(int ，单元的内力由当前构型的应力场的散度求得， )( ntH 为沙漏阻力。 把质量矩阵移到方程的右边，求得 nt 时刻的加速度为： [ ])()()()()( int1 nnnnn tUCtHtFtPMtU &&& −+−= − （2.5） 1+nt 时刻的速度和位移由下面公式求得： nnnn ttUtUtU ∆+= −+ )()()( 2 1 2 1 &&&& （2.6） 2 1 2 11 )()()( +++ ∆+= nnnn ttUtUtU & ，其中 2 1 2 1 + + ∆+∆ =∆ nn n tt t （2.7） 这样可以求得在 1+nt 时刻的位移，更新 nt 时刻的系统几何构型，得到 1+nt 时刻的系统新的几 何构型。由于采用集中质量矩阵M ，运动方程的求解是非耦合的，不需要组集成总体刚度 矩阵，并采用中心单点积分，因此大大节省存储空间和求解机时。 但是显式中心差分法是有条件稳定的，可以通过一个简单的线性自由弹簧系统来进行说 明，此时运动方程为： 0=+ KUUM && （2.8） 设φ为特征向量矩阵，则： 0=+ UKUM TT φφφφ && （2.9） 由于 IMT =φφ ， 2ωφφ =KT ，ω为圆频率，于是 nt 时刻运动方程为： 0)()( 2 =+ nn tUtU ω&& （2.10） 如果时间积分采用中心差分法，那么： t tUtU tU nnn ∆ − = −+ 2 )()( 11 ）（ & （2.11） 2 11 )()(2)()( t tUtUtU tU nnnn ∆ +− = −+&& （2.12） 其中 t∆ 为时间步。把 )( ntU&& 代入运动方程，可得： ( ) 0)()(2)( 1221 =+∆−− −+ nnn tUtUttU ω （2.13） 设 n ntU λ=)( ，代入方程就可把差分方程变为多项式方程： 9 ( ) 012 222 =+∆−− λωλ t （2.14） 当 ∞→n 时，若 )( ntU 是有界的，则可以得到稳定解，这就要求 1≤λ ，对于 1≤λ 的方 程所有的根，满足稳定条件的最大 t∆ 的值为： max 2 ω =∆ critt （2.15） 式中 maxω 为有限元网格的最大自然角频率。所以只有当 ∆ ∆t t crit≤ = 2 ω max （2.16） 此时，求解才是稳定的，所以显式算法采用很小的时间步来进行计算，一般只对瞬态问题有 效。 大家注意到，在运动方程 )()()()()( int nnnnn tUCtHtFtPtUM &&& −+−= 中有一项沙漏 阻力 )( ntH ，该项是人为加上的力，是为了防止沙漏变形，这是由于 LS-DYNA采用单点高 斯积分进行时间积分，会出现沙漏模式（零能模式），下面简单描述沙漏模式产生的原因。 单元内力由当前构型应力场的散度求得： ∑∫ Ω Ω= dBtF nTn σ)(int 单元计算 ∫ e T dVB σ 时，应力增量 t∆σ& 由应变率ε&根据材料本构关系求出，应变率ε&与 单元速度场 1x&， 2x&， 3x&有关，对于 8节点六面体实体等参元来说，单元内任意点的速度分 量为： ( ) ( ) )(,,,,, 8 1 txtx ki k ki && ∑ = = ζηξφζηξ （2.17） 式中形函数为： ( ) ( )( )( ) ( )ξηζζηξξζζξηζζηξηηξζζηηξξ ζζηηξξζηξφ kkkkkkkkkkkk kkkk +++++++= +++= 1 8 1 111 8 1 ,, ( )kkk ζηξ ,, ， 8,2,1 Κ=k 是节点的自然坐标值，见表 2.1 节点 1 2 3 4 5 6 7 8 ξ -1 1 1 -1 -1 1 1 -1 η -1 -1 1 1 -1 -1 1 1 ζ -1 -1 -1 -1 1 1 1 1 10 上式用矩阵形式表达为： ( ) ( ) ( ){ }txtx kiTTTTTTTTi && ξηζζξηζξηζηξζηξ 432132181,,, Γ+Γ+Γ+Γ+Λ+Λ+Λ+∑= （2.18） 式中： [ ]T11111111=∑ [ ] [ ]TT 11111111876543211 −−−−==Λ ξξξξξξξξ [ ] [ ]TT 11111111876543212 −−−−==Λ ηηηηηηηη [ ] [ ]TT 11111111876543213 −−−−==Λ ζζζζζζζζ [ ] [ ]T T 11111111 88776655443322111 −−−−= =Γ ηξηξηξηξηξηξηξηξ [ ] [ ]T T 11111111 88776655443322112 −−−−= =Γ ζηζηζηζηζηζηζηζη [ ] [ ]T T 11111111 88776655443322113 −−−−= =Γ ξζξζξζξζξζξζξζξζ [ ] [ ]T T 11111111 8887776665554443332221114 −−−−= =Γ ζηξζηξζηξζηξζηξζηξζηξζηξ 单元的速度场是由基矢量 321321 ,,,,,, ΓΓΓΛΛΛ∑ 和 4Γ 组成，其中基矢量∑反映单元 的刚体平移运动，基矢量 1Λ 反映单元的拉压变形，基矢量 2Λ 和 3Λ 反映单元的剪切变形。 321 ,, ΓΓΓ 和 4Γ 则称为沙漏基矢量,模态如图所示： Γ1 Γ2 Γ3 Γ4 11 当计算单元内任意点的应变率时，公式为： 1 1 11 x x ∂ ∂ = &&ε ， 2 2 22 x x ∂ ∂ = &&ε ， 3 3 33 x x ∂ ∂ = &&ε ΚΚ 。把 2.17式代入即可把速度对总体坐标的偏 导数转化为形函数对总体坐标的偏导数： ( ) 3,2,1,,,8 1 11 = ∂ ∂ = ∂ ∂ ∑ = ix xx x k i k ki && ζηξφ ( ) 3,2,1,,,8 1 22 = ∂ ∂ = ∂ ∂ ∑ = ix xx x k i k ki && ζηξφ ( ) 3,2,1,,,8 1 33 = ∂ ∂ = ∂ ∂ ∑ = ix xx x k i k ki && ζηξφ 由雅可比矩阵转换，又可把速度场对总体坐标的偏导数转化为对自然坐标的偏导数： [ ]       ∂ ∂ ∂ ∂ ∂ ∂ =      ∂ ∂ ∂ ∂ ∂ ∂ − ζ φ η φ ξ φφφφ kkk T kkk J xxx 1 321 这样把求应变率转化为形函数对自然坐标的偏导，如果在单元形心处 )0( === ζηξ 进行单点高斯积分，那么该处速度场对自然坐标的导数（参考公式 2.18）可以表达为： ( ) ( ) ( ){ }txtx kiTTTTi && 00081,,, 43110 ×Γ+×Γ+×Γ+Λ=∂ ∂ === ζηξξ ζηξ ΚΚ 正是由于采用了单元形心处 )0( === ζηξ 进行单点高斯积分，沙漏模态被丢失， 即它对单元应变能的计算没有影响，固称为零能模式，在动力响应计算时，沙漏模态将 不受控制，从而出现计算的数值振荡。 虽然单点积分会引起沙漏模式，但由于只进行单点的积分，相比 222 ×× 或 333 ×× 的高斯积分，单元的数据存储量和计算机时要降低到 1/8 或 1/27，可以大大节 省计算机时，同时也有利于大变形分析。所以依然采用单点积分，不过人为控制沙漏模 式。 LS-DYNA采用沙漏粘性阻尼方式或刚度方式来控制，目前 960版有 7种算法来进 行沙漏阻力的施加。 将各单元节点的沙漏阻力组集成总体结构沙漏阻力，于是就得到了公式 2.4 中的 )( ntH 。通过施加沙漏阻力，沙漏模态在运算中不断得到控制。沙漏的控制方法见第三 章的沙漏控制一节。 csu 矩形