《离散元与壳体有限元结合的多尺度方法及其应用》由会员分享，可在线阅读，更多相关《离散元与壳体有限元结合的多尺度方法及其应用（8页珍藏版）》请在金锄头文库上搜索。
1、离散元与壳体有限元结合的多尺度方法及其应用去唐志平，胥建龙( 中国科学技术大学力学和机械工程系，中科院材料力学行为和设计重点实验室合肥)麓要= 在深入研究复杂结构和菲均质材料冲击响应和破坏机理的过程中，往往遇到多尺度计算问题。本文尝试建立三维离散元与壳体有限元结合的多尺度方法用于处理圆柱壳问题，该方法采用三维离散元对感兴趣的局域进行局部模拟，利用平板壳体有限元进行宏观的模拟，采用一种特殊的过渡层使离散元区和有限元区能很好的衔接。我们将这一方法应用于激光辐照下充压柱壳的热力耦合冲击破坏响应，得到的模拟结果与文献报道有较好的吻合。关键词：离散元壳体有限元多尺度计算方法激光辐照1 引言对很多现实问题
2、的深入研究，例如结构和材料的损伤破坏问题，往往需要从多个尺度来考察，既要从整体上把握问题的宏观特性，又能够对感兴趣的局部进行细观的模拟，揭示局域破坏过程和机理。多尺度计算是当今数值计算积极研究和探索的主攻方向之一。近三十年发展起来的离散元方法，在模拟岩石边坡运动，颗粒介质力学行为，结构失效，热一力耦合以及化学反应，相变，结构冲击响应等问题中得到了广泛的应用 1 - 6 1 。有限元方法【7 J 作为一种比较成熟的数值模拟方法，在宏观尺度上能够很好的模拟连续介质和结构响应问题。虽然单纯用细观离散元方法也可以模拟某些简单的结构，但其所需要的内存和机时是惊人的，而有限元方法又无法像离散元那样方便的模
3、拟整个的破坏过程。为此我们曾提出离散元和2 维平面有限元结合的方法I s 】，并应用于激光辐照下预应力铝板的破坏响应，模拟结果与实验有较好的吻合。本文尝试将三维离散元与壳体有限元结合起来处理圆柱壳问题。三维离散元在处理问题时是十分灵活的，我们曾成功的利用三维离散元对激光辐照下柱壳的失效问题进行了模拟1 6 1 ，得到了比较满意的结果。但由于受元的数目的限制，原模拟只能取壳体的一部分并采用固定边界，这不可避免的带来边界效应。为了迸一步模拟柱壳动态响应问题，离散元与壳体有限元的结合也许是较好的方法之一。本文首先简要介绍离散元和有限元结合思想和过渡层概念，然后介绍和壳体有限元结合的方法，最后给出激光
4、辐照下充压柱壳失效的算例。2 离散元和有限元结合方法和过渡层离散元方法的基本思想是把介质离散为“元”，相邻的元之间存在某种或某几种作用力，元的运动受牛顿运动定律支2 1 3 , 5 1聊J l e 国家自然科学基金资助项目( 1 0 4 7 2 1 1 4 )4 9的乙净一K+一E秘“艏，庐矿Mv厶一Mv厶一=| l式中聊、J j 分别为元i 的质量和转动惯量，芦、万分别为元i 中心的位置和角矢量，户”是元f 的邻居面对它的作用力，牙”是该作用力到元f 中心的力臂，和分别为f 元所受的外力和外力矩，j 是元f 的邻居元数目，是系统中元的总数。根据初、边值条件求解上述2 N 个矢量常微分方程就得
5、到系统的演化规律。有限元方法是将连续介质划分为有限个单元，从弹性动力学的基本方程出发，可以得到系统的求解方程【7 】施= 足+ 吒( 2 )其中吒= 一C a K a ，石，a ，口分别为系统的节点加速度向量，节点速度向量和节点位移向量，M ，C ，K 分别为系统的质量矩阵，阻尼矩阵和刚度矩阵，艺为节点载荷向量。比较( 1 ) 式和( 2 ) 式，我们可以看到两种方法的求解方程是很相似的，这种相似性为两种方法的结合提供了可能。为了使两种方法有效地结合在同一个计算模型中，关键的问题是要使得离散元区域和有限元区域之间的各种力学、物理特性有一个平滑的过渡。我们引人过渡层的概念来实现两者之间的平滑过渡
6、，如图l 所示【8 】。在过渡层中同时用有限元和离散元将介质离散，过渡层的厚度取为两层离散元的厚度，可保证第一层中的离散元有完整的邻居元，与单纯用离散元模拟的情况相符。过渡层中的离散元与非过渡层中的离散元通过所定义的势函数相互作用。同时，过渡层还由有限元网格所构成，其节点与相应的离散元位置始终保持一致，离散元的位置随有限元节点位置的变化而变化。过渡层中离散元所受邻居元作用的合力，作为节点外力施加在相应的有限元节点上的。除了过渡层和力边界外，其余部分的节点外力均为零。为了消除过渡层材料密度的不匹配，让过渡层中有限元节点的质量与离散元的质量相等。在过渡层之外的离散元区域和有限元区域则分别按照各自已
7、有的方法求解。 - 离散兀区域 5 f 气o 誓 f U I 一乃兀一一 一一 趴扒12 | 过渡层图1 离散元和有限元通过过渡层相结合在上述的模型中可以看到，在每一步计算中，过渡层中的离散元和有限元相互为对方提供边界条件，由于过渡层中离散元的位置始终是与有限元节点保持一致的，过渡层上离散元所受的作用力作为节点外力旌加给有限元，保证了位移、应变和应力的连续性。对于求解步长的稳定性，离散元和有限元有不同的要求。为了得到足够的分辨率，离散元所取的临界步长比有限元的小得多，因此计算中采用离散元所确定的时间步长，同时也保证了有限元解的稳定性【8 】。3 平板壳体单元简介在用有限元进行壳体分析时，一般有
8、平板壳元和曲面壳元两种单元。曲面单元能较好的反映壳体的真实几何形状，在单元尺寸相同的情况下，通常可以得到比平板单元更好的结果，但是曲面单元插值函数的构造是一件很困难的事。平板单元的最主要的优点是表达格式简单，只要计算网格适当的加密，也可以得到较好的精度。本文为了和平面有限元工作衔接8 1 ，采用了三角形的平板壳体单元。平板壳体单元可以看成是平面应力单元和平板弯曲单元的组合，因此其刚度矩阵可以由这两种单元的刚度矩阵组合而成。如图2 所示，在局部坐标系x y z 中，三角形壳体单元的节点位移向量可以表示为：a ，= v jw 色吼包】r( 3 )在平面应力状态下，选用线性平面单元，其节点的位移向量
9、为【v f 】7 ，其刚度矩阵为：K ；，1 ) _ 小B W w D 州B y d x d y其中B y ，D 为平面单元的应变矩阵和弹性矩阵，加上标( 聊) ，表示属于薄膜应力状态。ZX图2 三角形平板壳体有限元的节点位移( 4 )本文中，板弯曲单元选用3 节点9 个位移参数的非协调单元，通过数值实验，这两种单元结合而成的壳体元可以得到较好的结果。在平板弯曲状态下，单元内的挠度w 可以写为：w = c N 。N ：N 。， 兰其中a j = 【w f色巳r ，( f = 1 ，2 ，3 ) 为节点向量，N 为插值函数，具体取如下的形式：小卧厶+ 置厶+ 葺厶一厶置一厶鹭6 2 ( 厶葺+
10、1 L , L ：g ) 一6 3 ( 葺厶+ 三厶厶厶)c ：( 厶葺+ 互I 厶厶厶) 一c 3 ( 葺厶+ 圭厶厶厶)( 5 )厶，厶，厶是三角形的面积坐标，6 2 ，6 3 ，c 2 ，巳是单元参数，可以由单元节点的坐标求得。N 2 ，N 3 可以通过下标轮换得到。由插值函数，根据有限元求解的基本步骤可以得到其刚度矩阵为：I q ) _ f f ( B “ D B 尹d x d y( 7 )其中B ，D ( 6 是平板弯曲单元的应变矩阵和弹性矩阵，上标( 6 ) 表示平板弯曲状态，B 夕，D Q 的具体表达式及推导过程可以参看相关的文献p ，9 1 们。这样，平板壳体单元的刚度矩阵可以
11、表示为：5 lK ，=K 扩，；。0 兰躺00 ；000 ：K 擘：000 ：000 000 0( 8 )得到单元刚度矩阵后，由于局部坐标与总体坐标不重合，还需要经过坐标变换才能集成总体刚度矩阵。设总体坐标为x Y Z ，则单元刚度矩阵的转换关系为：K 。= T K 。T T ，K 。= T T K 。T( 9 ) 其中T =L000L000L陬，L = l。0O l ，九= 九l丑，五。五，：乃，乃t y以t ：丸，也y乃。式中的丑，。= C O S ( X ，x ) 等是局部坐标轴在总体坐标系中的方向余弦。关于三角形壳体元局部坐标系的选取以及方向余弦的求解可以参看文献【7 】中相关部分的内
12、容。同样，单元节点的载苟向量也要通过坐标变换才能形成总的载荷向量，其做法与刚度矩阵的情形类似。在对静态圆柱壳在内压下的情况进行计算的结果表明【l 们，上述平板壳体单元在处理内压问题时可以得到较为满意的结果，因此本文将利用这种壳体单元与离散元结合来模拟充压柱壳在激光辐照下的失效行为的模拟。4 离散元与有限元结合在壳体问题中的应用在将壳单元与离散元结合时有几个问题需要注意：l 、在平面问题中节点自由度只有两个，只要保证位移和力平滑过渡即可，但在壳体元单元中由于节点的自由度增加，在过渡层中也要相应地增加角度，角速度，弯矩等量的平滑过渡。2 、在求质量矩阵时，采用集中质量阵，忽略了转动项。3 、阻尼矩
13、阵依然采用R a y l e i g h 阻尼。4 1 模型描述充压柱壳结构在连续激光辐照下的失效问题的物理模型如图3 a 所示，一半径为R 的圆柱形壳体，厚度为t ，圆柱上部受到一半径为r 的圆形激光光斑照射，柱壳内部受到均匀压力P 的作用。图3 b 和c 为计算模型，中间矩形为包含光斑的离散元部分，其余则为平板壳体元网格。为了验证方法的正确性，首先考察没有激光辐照时圆柱壳在内压作用下的薄膜应力，计算的结果如图4 所示。图4 可见A 、B 两点的薄膜应力基本上在理论解附近，存在误差是由于网格划分的比较粗( 1 3 6 个离散元，8 3 5 个有限元网格) 。A 点的值比B 点稍高，这与离散元
14、与有限元宏观应力的计算方法不同有一定的关系。开始由于动态加载，存在较大的振荡，但应力很快就平稳了，并趋向于理论解。这就说明了用离散元和有限元结合的方法处理薄壳问题的可行性。接着，进一步模拟激光辐照下的响应，为了便于和文献 6 1 q h 的结果相比较，选取直径R = 3 6 0 m m ， 壁厚t = S m m ，内压为P = 5 M P a 的3 0 C r M n S i A 钢和铝合金圆柱壳进行模拟。计算中两种材料的主要 参数如表l 所示( 表中，C 知，b 声为熔点温度、室温比热、熔点比热、热导率和热膨胀系数， O o ，m ，n 为离散元模型中的材料常数【6 8 】) 。5 2(a
15、)to)图3 激光辐照下充压柱壳的物理模型( a ) 和生成的计算模型( b )图4 圆柱壳薄膜应力随时间的发展表1 主要材料参数材料P o ( g c m 3 )暇)c p o ( J ( g K ) )C p m ( J ( g K ) )O t o ( G P 曲mnr ( W ( c m K ) )f l ( 1 0 - 6 K - 1 )F e7 91 6 6 30 4 7 31 9 21 0 4l2O 2 4 51 5A I2 79 9 30 9 0 ll ，6 5O 4 5232 3 32 9将光斑附近2 0 e m x 2 0 c m 的区域采用离散元划分，其余壳体采用有限元划
16、分。离散元直径5 m m ，共有1 3 4 5 个离散元，1 3 9 0 个有限元网格。激光仍采用高斯分布的圆形光斑。柱壳的两端取固支边界条件。4 2 计算结果激光参数( 功率密度而，光斑半径r o ) 和计算结果列于表2 ，并示于图5 至9 。表2 中列出了各种情况下的起始破坏时间乙，和对应的中心温度瓦( 破坏阈值) 等计算结果，表中d m 。是最大延伸度，超过这一极限，元之间将产生裂纹。图5 是损伤数D 随辐照时间的演化。图6 和7 所示的是3 0 C r M n S i A5 3钢的情况，图6 中可以看到光斑辐照区在热软化和内压的联合作用下向外鼓出，并在中心产生小破坏，达到一定时间后，在短时间内发生爆裂破坏。从图5 的损伤演化曲线也可以看出，损伤数D 先维持一个很小的值，然后突然升高。当提高而减小9 0 时，
