弹性应力波的时空多尺度计算方法

弹性应力波的时空多尺度计算方法

许多复杂的力学问题，如裂缝边缘的应力场问题和热耦合结构的局部破坏，很难使用单一的数值方法获得令人满意的模拟结果。发展多尺度计算,针对不同区域采用不同的数值计算方法能够有效解决这类问题,现已成为当今数值模拟研究的一个热点。在固体结构损伤和破坏的力学问题中,充分考虑离散元法和有限元法的优势和弊端,即离散元法能够方便地呈现局部细观非均匀性,但需要大量计算时间;有限元法善于快速求解整体结构宏观弹性变形,但难以准确处理局部结构细微损伤、断裂等力学现象,将这两种方法耦合实现多尺度计算,分别应用于宏观弹性变形区和局部损伤破坏区,能够准确、高效率地分析这些复杂的力学问题。目前的研究成果主要集中在空间多尺度方面,很少进行时间多尺度以及时间与空间兼顾的多尺度计算方法的研究。事实上,时间多尺度与空间多尺度对于提高多尺度计算效率具有同等重要作用。本文采用力学与动力学参数传递以及权值分配法,建立了离散元与有限元耦合的时空多尺度计算方法,既能全面实现空间多尺度,同时又兼顾时间多尺度,有效地提高了多尺度数值计算的效率。同时,作者将这一时空多尺度计算方法应用于激光辐照下受拉铝板破坏行为的数值模拟研究中,取得了较好的效果。1模型和方法的计算1.1有限元时间步长的限制条件离散元法与有限元法具有相似形式的求解方程,都可以归结为一般的牛顿运动定律形式,课题组曾根据这一理论基础,通过引入过渡层概念,对元与节点重叠的有限元和离散元耦合空间多尺度方法(以下简称空间多尺度法)进行了研究。研究发现,这种多尺度法只部分实现了空间多尺度,同时尚未实现时间多尺度计算,数值模拟仍需要较多计算时间,多尺度高效率优势尚未得到较好发挥。离散元法与有限元法的时间步长直接影响时间多尺度计算的实现。离散元法的时间步长存在几何和物理两方面的限制。几何方面,在一个时间步长内不允许一个元与其邻居元沿它们的中心轴线相互贯穿对方;物理方面,在一个时间步长内所产生的动量转移不应超过整个碰撞过程中全部动量的转移。离散元法计算的临界时间步长应当取两者中的最小值。另外,为了提高碰撞过程的时间分辨率,计算中应该选取大于1的时间步数,这样离散元数值模拟的时间步长限制为Δtd≤Δtdcr=1kmini,j{ΔtGdcr,ΔtPdcr}(1)Δtd≤Δtdcr=1kmini,j{ΔtdcrG,ΔtdcrΡ}(1)式中:Δtd为离散元时间步长;Δtdcr为离散元临界时间步长;k为时间步数,k>1;i和j为两个相邻元;ΔtGdcr为离散元几何限制临界时间步长;ΔtPdcr为离散元物理限制临界时间步长。由于离散元法采用显式差分格式求解,为了耦合两种方法实现多尺度计算,有限元法的计算采用显式求解法中的中心差分法。为了保证有限元求解方程的稳定性,时间步长需要满足条件Δtf≤Δtfcr=LminE(1−ν)ρ(1+ν)(1−2ν)√(2)Δtf≤Δtfcr=LminE(1-ν)ρ(1+ν)(1-2ν)(2)式中:Δtf为有限元时间步长;Δtfcr为有限元临界时间步长;Lmin为模型中网格单元最小尺寸;ρ为材料密度;E为杨氏模量;ν为泊松比。从离散元与有限元时间步长的限制条件可知,离散元时间步长取决于元尺寸、邻居元距离以及元性质,所以离散元时间步长是很小的。有限元时间步长不仅与结构或材料的性质有关,而且与模型的网格单元最小尺寸成正比。因此,为了实现时间多尺度计算,离散元采用固定较小时间步长的同时,有限元需要采用较大的时间步长,这时根据式(2)就需要模型具有较大的有限元网格单元最小尺寸。1.2有限元时间步长求解在空间多尺度法的基础上,基于以上思想作者建立了既能有效实现时间多尺度,同时又能完善空间多尺度的离散元与有限元耦合的时空多尺度计算方法,这种方法的计算模型如图1所示。从图中可以看出,在耦合过渡区中,有限元网格单元的最小尺寸不必受到元的大小及元间距的限制,这样网格单元的最小尺寸可以取得较大,从而增大有限元时间步长,而这时离散元可以采用较小的时间步长,完成时间多尺度计算。另外,该模型中整体有限元网格的数量受耦合过渡区网格尺度的影响较小,具备完善的空间多尺度计算能力。在计算模型的耦合过渡区中,离散元和有限元相互为对方提供边界条件,并采用权值分配法进行力学和动力学参数的传递。也就是说,离散元为有限元提供力边界条件采用式(3)进行计算,有限元为离散元提供的速度边界条件采用式(4)进行计算。当耦合过渡区内施加外载荷时,需要将外载荷加到式(3)和式(4)的网格节点或者元的相应变量中。F(D)FE=∑i=1nk(m)iFi(m)DE+∑p∑l<m,i=1qk(l)iFi(l)DE(3)FFE(D)=∑i=1nki(m)FDEi(m)+∑p∑l<m,i=1qki(l)FDEi(l)(3)式中:F(D)FEFE(D)为离散元提供力边界条件后的网格节点作用力;k(m)ii(m)为第m网格内第i元在该网格节点的权值;Fi(m)DEDEi(m)为第m网格内第i元受到的合力;n表示第m网格包含的元数目;k(l)ii(l)为第l(l<m)网格内第i元在该网格节点的权值;Fi(l)DEDEi(l)为第l(l<m)网格内第i元受到的合力;q为第l网格包含的元数目;p为包含该网格节点的所有l(l<m)网格的数目。V(F)DE=∑i=13kiViFE(4)VDE(F)=∑i=13kiVFEi(4)式中:V(F)DEDE(F)为有限元提供速度边界条件后的元速度;ki为该元所属网格的第i节点(三角形网格,i=1,2,3)在该元位置处的权值;ViFEFEi为该元所属网格单元的第i节点的速度。为了实现耦合过渡区材料密度的匹配性,需要使过渡区内网格节点的质量与相应元过渡到该节点的质量相等,计算方法为M(D)FE=∑i=1nk(m)iMi(m)DE+∑p∑l<m,i=1qk(l)iMi(l)DE(5)ΜFE(D)=∑i=1nki(m)ΜDEi(m)+∑p∑l<m,i=1qki(l)ΜDEi(l)(5)式中:M(D)FEFE(D)为元质量传递后的网格节点质量;Mi(m)DEDEi(m)为第m网格内第i元的质量;Mi(l)DEDEi(l)为第l(l<m)网格内第i元的质量。在时空多尺度数值计算中,离散元时间步长设为dt,有限元时间步长设为fdt,两者关系为fdt/dt=N。也就是说,在每个有限元时间步长求解计算中,执行N次离散元计算过程来完成时间多尺度计算。整个系统的计算流程如图2所示。2有限元数值模拟结果用弹性应力波传播算例来验证时空多尺度计算方法的有效性和优越性。选用一根长度、宽度和厚度分别为3.2×10-1m、2.4×10-3m和3.2×10-4m的细长杆。杆的右端固定,左端施加一个突然加载和卸载的分布载荷,载荷的大小为1.0×104N/m,脉宽为1.5×10-5s。杆的材料为铝,密度为2.7×103kg/m3,弹性模量为7.0×1015Pa。在时空多尺度数值计算中,将杆中部4.0×10-3m的区域用离散元模拟,而其余部分用有限元模拟。在细长杆上取3个不同位置(A、B、C)来观察应力波的传播情况,这3个位置分别距离杆左端7.80×10-2m,1.60×10-1m和2.42×10-1m,如图3所示。图4为使用单一有限元法、时空多尺度法以及解析法时应力波在细长杆A、B、C三个位置的传播情况。从图中可以看出,时空多尺度法与有限元法相比,两者在三个位置获得的应力波传播情况基本一致,应力波能够在有限元区(A和C位置)和离散元区(B位置)之间光滑过渡。同时,时空多尺度法计算的应力波到达三个位置的时间及峰值与解析解基本吻合。以上结果验证了离散元与有限元耦合的时空多尺度计算方法的有效性。由于引入了阻尼,波在传播过程中有一定的衰减。由于杆的横向惯性效应,波形将出现剧烈振荡,图4中的模拟结果很好地说明了这一点。为了研究时空多尺度方法的计算效率,将空间多尺度法与耦合过渡区内布设不同层离散元的时空多尺度法在弹性应力波传播算例中的计算时间进行了对比,如表1所示。根据式(2)可知,有限元临界时间步长与网格单元最小尺寸成正比。为了对比分析方便,在时空多尺度法中,耦合过渡区网格尺寸(不大于非耦合过渡区网格最小尺寸)为离散单元间距的倍数,令耦合过渡区内网格尺寸与元间距相等时,有限元与离散元的时间步长相等,随着元层数增加,有限元时间步长相应增倍。从表1中可以看出,当采用时空多尺度法时,随着耦合过渡区内离散元层数的增多以及有限元时间步长的增大,计算时间明显减小。并且,时空多尺度法所有情况下的计算时间均明显少于空间多尺度法的计算时间,当耦合过渡区内布置五层离散元时,采用时空多尺度法的计算时间仅相当于空间多尺度法计算时间的1/3。可见,时空多尺度法在节省计算时间方面具有较大的优越性,有效地发挥了多尺度数值计算的高效率优势。3受拉铝板微裂纹的有限元模拟研究激光辐照下预载金属结构的破坏机理对于工业生产和军事应用具有重要的科学意义。本文使用离散元与有限元耦合的时空多尺度计算方法对激光辐照下受拉铝板破坏行为进行数值模拟。选用长度、宽度和厚度分别为1.0×10-1、5.92×10-2和6.5×10-4m的矩形薄板,板的材料为铝合金,密度为2.7×103kg/m3,弹性模量为7.0×1015Pa,熔点为993K,热导率为233W/(m·K),热膨胀系数为2.9×10-7K-1。薄板沿y轴受到双向均匀预拉载荷,载荷力的大小为该板静态破坏载荷Fb(大小为1.4×105N/m)的1/4,同时薄板中部受到高斯分布的激光辐照,光斑半径为1.5×10-2m,激光束平均功率密度为2.5×106W/m2,辐照表面吸收系数为0.15。由于薄板中部受到热力耦合的作用,将发生损伤、断裂等非均匀破坏行为,采用离散元进行模拟,共有5289个离散元;而远离激光辐照的区域,其应力状态仍在弹性范围,采用有限元进行模拟,共有294个有限元网格,计算模型如图5所示。图6为激光辐照下受拉铝板破坏过程的时空多尺度数值模拟图,图中的黑色三角形标出微裂纹的位置。从图中可以看出,微裂纹首先出现在垂直于拉伸方向的光斑直径与光斑边缘的两个交点附近,随后以该微裂纹处作为裂源,在其周围萌生大量新微裂纹,并沿横向迅速向铝板边缘和光斑中心双向扩展。微裂纹的扩大将形成宏观裂纹,并在很短的时间内沿水平方向贯穿整个铝板,直至完全断裂,这些模拟结果与文献提供的实验结果基本一致。4多尺度计算方法数值模拟本文建立了离散元与有限元耦合的时空多尺度计算方法,通过弹性应力波传播算