1、汽车仿真分析基础汽车仿真分析基础汽车正面碰撞安全仿真汽车正面碰撞安全仿真汽车正面碰撞安全仿真汽车正面碰撞安全仿真可变性障碍壁汽车偏置碰撞安全仿真可变性障碍壁汽车偏置碰撞安全仿真可变性障碍壁汽车偏置碰撞安全仿真可变性障碍壁汽车偏置碰撞安全仿真汽车侧面碰撞安全仿真汽车侧面碰撞安全仿真安全气囊展开与前挡玻璃作用仿真安全气囊展开与前挡玻璃作用仿真安全气囊展开与驾驶员作用仿真安全气囊展开与驾驶员作用仿真行人头部保护安全仿真行人头部保护安全仿真教学目的和要求教学目的和要求 基本掌握汽车碰撞安全仿真分析软件基本掌握汽车碰撞安全仿真分析软件LS-DYNALS-DYNA的使用方法的使用方法了解隐式、特别是显式有
2、限元方法的相关基础知识了解隐式、特别是显式有限元方法的相关基础知识 理论教学内容理论教学内容 LS-DYNA简介：简介：LS-DYNA软件的发展历程、数据结构、基本使软件的发展历程、数据结构、基本使用方法。用方法。 隐式、显式有限元方法：特别是显式有限元分析的一般方法、时隐式、显式有限元方法：特别是显式有限元分析的一般方法、时间步长的确定、沙漏产生原因及控制方法。间步长的确定、沙漏产生原因及控制方法。 材料模型：包含弹塑性材料、刚体材料在内的碰撞仿真分析常用材料模型：包含弹塑性材料、刚体材料在内的碰撞仿真分析常用材料模型。材料模型。 接触方法：点面接触、面面接触的基本方法和注意事项。接触方法：
3、点面接触、面面接触的基本方法和注意事项。 分析实例：介绍汽车及零部件冲击问题有限元分析实例介绍。分析实例：介绍汽车及零部件冲击问题有限元分析实例介绍。实验教学内容实验教学内容 构造一个钢球撞击一块板的有限元分析模型并进行分析和结果构造一个钢球撞击一块板的有限元分析模型并进行分析和结果处理。了解基于处理。了解基于Hypermesh的有限元网格划分、模型建立、仿真计的有限元网格划分、模型建立、仿真计算和结果处理的基本方法。算和结果处理的基本方法。 平衡方程：不依赖于变形状态的线性方程；平衡方程：不依赖于变形状态的线性方程；几何方程：应变和位移的关系呈线性；几何方程：应变和位移的关系呈线性；物性方程
4、：应力和应变的关系呈线性；物性方程：应力和应变的关系呈线性；边界条件：力边界上的外力和位移边界上的位移是独立的，或边界条件：力边界上的外力和位移边界上的位移是独立的，或线性依赖于变形状态。线性依赖于变形状态。上述方程或边界条件中的任何一个不符合所述特点，就成为非上述方程或边界条件中的任何一个不符合所述特点，就成为非线性问题。线性问题。非线性有限元非线性有限元 非线性有限元非线性有限元 非线性问题：非线性问题：材料非线性材料非线性 塑性域的发生，高温条件下结构的蠕变塑性域的发生，高温条件下结构的蠕变几何非线性几何非线性 大变形（板壳结构的大挠度，屈曲）大变形（板壳结构的大挠度，屈曲）边界非线性边
5、界非线性 接触和碰撞接触和碰撞非线性问题分析步骤：非线性问题分析步骤：1. 明确分析目的，决定模型化方针明确分析目的，决定模型化方针明确分析目的：确认明确分析目的：确认“什么什么” 评价标准（设计规范、安全标准）、评价部位（冲击部位、焊接部位）选择求解器：关键是静态分析还是动态分析？选择求解器：关键是静态分析还是动态分析？ 决定模型形状：形状省略、对称性利用决定模型形状：形状省略、对称性利用选择单元类型选择单元类型材料物性选择材料物性选择确定使用条件确定使用条件选择单位选择单位非线性问题分析步骤：非线性问题分析步骤：2. 构造分析数据（前处理器）构造分析数据（前处理器）网格划分网格划分分析条件
6、设定（载荷、边界、约束和接触等）分析条件设定（载荷、边界、约束和接触等）材料定义（物性和单元断面定义等）材料定义（物性和单元断面定义等）分析控制参数确定（分析结束时间等）分析控制参数确定（分析结束时间等）输出选择设定输出选择设定非线性问题分析步骤：非线性问题分析步骤：3. 执行计算（求解器）执行计算（求解器）分析计算开始分析计算开始计算数据文件确认计算数据文件确认分析计算过程确认分析计算过程确认分析计算结束分析计算结束非线性问题分析步骤：非线性问题分析步骤：4. 分析计算结果研究（后处理器）分析计算结果研究（后处理器）显式有限元分析软件显式有限元分析软件-LS-DYNA John O. Hal
7、lquistPresident and FounderLawrence Software Technology Corporation John Hallquist earned his B.S. in Industrial Engineering, from Western Michigan University in 1970. He received an M.S. degree in Engineering Mechanics from MTU in 1972，and earned a Ph.D. degree in Mechanical Engineering and Enginee
8、ring Mechanics in 1974. John joined the weapons laboratory at Lawrence Livermore National Laboratory (LLNL). 显式有限元分析软件显式有限元分析软件-LSDYNA LS-DYNA的输入数据结构的输入数据结构 KEYWORD数值卡数值卡节点定义节点定义立体单元定义立体单元定义初速度定义初速度定义 面的定义面的定义 定义顺序可任意变动定义顺序可任意变动 *KEYWORD*ENDLS-DYNA的输入数据结构的输入数据结构 材料定义材料定义 单元类型单元类型-体单元体单元 定义体单元关键字需要的参
9、数：定义体单元关键字需要的参数：*ELEMENT_SOLID体单元的体单元的ID号。号。PART的的ID号号8个节点的个节点的ID号号单元公式的指定，通过单元公式的指定，通过*SECTION_SOLIDSECTION的的ID号号体单元计算公式。体单元计算公式。*0专门用于蜂窝材料专门用于蜂窝材料*MAT_MODIFIED_HONEYCOMB的体单元的体单元*1-缺省的常应力单元缺省的常应力单元*2-全积分全积分S/R体单元体单元*ELEMENT_SOLID$# eid pid n1 n2 n3 n4 n5 n6 n7 n8 2843 7 285 284 289 287 6183 6186 61
10、85 6184 2844 7 284 268 267 289 6186 6188 6187 6185 2845 7 287 289 290 288 6184 6185 6190 6189单元类型单元类型-体单元体单元 1-缺省的常应力单元公式，也称为缺省的常应力单元公式，也称为1点积分单元点积分单元积分点在单元的中心，需要控制沙漏。积分点在单元的中心，需要控制沙漏。通常用于通常用于8节点节点6面体单元，但也可用于退化后的面体单元，但也可用于退化后的4节点、节点、6节点单元，但精度有所下降。节点单元，但精度有所下降。对于弯曲板件，厚度方向至少划分两个单元（因为该单元为常应力单元，表示拉压对于弯曲
11、板件，厚度方向至少划分两个单元（因为该单元为常应力单元，表示拉压应力状态至少需要两个单元应力状态至少需要两个单元2-全积分全积分S/R体单元体单元单元内单元内8点高斯积分，没有沙漏发生；但需要点高斯积分，没有沙漏发生；但需要1点积分单元点积分单元4倍的计算时间。倍的计算时间。没有体积锁死（由于采用了选择性缩减积分，且各积分点压力成分相同）没有体积锁死（由于采用了选择性缩减积分，且各积分点压力成分相同）大变形问题时不是很稳定大变形问题时不是很稳定单元类型单元类型-壳单元壳单元 定义壳单元关键字所需要的参数：定义壳单元关键字所需要的参数：*ELEMENT_SHELL壳单元的壳单元的ID号。号。PA
12、RT的的ID号号4个节点的个节点的ID号号单元公式的指定，通过单元公式的指定，通过*SECTION_SHELLSECTION的的ID号号体单元计算公式。体单元计算公式。*1-Hughe-Liu壳单元壳单元*2-Belyschko-Tsay单元单元*16-全积分壳单元（非常快）全积分壳单元（非常快）*ELEMENT_SHELL$# eid pid n1 n2 n3 n4 3639 8 7781 7783 7785 7780 3640 8 7780 7785 7776 7777 3641 8 7769 7786 7784 7768*SECTION_SHELL$HMNAME PROPS 1glass
13、 $# secid elform shrf nip propt qr/irid icomp setyp 1 2$# t1 t2 t3 t4 nloc marea 2.000000 2.000000 2.000000 2.000000单元类型单元类型-壳单元壳单元 节点的逆时针方向确定壳单元正面。节点的逆时针方向确定壳单元正面。板厚方向板厚方向1-5点积分指定可能。点积分指定可能。单元中使用平面应力单元计算公式。单元中使用平面应力单元计算公式。应力应变的输出为积分点位置的值。应力应变的输出为积分点位置的值。单元类型单元类型-壳单元壳单元 1-Hughes-Liu壳单元壳单元LS-DYNA中最早的
14、壳单元，定义来自体单元的退缩形。中最早的壳单元，定义来自体单元的退缩形。面内一点积分，采用面内一点积分，采用Jaumann应力更新，对翘曲的几何体有效。应力更新，对翘曲的几何体有效。能对壳单元的参考面进行偏置（中、内、外）能对壳单元的参考面进行偏置（中、内、外）计算时间是计算时间是BT壳单元的两倍壳单元的两倍单元类型单元类型-壳单元壳单元 2-Belytschko-Tsay壳单元壳单元缺省的壳单元公式，面内单点积分，计算速度很快，对大变形问题是最稳定有效的缺省的壳单元公式，面内单点积分，计算速度很快，对大变形问题是最稳定有效的方法。方法。采用采用Co-rotational应力更新、单元坐标系置
15、于单元中心、基于平面单元假定，所以对应力更新、单元坐标系置于单元中心、基于平面单元假定，所以对于翘曲的几何体不适用。于翘曲的几何体不适用。单元类型单元类型-壳单元壳单元 16-全积分壳单元全积分壳单元以以BT壳单元为基础，面内壳单元为基础，面内4点积分。点积分。只比只比BT壳多壳多2-3倍计算时间。倍计算时间。可以处理翘曲的几何体问题可以处理翘曲的几何体问题,此时需要激活第种沙漏控制公式。此时需要激活第种沙漏控制公式。单元类型单元类型-梁单元梁单元 定义梁单元关键字既需要的参数：定义梁单元关键字既需要的参数：*ELEMENT_BEAM梁单元的梁单元的ID号。号。PART的的ID号号3个节点的个
16、节点的ID号号单元公式的指定，通过单元公式的指定，通过*SECTION_BEAMSECTION的的ID号号梁单元计算公式。梁单元计算公式。1-Hughe-Liu断面积分梁断面积分梁2-Belytschko-Schwer合力梁合力梁3-杆单元杆单元9-可变形焊点梁单元可变形焊点梁单元*ELEMENT_BEAM$# eid pid n1 n2 n3 3639 8 7781 7783 7785 3640 8 7780 7785 7776 3641 8 7769 7786 7784*SECTION_BEAM (Material)$ secid elform shrf qr/irid cst scoor
17、 3 9 0.0000000 0.0000000 1.0000000 0.0000000$ ts1 ts2 tt1 tt2 nsloc ntloc 6.5000000 6.5000000 0.0000000 0.0000000 0.0000000 0.0000000单元类型单元类型-梁单元梁单元 1-Hughe-Liu断面积分梁断面积分梁N1,N2N1,N2两节点确定两节点确定r r轴方向，第三个节点确定断轴方向，第三个节点确定断面内惯性主轴。面内惯性主轴。节点立体单元节点退缩为节点节点立体单元节点退缩为节点单元类型单元类型-梁单元梁单元 1-Hughe-Liu断面积分梁断面积分梁断面内积分点
18、数和断面形状需要定义断面内积分点数和断面形状需要定义断面形状为矩形和圆管可在断面形状为矩形和圆管可在SECTION_SHELLSECTION_SHELL下定义，下定义，工字梁等的其他断面，使用工字梁等的其他断面，使用* *INTEGRATIONINTEGRATION定义定义与轴垂直方向上的应力假定为平面应力状态与轴垂直方向上的应力假定为平面应力状态轴向积分点为单元中心一点，并在此截面计算应力轴向积分点为单元中心一点，并在此截面计算应力单元轴向弯矩一定单元轴向弯矩一定应力输出为局部坐标系。应力输出为局部坐标系。2-Belytschko-Schwer合力梁合力梁不计算梁单元应力，仅计算节点上的力和
19、力矩不计算梁单元应力，仅计算节点上的力和力矩3-杆单元杆单元不需要第三个节点；断面内一点积分不需要第三个节点；断面内一点积分单元类型单元类型-梁单元梁单元 单元类型单元类型-梁单元梁单元 9-可变形点焊梁单元可变形点焊梁单元需要需要*MAT_SPOTWELD材料定义，包括塑性和失效功能材料定义，包括塑性和失效功能单元类型单元类型-梁单元梁单元 9-可变形点焊梁单元可变形点焊梁单元单元类型单元类型-梁单元梁单元 9-可变形点焊梁单元可变形点焊梁单元单元类型单元类型-梁单元梁单元 9-可变形点焊梁单元可变形点焊梁单元单元类型单元类型-弹簧阻尼单元弹簧阻尼单元 *使用两个节点和离散的材料来定义无质量
20、的弹簧阻尼单元使用两个节点和离散的材料来定义无质量的弹簧阻尼单元*能与其他所有显式单元连接能与其他所有显式单元连接*具有平动和转动自由度具有平动和转动自由度*能定义复杂的力能定义复杂的力-位移关系位移关系*由于只能同时定义一个弹簧或阻尼，所以定义弹簧由于只能同时定义一个弹簧或阻尼，所以定义弹簧-阻尼集合时需要重叠定义阻尼集合时需要重叠定义两个单元两个单元单元类型单元类型-弹簧阻尼单元弹簧阻尼单元 求解方法求解方法 l 静态线性问题静态线性问题l 特征值问题特征值问题l 瞬态分析方法：隐式和显式瞬态分析方法：隐式和显式静态线性问题静态线性问题 Solve KD=F for D Gauss eli
21、mination LU decomposition Etc.自由振动分析（特征值分析）自由振动分析（特征值分析） 0DMKD (Homogeneous equation, F = 0)Assume)exp(ti D02 MKLet20 MK0detMKMK K i M i = 0(Eigenvector：特征向量特征向量)(Roots of equation are the eigenvalues:特征值特征值)特征值方程特征值方程自由振动分析（特征值分析）自由振动分析（特征值分析） Methods of solving eigenvalue equation Jacobis method G
22、ivens method and Householders method The bisection method (Sturm sequences) Inverse iteration QR method Subspace iteration Lanczos method瞬态分析方法瞬态分析方法 l Structure systems are very often subjected to transient excitation. l A transient excitation is a highly dynamic time dependent force exerted on the
23、 structure, such as earthquake, impact, and shocks. l The discrete governing equation system usually requires a different solver from that of eigenvalue analysis. l The widely used method is the so-called direct integration method.瞬态分析方法瞬态分析方法 l The direct integration method is basically using the f
24、inite difference method for time stepping.l There are mainly two types of direct integration method; one is implicit and the other is explicit.l Implicit method (e.g. Newmarks method) is more efficient for relatively slow phenomenal Explicit method (e.g. central differencing method) is more efficien
25、t for very fast phenomena, such as impact and explosion.Assume that212ttttttttt DDDDD 1ttttttt DDDDKDCDMDFSubstitute into 212 1tttttttttttttttt K DDDDC DDDMDFNewmarks method (Implicit) 、：常数常数Newmarks method (Implicit)residualcmttttKDFwhere2cmtt KKCM2residual112tttttttttttt FFK DDDC DDTherefore,1cmre
26、sidualttttDKF初始条件：初始条件：000DDD、显式时间积分显式时间积分 ( )mucukup t( )( )mucuk u up t考虑简单的单自由度线形弹簧阻尼系统，根据迖朗贝尔动力学原理：考虑简单的单自由度线形弹簧阻尼系统，根据迖朗贝尔动力学原理：( )MUCUKUP t可以用解析方法求解该常微分方程。若为非线性问题，如可以用解析方法求解该常微分方程。若为非线性问题，如k为为u的函数的函数此时用解析方法一般很难求解，只好应用数值方法来求解，常用的为此时用解析方法一般很难求解，只好应用数值方法来求解，常用的为有限差分法有限差分法和有限元法。和有限元法。上面的公式具有普遍意义，对
27、于有限元法而言，上述方程的矩阵形上面的公式具有普遍意义，对于有限元法而言，上述方程的矩阵形式为：式为：求解该方程使用两种方法：求解该方程使用两种方法：振型叠加法；逐步积分法振型叠加法；逐步积分法。逐步积分法又可分为增量法、迭代法和混合法。逐步积分法又可分为增量法、迭代法和混合法。隐式求解方法一般采用增量迭代法，需要转置换刚度矩阵，通过一系列线性逼近求隐式求解方法一般采用增量迭代法，需要转置换刚度矩阵，通过一系列线性逼近求解，对于存在内部接触这样的高度非线性动力学问题往往无法保证收敛。解，对于存在内部接触这样的高度非线性动力学问题往往无法保证收敛。有限中心差分方程有限中心差分方程显式时间积分显式
28、时间积分 ( )mucukup t( )( )mucuk u up t考虑简单的单自由度线形弹簧阻尼系统，根据迖朗贝尔动力学原理：考虑简单的单自由度线形弹簧阻尼系统，根据迖朗贝尔动力学原理：( )MUCUKUP t可以用解析方法求解该常微分方程。若为非线性问题，如可以用解析方法求解该常微分方程。若为非线性问题，如k为为u的函数的函数此时用解析方法一般很难求解，只好应用数值方法来求解，常用的为此时用解析方法一般很难求解，只好应用数值方法来求解，常用的为有限差分法有限差分法和有限元法。和有限元法。上面的公式具有普遍意义，对于有限元法而言，上述方程的矩阵形上面的公式具有普遍意义，对于有限元法而言，上
29、述方程的矩阵形式为：式为：求解该方程使用两种方法：求解该方程使用两种方法：振型叠加法；逐步积分法振型叠加法；逐步积分法。逐步积分法又可分为增量法、迭代法和混合法。逐步积分法又可分为增量法、迭代法和混合法。隐式求解方法一般采用增量迭代法，需要转置换刚度矩阵，通过一系列线性逼近求隐式求解方法一般采用增量迭代法，需要转置换刚度矩阵，通过一系列线性逼近求解，对于存在内部接触这样的高度非线性动力学问题往往无法保证收敛。解，对于存在内部接触这样的高度非线性动力学问题往往无法保证收敛。位移位移速度速度时间时间12显式时间积分显式时间积分 ( )( )( )( )( )nnnnnMU tP tF tH tCU
30、 t( )nP t( )nCU tLS-DYNA采用显式中心差分法来进行时间积分。在已知采用显式中心差分法来进行时间积分。在已知0,tn时间步阶时间步阶的情况下，求解的情况下，求解tn+1时间步的解，运动方程为：时间步的解，运动方程为：( )nH t( )nF t：外力矢量：外力矢量：内力矢量：内力矢量：HOURGLASS（沙漏）阻力矢量（沙漏）阻力矢量：粘性阻力矢量：粘性阻力矢量位移位移速度速度时间时间1234567显式时间积分显式时间积分 ctt 4/3solidKGcmin/ctlc显式中心差分法的稳定是有条件的（库仑条件）：显式中心差分法的稳定是有条件的（库仑条件）：使用中心差分法进行
31、时间积分时的时间增量，必须比在介质中传播的应力使用中心差分法进行时间积分时的时间增量，必须比在介质中传播的应力波穿过任何单元所需要的时间短一些。波穿过任何单元所需要的时间短一些。不满足库仑条件，就不能保证解的数值稳定性。设构成模型的各单元应力波通过时不满足库仑条件，就不能保证解的数值稳定性。设构成模型的各单元应力波通过时间的最小值为间的最小值为ct应力波的传播速度计算公式依单元形式而异：应力波的传播速度计算公式依单元形式而异：beamEc2(1)shellEc显式时间积分显式时间积分 在运动方程中有一项沙漏阻力：在运动方程中有一项沙漏阻力：该项阻力是人为加上的力，为了防止沙漏变形。沙漏变形是由
32、于该项阻力是人为加上的力，为了防止沙漏变形。沙漏变形是由于LS-DYNA采用单点采用单点高斯积分进行时间积分而导致的高斯积分进行时间积分而导致的零能变形模式零能变形模式。如果使用一点积分，在计算应变成份时即使节点发生位移应变值为如果使用一点积分，在计算应变成份时即使节点发生位移应变值为0，因而不产生，因而不产生应力的现象。这种现象如果任其泛滥，除变形失真外，单元代表长度变短，时间步应力的现象。这种现象如果任其泛滥，除变形失真外，单元代表长度变短，时间步长变小；单元形状发生畸变，计算精度降低。长变小；单元形状发生畸变，计算精度降低。( )( )( )( )nnnnnMU tP tF tCH tU
33、 t显式时间积分显式时间积分 以立体单元为研究对象：以立体单元为研究对象： 81itiuN uNu 88111, ,21212jiijiixyzjiijjittijjiuuxx y zuu u uxxNNuuxxBuBu 123,.kiiiiiNBBBBxiu：节点：节点 在在i方向上的位移方向上的位移显式时间积分显式时间积分 用中心一点积分进行应变评价时，下式成立：用中心一点积分进行应变评价时，下式成立：17352846,iiiiiiiiBBBBBBBB 由于对角线上的节点位移相等，应变为零由于对角线上的节点位移相等，应变为零显式时间积分显式时间积分 显式时间积分显式时间积分 显式时间积分显
34、式时间积分 接触问题接触问题 分析对象的一部分与另一部分接触而发生相互作用时，需要在解除部分定分析对象的一部分与另一部分接触而发生相互作用时，需要在解除部分定义接触面，不然的话，原本相接触的各物体就会毫无阻力地贯穿，而不发义接触面，不然的话，原本相接触的各物体就会毫无阻力地贯穿，而不发生相互作用。生相互作用。接触面的定义：一面为主面（接触面的定义：一面为主面（Master），），另一面为从面（另一面为从面（Slave）。）。每一每一接触面分别由三节点或四节点构成的接触面分别由三节点或四节点构成的segment组成，主面称为组成，主面称为Master segment，从面称为从面称为slave
35、segment。构成构成segment的节点分别称为的节点分别称为slave node和和master node。接触判断是以节点为单位进行。接触判断是以节点为单位进行。同一物体的不同部分相接触的接触面，成为单面接触，只要定义从面。同一物体的不同部分相接触的接触面，成为单面接触，只要定义从面。在被定义的接触面上，在被判断为接触后的节点和在被定义的接触面上，在被判断为接触后的节点和segment上，将作用由上，将作用由单元刚度计算的接触反力。同时还可以考虑摩擦力。单元刚度计算的接触反力。同时还可以考虑摩擦力。接触问题接触问题 在在LS-DYNA中，有四十多种接触处理方法。根据分析模型的不同，有多
36、中，有四十多种接触处理方法。根据分析模型的不同，有多种处理主从面关系的方法。但归纳起来，有以下种处理主从面关系的方法。但归纳起来，有以下4种。种。滑动接触（滑动接触（sliding only）：）：接触面不分离，无摩擦滑动。接触面不分离，无摩擦滑动。固结接触（固结接触（tied interface）：）：接触面被刚接，没有相对位置变化。接触面被刚接，没有相对位置变化。一般接触（一般接触（slide & void）：）：接触面时而接触，时而分离。接触面时而接触，时而分离。单面接触（单面接触（single surface）：）：一个面不同部分的相互接触。一个面不同部分的相互接触。还有将接触面上发生
37、严重畸变的单元自动剔除的接触类型。还有将接触面上发生严重畸变的单元自动剔除的接触类型。接触问题接触问题 接触面的处理方法，大致可分为如下三类。接触面的处理方法，大致可分为如下三类。Kinematic Constraint Method：节点约束法，用于节点约束法，用于tied interface。Penalty Method：惩罚法。惩罚法。Distributed Parameter Method：分布参数法，用于滑动接触类型。分布参数法，用于滑动接触类型。接触问题接触问题 接触计算过程中花费时间最多的是接触点的搜索（接触计算过程中花费时间最多的是接触点的搜索（slave search）。）。
38、从输从输入的大量入的大量segment和节点中，每一循环都必须调查相接触的节点和和节点中，每一循环都必须调查相接触的节点和segment。所以计算效率的提高非常重要。所以计算效率的提高非常重要。Ls-dyna中，根据接触类型的不同，有如中，根据接触类型的不同，有如下两种节点搜索方法。下两种节点搜索方法。利用利用segment形状连续性的方法（形状连续性的方法（connectivity）。）。利用调查空间节点位置的利用调查空间节点位置的bucket sort方法。方法。利用利用bucket sort的方法，的方法，segment的方向性不需考虑时，计算效率高。的方向性不需考虑时，计算效率高。接触
39、类型分类接触类型分类 接触计算过程包含以下三步：接触计算过程包含以下三步：接触计算用接触计算用segment数据定义数据定义接触面搜索接触面搜索被定义构成接触面的节点，或被定义构成接触面的节点，或segment各节点与哪个各节点与哪个segment接触？接触？接触反力计算，消除接触面渗透接触反力计算，消除接触面渗透接触类型分类接触类型分类 接触类型可以按下述方法分类接触类型可以按下述方法分类按接触点搜索方法分类：按接触点搜索方法分类：connectivity, bucket sort接触反力计算方法：接触反力计算方法：penalty, constraint接触面的定义方法：面接触面的定义方法：
40、面-面接触面接触(symmetry)，点点-面接触面接触(non-symmetry)将接触点搜索方法，反力计算和接触面定义组合起来，形成接触类型：将接触点搜索方法，反力计算和接触面定义组合起来，形成接触类型：SURFACE_TO_SURFACE（C, P,S）NODE_TO_SURFACE（C, P,N）AUTOMATIC_SURFACE_TO_SURFACE（B, P, S）CONSTRAINT_SURFACE_TO_SURFACE（B,C, S）AUTOMATIC_SINGLE_SURFACE（B, P, N）使用使用Soft constraint，penalty接触类型使用接触类型使用c
41、onstraint方法计算反力。方法计算反力。1.1.接触面搜索接触面搜索 接触处理是计算中最耗时的工作，而其中最主要的部分是接触面搜索。接触处理是计算中最耗时的工作，而其中最主要的部分是接触面搜索。 有接触可能的零件，有接触可能的零件，segmentsegment在输入文件中已经定义，接触面搜索的工作在于确定哪在输入文件中已经定义，接触面搜索的工作在于确定哪个部分与哪个部分相接触。在程序内部，调查某一特定的个部分与哪个部分相接触。在程序内部，调查某一特定的slaveslave节点与哪个节点与哪个master master segmentsegment接触。接触。 为了选出相接触的为了选出相接
42、触的segmentsegment，需要找到与，需要找到与slaveslave节点最近的节点最近的mastermaster节点（构成节点（构成master master segmentsegment的节点）。这个工作是单纯地计算的节点）。这个工作是单纯地计算slaveslave节点和各节点和各mastermaster节点间的距离，选出最节点间的距离，选出最近的节点，但这需要花费许多时间，因此，需设法节约计算时间。近的节点，但这需要花费许多时间，因此，需设法节约计算时间。 接触面搜索方法特征接触面搜索方法特征1.1.接触面搜索接触面搜索 方法之一是利用单元形状的连续性。将上一个时间步接触的方法之一
43、是利用单元形状的连续性。将上一个时间步接触的master segmentmaster segment事先记事先记忆下来，首先确认该忆下来，首先确认该segmentsegment是否接触。如果没有可能接触，就从与是否接触。如果没有可能接触，就从与slaveslave节点最靠近的、节点最靠近的、构成构成master segmentmaster segment的节点开始顺序搜索。的节点开始顺序搜索。 方法之二是利用方法之二是利用bucket sortbucket sort方法。这个方法是将节点分成组，在每个组内计算节点方法。这个方法是将节点分成组，在每个组内计算节点距离，达到减少计算工作量的目的。距
44、离，达到减少计算工作量的目的。 无论使用哪种方法，第一步的工作都是要找到最近的无论使用哪种方法，第一步的工作都是要找到最近的mastermaster节点。节点。 接触面搜索方法特征接触面搜索方法特征2.2.利用单元形状的连续性利用单元形状的连续性 利用利用master segmentmaster segment形状的连续性计算接触点。在接触范围限定的情况下，效率高，形状的连续性计算接触点。在接触范围限定的情况下，效率高，但当网格不连续时，可能带来接触不安定。但当网格不连续时，可能带来接触不安定。接触面搜索方法特征接触面搜索方法特征3.3.利用三维节点的空间布置的方法利用三维节点的空间布置的方法
45、 特征：接触范围广时，搜索高效，不需考虑接触面的方向。特征：接触范围广时，搜索高效，不需考虑接触面的方向。 说明：将包含网格的空间分为许多立方体（说明：将包含网格的空间分为许多立方体（MAX5000MAX5000）；搜索包含）；搜索包含mastermaster，slaveslave的的立方体。立方体。 分析错误：分析错误：立方体最大为立方体最大为50005000，立方体大小依存于网格。如果有单元发散而产生大，立方体大小依存于网格。如果有单元发散而产生大变形时，一个立方体内包含变形时，一个立方体内包含segmentsegment将超过将超过20002000，程序会给出错误信息，程序会给出错误信息
46、（more than more than 2000 nodes in bucket)2000 nodes in bucket)。Bucket sortBucket sort并非每个时间步都做，并非每个时间步都做，single_surfacesingle_surface（2525），），automatic_single_surface(100)automatic_single_surface(100)。这种搜索方法与。这种搜索方法与segmentsegment指向无关，壳单元必须考虑指向无关，壳单元必须考虑板厚。板厚。接触面搜索方法特征接触面搜索方法特征3.3.接触面搜索方法不同的影响接触面搜索
47、方法不同的影响 接触面搜索方法特征接触面搜索方法特征3.3.接触面搜索方法不同的影响接触面搜索方法不同的影响 接触面搜索方法特征接触面搜索方法特征由接触面搜索，确定与由接触面搜索，确定与slave节点接触的节点接触的segment后，下面计算后，下面计算slave节点与该节点与该segment的接触位置。的接触位置。接触反力计算方法特征接触反力计算方法特征确定接触位置后，给确定接触位置后，给slave节点和节点和master节点加上接触反力。接触反力的计算方法有节点加上接触反力。接触反力的计算方法有两种。两种。 接触面反力计算方法特征接触面反力计算方法特征(1)penalty方法方法(2) C
48、onstraint (2) Constraint 方法方法2sifklnf KAkV2slfm ndtk:接触弹簧刚度，接触弹簧刚度，l:接触面穿透深度，接触面穿透深度，n:接触面法向矢量，接触面法向矢量，fsi: penalty系数，系数，K:体积弹性模量，体积弹性模量，A:接触面积，接触面积，V:接触单元体积，接触单元体积，dt:时间增量，时间增量，ms:slave节点集中质量节点集中质量(1) Penalty(1) Penalty方法的特征：方法的特征： 该方法的关键在于在接触节点和接触面上设定有垂直弹簧。弹簧刚度以接触面上的该方法的关键在于在接触节点和接触面上设定有垂直弹簧。弹簧刚度以
49、接触面上的单元为单位计算，而且用户可以进行刚性参数的调节。单元为单位计算，而且用户可以进行刚性参数的调节。最初，确认最初，确认slaveslave节点是否穿透主面，如果没穿透，当然不进行接触反力计算。如有节点是否穿透主面，如果没穿透，当然不进行接触反力计算。如有穿透，接触反力被加载在从节点和主面节点之间，接触反力的大小与穿透量成反比。穿透，接触反力被加载在从节点和主面节点之间，接触反力的大小与穿透量成反比。接触反力计算方法特征接触反力计算方法特征( , )iln tr ( , )/iirrrrnn (1) Penalty(1) Penalty方法的特征：方法的特征： 接触反力计算方法特征接触反
50、力计算方法特征siiflk n ( , )imisff If l02siiiiisiiiiif K AkVf K AklIf l0立体单元立体单元壳单元壳单元li:壳单元对角线长，壳单元对角线长，fsi: penalty系数系数（通常取（通常取0.1），），Ki:体积弹性模量，体积弹性模量，Ai:接触接触segment面面积，积，Vi:接触单元体积接触单元体积接触面定义接触面定义面面-面接触面接触(symmetry)，点点-面接触面接触(non-symmetry)SURFACE_TO_SURFACENODE_TO_SURFACEAUTOMATIC_SURFACE_TO_SURFACECONST
