abaqus显示和隐式算法的差别

1.1.弹性动力学有限元基本解法结构系统的通用运动学方程为：tRKUUCUM（1）求解该动力学振动响应主要有三类方法：（1）时域法（2）频域法（3）响应谱法时域法又可分为：（1）直接积分法，（2）模态叠加法。直接积分法又可分为中心差分法（显式），Wilson（隐式）法以及Newmark（隐式）法等。本文介绍中心差分法（显式）与Newmark（隐式）法。1中心差分法（显式）假定0，1t，2t，…，nt时刻的节点位移，速度与加速度均为已知，现求解)(tttn时刻的结构响应。中心差分法对加速度，速度的导数采用中心差分代替，即为：)2(12ttttttUUUtU)(21tttttUUtU(2)将（2）式代入（1）式后整理得到tttRUMˆˆ（3）式（3）中CtMtM211ˆ2tttttUCtMtUMtKRR)211()2(ˆ22分别称为有效质量矩阵，有效载荷矢量。R，M，C，K为结构载荷，质量，阻尼，刚度矩阵。求解线性方程组（3），即可获得tt时刻的节点位移向量ttU，将ttU代回几何方程与物理方程，可得tt时刻的单元应力和应变。中心差分法在求解tt瞬时的位移ttU时，只需tt时刻以前的状态变量tU和ttU，然后计算出有效质量矩阵Mˆ，有效载荷矢量tRˆ，即可求出ttU，故称此解法为显式算法。中心差分法，在开始计算时，需要仔细处理。t=0时，要计算tU，需要知道tU的值。因此应该有一个起始技术，因而该算法不是自动起步的。由于0U，0U，0U是已知的，由t=0时的（2）式可知：02002UtUtUUt中心差分法中时间步长t的选择涉及两个方面的约束：数值算法的稳定性和计算时间。中心差分法的实质是用差分代替微分，并且对位移和加速度的导数采用线性外插，这限制了t的取值不可过大，否则结果可能失真过大。可以证明：中心差分法是条件稳定的。即当时间步长t必须小于由该问题求解方程性质所决定的一个时间步长的临界值。LS-DYNA中，采用“变时间步长法”，即每一时刻的步长t由当前结构的稳定性条件来控制。具体算法为：计算每一个单元的极限时间步长eit，i=1,2…，取)min(eitt为下一个时刻的时间步长。各种单元的eit计算方法如下。（1）1D杆，梁单元CLte其中为时间步长因子，系统默认为0.9。L为杆，梁单元的长度。EC材料声速。（2）2D板，壳单元CLtemin其中为时间步长因子，系统默认为0.9。minL为壳单元的最小单元边长度。)1(2EC为材料的声速。（3）3D单元2/122)(CQQLtee其中)0(0)0(01kkkkkkeLCCCQ0C和1C为无量纲常数，默认0C=1.5，1C=0.06.节点体单元）（对节点体单元）对48(minmaxLAVLeeeeL为单元等效长度，eV为单元体积，maxeA为单元最大侧面积。）（）（2-1)1(-1EC为材料声速。时间步因子可由用户设置，减小相当于减少时间步长。设置时间步长因子的关键字为*CONTROL_TIMESTEP，控制参数为TSSFAC。另外，质量缩放可以人为控制时间步长。即调整单元密度来改变时间步长。以壳单元为例说明质量缩放改变时间步长。CLtemin（假如=1）)1(2EC由ELtispecified)1()(22min得到)1()(22min2LEtspecifiediLS-DYNA中有2种质量缩放方案，修改*CONTROL_TIMESTEP中的参数。方案1：DT2MS为正的时间步通过调整单元密度，使所有单元时间步相同，只用于惯性效应不重要的情况。方案2：DT2MS为负的时间步质量缩放只用于小于指定时间步长DT的单元。惯性效应应该通过变形体的动能与内能的比例进行衡量（一般应该小于10%等）。2Newmark法（隐式）Newmark假定在时间间隔],[ttt内，加速度线性变化，即采用如下的加速度，速度公式：tUUUUtttttt])1[(2])21[(tUUtUUUttttttt(4)式中，为按积分的精度和稳定性要求可以调整的参数。根据（4）式可给出ttU和ttU用ttU，tU，tU表示的表达式，代入（1）式中整理得到ttttRUKˆˆ（5）其中KCtMtK21ˆ])12()1([])121(11[ˆ2ttttttttttUtUUtCUUtUtMRR称之为有效刚度矩阵和有效载荷矢量。由上式可以看出求解当前ttU，需要用到当前时刻的ttR，因此该算法为隐式算法。当载荷历史全部已知时，ttF为已知量，求解需要迭代实现。可以证明，当参数5.0，2)5.0(25.0时，Newmark法是无条件稳定的，即t的大小不影响数值稳定性。此时时间步长t的选择主要根据解得精度确定。一般，Newmark法可以比中心差分法的时间步长大得多。3．结论比较两种算法，显式中心差分法非常适合研究波的传播问题，如碰撞、高速冲击、爆炸等。分析式（3）发现，显式中心差分法的M与C矩阵是对角阵，如给定某些有限元节点以初始扰动，在经过一个时间步长后，和它相关的节点进入运动，即U中这些节点对应的分量成为非零量，此特点正好和波的传播特点相一致。另一方面，研究波传播的过程需要微小的时间步长，这也正是中心差分法的特点。而Newmark法更加适合于计算低频占主导的动力问题，从计算精度考虑，允许采用较大的时间步长以节省计算时间，同时较大的时间步长还可以过滤掉高阶不精确特征值对系统响应的影响。隐式方法要转置刚度矩阵，增量迭代，通过一系列线性逼近（Newton-Raphson）来求解。正因为隐式算法要对刚度矩阵求逆，所以计算时要求整体刚度矩阵不能奇异，对于一些接触高度非线性问题，有时无法保证收敛。1.2.１显式算法显式算法基本假定为：在一微小时间段内，模型任意点速度、加速度为常数。ＡＢＡＱＵＳ软件Ｅｘｐｌｉｃｉｔ模块应用中心差分法对运动方程进行显式时间积分，运动方程的解为¨ｕ（ｉ）＝Ｍ－１·（Ｆ（ｉ）－Ｉ（ｉ））（１）式中：Ｍ为集中质量矩阵；Ｆ为外荷载向量；Ｉ为单元内力向量。由于显式算法中不需要对刚度矩阵求逆，集中质量矩阵为对角矩阵，求逆简便，使显式算法并行计算数据传输量较小；且显式算法刚度矩阵大小与自由度数成线性关系，因此显式算法用于自由度数庞大的数值计算时具有很大优势。1.3.2隐式算法隐式算法含义为：ｔ＋Δｔ时刻状态不仅与ｔ时刻状态有关，且与ｔ＋Δｔ时刻某些量有关。ＡＢＡＱＵＳ软件Ｓｔａｎｄａｒｄ应用Ｈｉｌｂｅｒ－Ｈｕｇｈｅｓ－Ｔａｙｌｏｒ隐式算法、Ｎｅｗｔｏｎ－Ｒａｐｈｓｏｎ迭代法进行动力方程求解［１０］。运动方程解为Δｕ（ｉ＋１）＝Δｕ（ｉ）＋Ｋt－１·（Ｆ（ｉ）－Ｉ（ｉ））（２）式中：Ｋｔ为当前切线刚度矩阵：Δｕ为位移增量。求解方程位移增量Δｕ（ｉ＋１）时，必须对刚度矩阵Ｋ求逆。当自由度数非常庞大时，这项计算消耗资源。对Ｋ矩阵的求逆计算，计算机之间数据传输量非常大，随着自由度数增加，刚度矩阵Ｋ大小成指数增长。因此，隐式算法用于自由度数庞大的数值计算时，优势不明显，甚至会降低计算效率。根据地铁地下结构抗震研究需要，对有限元并行计算显式算法和隐式算法计算精度和效率进行比较，主要结论如下：（１）有限元并行计算中心差分显式算法与Ｈｉｌｂｅｒ－Ｈｕｇｈｅｓ－Ｔａｙｌｏｒ隐式算法计算精度相当，显式算法计算效率远高于隐式算法。（２）黏弹性人工边界在显式算法和隐式算法中，都能起到很好的模拟效果。（３）显式算法和隐式算法计算相对位移、相对速度时程基本一致，两种算法峰值应力差在９％以内。（４）地铁地下结构抗震分析中，对于自由度很大的三维结构，并行计算显式算法计算效率较高；而对于自由度较小的二维结构，并行计算隐式算法计算效率较高。（５）显式算法适合多处理器并行运算，对庞大自由度的地下结构三维非线性抗震分析，具有较好的并行计算效率。对于ＣＰＵ具体使用数和计算模型自由度的关系，尚需进一步研究。