第39卷 第2期2009年3月
吉林大学学报(工学版)
JournalofJilinUniversity(EngineeringandTechnologyEdition)
Vol.39 No.2 Mar.2009
离散元与有限元耦合的时空多尺度计算方法
张 锐1,唐志平2,郑 航2
(1.吉林大学地面机械仿生技术教育部重点实验室长春130022;2.中国科学技术大学中国科学院材料力学
行为和设计重点实验室合肥230027)
摘 要:采用力学与动力学参数传递以及权值分配法度计算方法。度双重优势,,,将这一多尺度,得到的模拟结果与实验结果基本一致关键词:;时空多尺度;离散元与有限元耦合;数值计算中图分类号:O242.1 文献标识码:A 文章编号:167125497(2009)0220408205
Timeandspacemultiscalenumericalmethodbycoupling
discreteelementmethodandfiniteelementmethod
ZHANGRui1,TANGZhi2ping2,ZHENGHang2
(1.KeyLaboratoryofTerrain2machineBionicsEngineering,MinistryofEducation,JilinUniversity,Changchun
130022,China;2.KeyLaboratoryforMechanicalBehaviorandDesignofMaterials(LMBD)ofCAS,UniversityofScienceandTechnologyofChina,Hefei230027,China)
Abstract:Throughmechanicalanddynamicparameterstransitionandweightvaluedistribution,atimeandspacemultiscalenumericalmethodbycouplingDiscreteElementMethod(DEM)andFiniteElementMethod(FEM)wasproposed.Simulationofanexampleofelasticstress2wavepropagationwascarriedout.
Theresultsshowthatthemultiscalenumericalmethodpossessesthedual
superioritiesofmultiscaletimeandmultiscalespace,whichsignificantlyimprovetheefficiencyofsimulation.Itisshownthatthesimulationresultsbythemethodsatisfytherequirementofaccuracy.Themethodisalsoemployedtosimulatethedamagebehaviorsofatensilealuminumplateunderlaserirradiation,andtheresultsareaccordantwiththeexperimentalresults.
Keywords:numericalmathematics;multiscaletimeandspace;couplingdiscreteelementmethodandfiniteelementmethod;numericalcomputation
许多复杂力学问题,例如裂纹尖端应力场问题、结构在热力耦合载荷作用下局部破坏问题等,
收稿日期:2007208222.
难以用单一尺度数值方法得到满意的模拟结果。发展多尺度计算,针对不同区域采用不同的数值
基金项目:国家自然科学基金项目(50805064,10472114);吉林省科技发展计划项目(20080146);2006年度中国科学
院王宽诚博士后工作奖励基金项目.
作者简介:张锐(19752),男,讲师,博士.研究方向:数值计算.E2mail:[email protected]
第2期张 锐,等:离散元与有限元耦合的时空多尺度计算方法
・409・
(1)
计算方法能够有效解决这类问题,现已成为当今数值模拟研究的一个热点
[1]
。在固体结构损伤和
GP
Δtd≤Δtdcr=min{Δtdcr,Δtdcr}
k
i,j
破坏的力学问题中,充分考虑离散元法和有限元法的优势和弊端,即离散元法能够方便地呈现局部细观非均匀性,但需要大量计算时间;有限元法善于快速求解整体结构宏观弹性变形,但难以准确处理局部结构细微损伤、断裂等力学现象,将这两种方法耦合实现多尺度计算,分别应用于宏观弹性变形区和局部损伤破坏区,能够准确、高效率地分析这些复杂的力学问题[223]究[426]。,。本文采用力学与动力学参数传递以及权值分配法,建立了离散元与有限元耦合的时空多尺度计算方法,既能全面实现空间多尺度,同时又兼顾时间多尺度,有效地提高了多尺度数值计算的效率。同时,作者将这一时空多尺度计算方法应用于激光辐照下受拉铝板破坏行为的数值模拟研究中,取得了较好的效果。
式中:Δtd为离散元时间步长;Δtdcr为离散元临界时间步长;k为时间步数,k>1;i和j为两个相邻GP
元;Δtdcr为离散元几何限制临界时间步长;Δtdcr为离散元物理限制临界时间步长。
由于离散元法采用显式差分格式求解,为了耦合两种方法实现多尺度计算为了保证有限,[10]
ftfcr=
(1+ν)(1-2ν)
(2)
式中:Δtf为有限元时间步长;Δtfcr为有限元临界时间步长;Lmin为模型中网格单元最小尺寸;ρ为
材料密度;E为杨氏模量;ν为泊松比。
从离散元与有限元时间步长的限制条件可知,离散元时间步长取决于元尺寸、邻居元距离以及元性质,所以离散元时间步长是很小的。有限元时间步长不仅与结构或材料的性质有关,而且与模型的网格单元最小尺寸成正比。因此,为了实现时间多尺度计算,离散元采用固定较小时间步长的同时,有限元需要采用较大的时间步长,这时根据式(2
)就需要模型具有较大的有限元网格单元最小尺寸。
1.2 时空多尺度计算模型
1 计算模型与方法
1.1 基本理论
离散元法与有限元法具有相似形式的求解方程,都可以归结为一般的牛顿运动定律形式,课题组曾根据这一理论基础,通过引入过渡层概念,对元与节点重叠的有限元和离散元耦合空间多尺度
方法(以下简称空间多尺度法)进行了研究[728]。研究发现,这种多尺度法只部分实现了空间多尺度,同时尚未实现时间多尺度计算,数值模拟仍需要较多计算时间,多尺度高效率优势尚未得到较好发挥。
离散元法与有限元法的时间步长直接影响时间多尺度计算的实现。离散元法的时间步长存在几何和物理两方面的限制[9]。几何方面,在一个时间步长内不允许一个元与其邻居元沿它们的中心轴线相互贯穿对方;物理方面,在一个时间步长内所产生的动量转移不应超过整个碰撞过程中全部动量的转移。离散元法计算的临界时间步长应当取两者中的最小值。另外,为了提高碰撞过程的时间分辨率,计算中应该选取大于1的时间步数,这样离散元数值模拟的时间步长限制为
在空间多尺度法的基础上,基于以上思想作者建立了既能有效实现时间多尺度,同时又能完善空间多尺度的离散元与有限元耦合的时空多尺度计算方法,这种方法的计算模型如图1所示。
从图中可以看出,在耦合过渡区中,有限元网格单元的最小尺寸不必受到元的大小及元间距的
图1 离散元与有限元耦合的时空多尺度计算模型
Fig.1 Multiscalenumericalmodelofspaceand
timebycouplingDEMandFEM
・410・
吉林大学学报(工学版)第39卷
限制,这样网格单元的最小尺寸可以取得较大,从而增大有限元时间步长,而这时离散元可以采用较小的时间步长,完成时间多尺度计算。另外,该模型中整体有限元网格的数量受耦合过渡区网格尺度的影响较小,具备完善的空间多尺度计算能力。
在计算模型的耦合过渡区中,离散元和有限元相互为对方提供边界条件,并采用权值分配法进行力学和动力学参数的传递。也就是说,离散元为有限元提供力边界条件采用式(3)进行计算,)进行计算。,)相应变量中n
q
(m)i
F
(D)FE(D)FE
=
i=1
∑k
F
i(m)DE
+
∑∑k
p
l
(l)i
FDE
i(l)
(3)
式中:F
为离散元提供力边界条件后的网格节
i(m)DE
点作用力;ki(m)为第m网格内第i元在该网格节点的权值;F
为第m网格内第i元受到的合
(l)i
图2 时空多尺度算法流程图
Fig.2 Flowchartofthemultiscalealgorithm
oftimeandspace
力;n表示第m网格包含的元数目;k
为第l(l
i(l)
m)网格内第i元在该网格节点的权值;FDE为第l(l
2 模型验证和有效性评价
用弹性应力波传播算例来验证时空多尺度计算方法的有效性和优越性。选用一根长度、宽度和厚度分别为3.2×10-1m、2.4×10-3m和3.2×10-4m的细长杆。杆的右端固定,左端施加一个突然加载和卸载的分布载荷,载荷的大小为1.0×104N/m,脉宽为1.5×10-5s。杆的材料为
含的元数目;p为包含该网格节点的所有l(l
m)网格的数目。
3
V
(F)DE
=
i=1
∑kV
i
iFE(4)
(F)
式中:VDE为有限元提供速度边界条件后的元速
度;ki为该元所属网格的第i节点(三角形网格,i=1,2,3)在该元位置处的权值;ViFE为该元所属
网格单元的第i节点的速度。
为了实现耦合过渡区材料密度的匹配性,需要使过渡区内网格节点的质量与相应元过渡到该节点的质量相等,计算方法为
n
q
铝,密度为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所示。
M
(D)
FE
=
i=1
∑
k
(m)i
M
i(m)DE
+
∑∑
p
(l)i(l)
kiMDE(5)
l
式中:M
M
i(m)DE
(D)FE为元质量传递后的网格节点质量;
i(l)DE
为第m网格内第i元的质量;M为第l(l
在时空多尺度数值计算中,离散元时间步长设为dt,有限元时间步长设为fdt,两者关系为
fdt/dt=N。也就是说,在每个有限元时间步长
图4为使用单一有限元法、时空多尺度法以及解析法时应力波在细长杆A、B、C三个位置的
求解计算中,执行N次离散元计算过程来完成时
间多尺度计算。整个系统的计算流程如图2所示
。
图3 细长杆的时空多尺度计算模型
Fig.3 Multiscalenumericalmodeloftimeand
spaceoftheslenderpole
第2期张 锐,等:
离散元与有限元耦合的时空多尺度计算方法
・411・
随着耦合过渡区内离散元层数的增多以及有限元时间步长的增大,计算时间明显减小。并且,时空
多尺度法所有情况下的计算时间均明显少于空间多尺度法的计算时间,当耦合过渡区内布置五层离散元时,采用时空多尺度法的计算时间仅相当于空间多尺度法计算时间的1/3。可见,时空多尺度法在节省计算时间方面具有较大的优越性,。
图4 不同方法计算的细长杆应力波传播情况
Fig.4 Stress2waveby[11]。本文使用离散元与有限元耦合的时空多尺度计算方法对激光辐照下受拉铝板破坏行为进行数值模拟。选用长度、宽度和厚度分别为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为激光辐照下受拉铝板破坏过程的时空多尺度数值模拟图,图中的黑色三角形标出微裂纹的位置。从图中可以看出,微裂纹首先出现在
传播情况,限元法相比,情况基本一致,应力波能够在有限元区(A和C位置)和离散元区(B位置)之间光滑过渡。同时,时空多尺度法计算的应力波到达三个位置的时间及峰值与解析解基本吻合。以上结果验证了离散元与有限元耦合的时空多尺度计算方法的有效性。由于引入了阻尼,波在传播过程中有一定的衰减。由于杆的横向惯性效应,波形将出现剧烈振荡,图4中的模拟结果很好地说明了这一点。
为了研究时空多尺度方法的计算效率,将空间多尺度法与耦合过渡区内布设不同层离散元的时空多尺度法在弹性应力波传播算例中的计算时间进行了对比,如表1所示。根据式(2)可知,有限元临界时间步长与网格单元最小尺寸成正比。为了对比分析方便,在时空多尺度法中,耦合过渡区网格尺寸(不大于非耦合过渡区网格最小尺寸)为离散单元间距的倍数,令耦合过渡区内网格尺寸与元间距相等时,有限元与离散元的时间步长相等,随着元层数增加,有限元时间步长相应增倍。从表1中可以看出,当采用时空多尺度法时,
表1 不同多尺度法的计算时间
Table1 Calculationtimebydifferentmultiscalemethods
空间多尺度
离散元数离散元时间步长/s计算时间步数有限元网格数
耦合过渡区内离散元层数有限元与离散元时间步长比计算时间/s
最小与最大时间比值
845.0×10-921000536
时空多尺度
845.0×10-921000472
二
1127.210.36
三
256.97
四
350.39
五
445.83
图5 激光辐照受拉铝板的时空多尺度数值计算模型
Fig.5 Multiscalenumericalmodeloftimeandspaceof
atensilealuminumplateunderlaserirradiation
・412・
吉林大学学报(工学版)第39
卷
[3]MunjizaA,BangashT,JohnNWM.Thecombined
finite2discreteelementmethodforstructuralfailureandcollapse[J].EngineeringFractureMechanics,2004,71:4692483.
[4]傅华,刘仓理,王文强,等.冲击动力学中离散元与
有限元相结合的计算方法研究[J].高压物理学报,
2006,20(4):3792385.
FuHua,LiuCang2li,Wang2qiang,etal.Adiscrete/elementshockdynam2
图6 激光辐照受拉铝板破坏过程的时空
多尺度数值模拟图
Fig.6 Multiscaleindamage
of].ofPressurePhysics,,20(:379[YT.Parallelisedfinite/discrete
simulationofmulti2fracturingsolidsanddis2cretesystems[J].EngineeringComputations,2001,18(3/4):5572576.
[6]MunjizaA,JohnNWM.Meshsizesensitivityof
thecombinedFEM/DEMfractureandfragmentationalgorithms[J].EngineeringFractureMechanics,2002,69:2812295.
[7]胥建龙,唐志平.离散元与有限元结合的多尺度方
点附近,随后以该微裂纹处作为裂源,在其周围萌生大量新微裂纹,并沿横向迅速向铝板边缘和光斑中心双向扩展。微裂纹的扩大将形成宏观裂纹,并在很短的时间内沿水平方向贯穿整个铝板,直至完全断裂,这些模拟结果与文献[11]提供的实验结果基本一致。
法及其应用[J].计算物理,2003,20(6):4772482.
XuJian2long,TangZhi2ping.Combineddiscrete/fi2niteelementmultiscalenumericalmethodanditsap2plication[J].
ChineseJournalofComputational
Physics,2003,20(6):4772482.
[8]TangZhi2ping,XuJian2long.AcombinedDEM/
FEMmultiscalemethodandstructurefailuresimu2lationunderlaserirradiation[C]∥ShockCompres2sionofCondensedMatter22005:ProceedingsoftheConferenceoftheAmericanPhysicalSocietyTopicalGrouponShockCompressionofCondensedMatter,Balitimore,Maryland,USA,2005.
[9]唐志平.三维离散元方法及其在冲击力学中的应用
[J].中国科学(E辑),2003,33(11):9892998.TangZhing2ping.32DimensionDEMtheoryanditsapplicationtoimpactmechanics[J].ScienceinChina(SeriesE),2003,33(11):9892998.
[10]钱东升,华林,左治江,等.环件轧制三维有限元模
4 结束语
本文建立了离散元与有限元耦合的时空多尺度计算方法,通过弹性应力波传播算例表明,该方法能够满足精确性要求,并能并行实现时间多尺度与空间多尺度计算,有效地提高了多尺度数值计算效率。同时,本文将该多尺度方法应用于激光辐照下受拉铝板破坏行为的数值模拟中,模拟结果与实验结果基本一致。时空多尺度计算方法通过有限元区域采用较大时间步长来而离散元区域采用较小时间步长来实现时间多尺度,并通过增大耦合过渡区内有限元网格尺寸来优化空间多尺度。因此,对于包含大型有限元区域的复杂结构力学问题,该方法更能发挥高效率计算的优势,具有重要的工程应用前景。参考文献:
[1]ParkHS,LiuWK.Anintroductionandtutorialon
multiple2scaleanalysisinsolids[J].2004,193:173321772.
[2]AzevedoNM,LemosJV.Hybriddiscreteele2
ment/finiteelementmethodforfractureanalysis[J].ComputerMethodsinAppliedMechanicsandEngi2neering,2006,195:457924593.
Computer
MethodsinAppliedMechanicsandEngineering,
拟中质量缩放方法的运用[J].塑性工程学报,
2005,12(5):86291,100.
QianDong2sheng,HuaLin,ZuoZhi2jiang,etal.Applicationofmassscalinginsimulationofringrollingbythree2dimensionalfiniteelementmethod[J].JournalofPlasticityEngineering,2005,12(5):86291,100.
[11]孙承伟.激光辐照效应[M].北京:国防工业出版
社,2002.
第39卷 第2期2009年3月
吉林大学学报(工学版)
JournalofJilinUniversity(EngineeringandTechnologyEdition)
Vol.39 No.2 Mar.2009
离散元与有限元耦合的时空多尺度计算方法
张 锐1,唐志平2,郑 航2
(1.吉林大学地面机械仿生技术教育部重点实验室长春130022;2.中国科学技术大学中国科学院材料力学
行为和设计重点实验室合肥230027)
摘 要:采用力学与动力学参数传递以及权值分配法度计算方法。度双重优势,,,将这一多尺度,得到的模拟结果与实验结果基本一致关键词:;时空多尺度;离散元与有限元耦合;数值计算中图分类号:O242.1 文献标识码:A 文章编号:167125497(2009)0220408205
Timeandspacemultiscalenumericalmethodbycoupling
discreteelementmethodandfiniteelementmethod
ZHANGRui1,TANGZhi2ping2,ZHENGHang2
(1.KeyLaboratoryofTerrain2machineBionicsEngineering,MinistryofEducation,JilinUniversity,Changchun
130022,China;2.KeyLaboratoryforMechanicalBehaviorandDesignofMaterials(LMBD)ofCAS,UniversityofScienceandTechnologyofChina,Hefei230027,China)
Abstract:Throughmechanicalanddynamicparameterstransitionandweightvaluedistribution,atimeandspacemultiscalenumericalmethodbycouplingDiscreteElementMethod(DEM)andFiniteElementMethod(FEM)wasproposed.Simulationofanexampleofelasticstress2wavepropagationwascarriedout.
Theresultsshowthatthemultiscalenumericalmethodpossessesthedual
superioritiesofmultiscaletimeandmultiscalespace,whichsignificantlyimprovetheefficiencyofsimulation.Itisshownthatthesimulationresultsbythemethodsatisfytherequirementofaccuracy.Themethodisalsoemployedtosimulatethedamagebehaviorsofatensilealuminumplateunderlaserirradiation,andtheresultsareaccordantwiththeexperimentalresults.
Keywords:numericalmathematics;multiscaletimeandspace;couplingdiscreteelementmethodandfiniteelementmethod;numericalcomputation
许多复杂力学问题,例如裂纹尖端应力场问题、结构在热力耦合载荷作用下局部破坏问题等,
收稿日期:2007208222.
难以用单一尺度数值方法得到满意的模拟结果。发展多尺度计算,针对不同区域采用不同的数值
基金项目:国家自然科学基金项目(50805064,10472114);吉林省科技发展计划项目(20080146);2006年度中国科学
院王宽诚博士后工作奖励基金项目.
作者简介:张锐(19752),男,讲师,博士.研究方向:数值计算.E2mail:[email protected]
第2期张 锐,等:离散元与有限元耦合的时空多尺度计算方法
・409・
(1)
计算方法能够有效解决这类问题,现已成为当今数值模拟研究的一个热点
[1]
。在固体结构损伤和
GP
Δtd≤Δtdcr=min{Δtdcr,Δtdcr}
k
i,j
破坏的力学问题中,充分考虑离散元法和有限元法的优势和弊端,即离散元法能够方便地呈现局部细观非均匀性,但需要大量计算时间;有限元法善于快速求解整体结构宏观弹性变形,但难以准确处理局部结构细微损伤、断裂等力学现象,将这两种方法耦合实现多尺度计算,分别应用于宏观弹性变形区和局部损伤破坏区,能够准确、高效率地分析这些复杂的力学问题[223]究[426]。,。本文采用力学与动力学参数传递以及权值分配法,建立了离散元与有限元耦合的时空多尺度计算方法,既能全面实现空间多尺度,同时又兼顾时间多尺度,有效地提高了多尺度数值计算的效率。同时,作者将这一时空多尺度计算方法应用于激光辐照下受拉铝板破坏行为的数值模拟研究中,取得了较好的效果。
式中:Δtd为离散元时间步长;Δtdcr为离散元临界时间步长;k为时间步数,k>1;i和j为两个相邻GP
元;Δtdcr为离散元几何限制临界时间步长;Δtdcr为离散元物理限制临界时间步长。
由于离散元法采用显式差分格式求解,为了耦合两种方法实现多尺度计算为了保证有限,[10]
ftfcr=
(1+ν)(1-2ν)
(2)
式中:Δtf为有限元时间步长;Δtfcr为有限元临界时间步长;Lmin为模型中网格单元最小尺寸;ρ为
材料密度;E为杨氏模量;ν为泊松比。
从离散元与有限元时间步长的限制条件可知,离散元时间步长取决于元尺寸、邻居元距离以及元性质,所以离散元时间步长是很小的。有限元时间步长不仅与结构或材料的性质有关,而且与模型的网格单元最小尺寸成正比。因此,为了实现时间多尺度计算,离散元采用固定较小时间步长的同时,有限元需要采用较大的时间步长,这时根据式(2
)就需要模型具有较大的有限元网格单元最小尺寸。
1.2 时空多尺度计算模型
1 计算模型与方法
1.1 基本理论
离散元法与有限元法具有相似形式的求解方程,都可以归结为一般的牛顿运动定律形式,课题组曾根据这一理论基础,通过引入过渡层概念,对元与节点重叠的有限元和离散元耦合空间多尺度
方法(以下简称空间多尺度法)进行了研究[728]。研究发现,这种多尺度法只部分实现了空间多尺度,同时尚未实现时间多尺度计算,数值模拟仍需要较多计算时间,多尺度高效率优势尚未得到较好发挥。
离散元法与有限元法的时间步长直接影响时间多尺度计算的实现。离散元法的时间步长存在几何和物理两方面的限制[9]。几何方面,在一个时间步长内不允许一个元与其邻居元沿它们的中心轴线相互贯穿对方;物理方面,在一个时间步长内所产生的动量转移不应超过整个碰撞过程中全部动量的转移。离散元法计算的临界时间步长应当取两者中的最小值。另外,为了提高碰撞过程的时间分辨率,计算中应该选取大于1的时间步数,这样离散元数值模拟的时间步长限制为
在空间多尺度法的基础上,基于以上思想作者建立了既能有效实现时间多尺度,同时又能完善空间多尺度的离散元与有限元耦合的时空多尺度计算方法,这种方法的计算模型如图1所示。
从图中可以看出,在耦合过渡区中,有限元网格单元的最小尺寸不必受到元的大小及元间距的
图1 离散元与有限元耦合的时空多尺度计算模型
Fig.1 Multiscalenumericalmodelofspaceand
timebycouplingDEMandFEM
・410・
吉林大学学报(工学版)第39卷
限制,这样网格单元的最小尺寸可以取得较大,从而增大有限元时间步长,而这时离散元可以采用较小的时间步长,完成时间多尺度计算。另外,该模型中整体有限元网格的数量受耦合过渡区网格尺度的影响较小,具备完善的空间多尺度计算能力。
在计算模型的耦合过渡区中,离散元和有限元相互为对方提供边界条件,并采用权值分配法进行力学和动力学参数的传递。也就是说,离散元为有限元提供力边界条件采用式(3)进行计算,)进行计算。,)相应变量中n
q
(m)i
F
(D)FE(D)FE
=
i=1
∑k
F
i(m)DE
+
∑∑k
p
l
(l)i
FDE
i(l)
(3)
式中:F
为离散元提供力边界条件后的网格节
i(m)DE
点作用力;ki(m)为第m网格内第i元在该网格节点的权值;F
为第m网格内第i元受到的合
(l)i
图2 时空多尺度算法流程图
Fig.2 Flowchartofthemultiscalealgorithm
oftimeandspace
力;n表示第m网格包含的元数目;k
为第l(l
i(l)
m)网格内第i元在该网格节点的权值;FDE为第l(l
2 模型验证和有效性评价
用弹性应力波传播算例来验证时空多尺度计算方法的有效性和优越性。选用一根长度、宽度和厚度分别为3.2×10-1m、2.4×10-3m和3.2×10-4m的细长杆。杆的右端固定,左端施加一个突然加载和卸载的分布载荷,载荷的大小为1.0×104N/m,脉宽为1.5×10-5s。杆的材料为
含的元数目;p为包含该网格节点的所有l(l
m)网格的数目。
3
V
(F)DE
=
i=1
∑kV
i
iFE(4)
(F)
式中:VDE为有限元提供速度边界条件后的元速
度;ki为该元所属网格的第i节点(三角形网格,i=1,2,3)在该元位置处的权值;ViFE为该元所属
网格单元的第i节点的速度。
为了实现耦合过渡区材料密度的匹配性,需要使过渡区内网格节点的质量与相应元过渡到该节点的质量相等,计算方法为
n
q
铝,密度为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所示。
M
(D)
FE
=
i=1
∑
k
(m)i
M
i(m)DE
+
∑∑
p
(l)i(l)
kiMDE(5)
l
式中:M
M
i(m)DE
(D)FE为元质量传递后的网格节点质量;
i(l)DE
为第m网格内第i元的质量;M为第l(l
在时空多尺度数值计算中,离散元时间步长设为dt,有限元时间步长设为fdt,两者关系为
fdt/dt=N。也就是说,在每个有限元时间步长
图4为使用单一有限元法、时空多尺度法以及解析法时应力波在细长杆A、B、C三个位置的
求解计算中,执行N次离散元计算过程来完成时
间多尺度计算。整个系统的计算流程如图2所示
。
图3 细长杆的时空多尺度计算模型
Fig.3 Multiscalenumericalmodeloftimeand
spaceoftheslenderpole
第2期张 锐,等:
离散元与有限元耦合的时空多尺度计算方法
・411・
随着耦合过渡区内离散元层数的增多以及有限元时间步长的增大,计算时间明显减小。并且,时空
多尺度法所有情况下的计算时间均明显少于空间多尺度法的计算时间,当耦合过渡区内布置五层离散元时,采用时空多尺度法的计算时间仅相当于空间多尺度法计算时间的1/3。可见,时空多尺度法在节省计算时间方面具有较大的优越性,。
图4 不同方法计算的细长杆应力波传播情况
Fig.4 Stress2waveby[11]。本文使用离散元与有限元耦合的时空多尺度计算方法对激光辐照下受拉铝板破坏行为进行数值模拟。选用长度、宽度和厚度分别为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为激光辐照下受拉铝板破坏过程的时空多尺度数值模拟图,图中的黑色三角形标出微裂纹的位置。从图中可以看出,微裂纹首先出现在
传播情况,限元法相比,情况基本一致,应力波能够在有限元区(A和C位置)和离散元区(B位置)之间光滑过渡。同时,时空多尺度法计算的应力波到达三个位置的时间及峰值与解析解基本吻合。以上结果验证了离散元与有限元耦合的时空多尺度计算方法的有效性。由于引入了阻尼,波在传播过程中有一定的衰减。由于杆的横向惯性效应,波形将出现剧烈振荡,图4中的模拟结果很好地说明了这一点。
为了研究时空多尺度方法的计算效率,将空间多尺度法与耦合过渡区内布设不同层离散元的时空多尺度法在弹性应力波传播算例中的计算时间进行了对比,如表1所示。根据式(2)可知,有限元临界时间步长与网格单元最小尺寸成正比。为了对比分析方便,在时空多尺度法中,耦合过渡区网格尺寸(不大于非耦合过渡区网格最小尺寸)为离散单元间距的倍数,令耦合过渡区内网格尺寸与元间距相等时,有限元与离散元的时间步长相等,随着元层数增加,有限元时间步长相应增倍。从表1中可以看出,当采用时空多尺度法时,
表1 不同多尺度法的计算时间
Table1 Calculationtimebydifferentmultiscalemethods
空间多尺度
离散元数离散元时间步长/s计算时间步数有限元网格数
耦合过渡区内离散元层数有限元与离散元时间步长比计算时间/s
最小与最大时间比值
845.0×10-921000536
时空多尺度
845.0×10-921000472
二
1127.210.36
三
256.97
四
350.39
五
445.83
图5 激光辐照受拉铝板的时空多尺度数值计算模型
Fig.5 Multiscalenumericalmodeloftimeandspaceof
atensilealuminumplateunderlaserirradiation
・412・
吉林大学学报(工学版)第39
卷
[3]MunjizaA,BangashT,JohnNWM.Thecombined
finite2discreteelementmethodforstructuralfailureandcollapse[J].EngineeringFractureMechanics,2004,71:4692483.
[4]傅华,刘仓理,王文强,等.冲击动力学中离散元与
有限元相结合的计算方法研究[J].高压物理学报,
2006,20(4):3792385.
FuHua,LiuCang2li,Wang2qiang,etal.Adiscrete/elementshockdynam2
图6 激光辐照受拉铝板破坏过程的时空
多尺度数值模拟图
Fig.6 Multiscaleindamage
of].ofPressurePhysics,,20(:379[YT.Parallelisedfinite/discrete
simulationofmulti2fracturingsolidsanddis2cretesystems[J].EngineeringComputations,2001,18(3/4):5572576.
[6]MunjizaA,JohnNWM.Meshsizesensitivityof
thecombinedFEM/DEMfractureandfragmentationalgorithms[J].EngineeringFractureMechanics,2002,69:2812295.
[7]胥建龙,唐志平.离散元与有限元结合的多尺度方
点附近,随后以该微裂纹处作为裂源,在其周围萌生大量新微裂纹,并沿横向迅速向铝板边缘和光斑中心双向扩展。微裂纹的扩大将形成宏观裂纹,并在很短的时间内沿水平方向贯穿整个铝板,直至完全断裂,这些模拟结果与文献[11]提供的实验结果基本一致。
法及其应用[J].计算物理,2003,20(6):4772482.
XuJian2long,TangZhi2ping.Combineddiscrete/fi2niteelementmultiscalenumericalmethodanditsap2plication[J].
ChineseJournalofComputational
Physics,2003,20(6):4772482.
[8]TangZhi2ping,XuJian2long.AcombinedDEM/
FEMmultiscalemethodandstructurefailuresimu2lationunderlaserirradiation[C]∥ShockCompres2sionofCondensedMatter22005:ProceedingsoftheConferenceoftheAmericanPhysicalSocietyTopicalGrouponShockCompressionofCondensedMatter,Balitimore,Maryland,USA,2005.
[9]唐志平.三维离散元方法及其在冲击力学中的应用
[J].中国科学(E辑),2003,33(11):9892998.TangZhing2ping.32DimensionDEMtheoryanditsapplicationtoimpactmechanics[J].ScienceinChina(SeriesE),2003,33(11):9892998.
[10]钱东升,华林,左治江,等.环件轧制三维有限元模
4 结束语
本文建立了离散元与有限元耦合的时空多尺度计算方法,通过弹性应力波传播算例表明,该方法能够满足精确性要求,并能并行实现时间多尺度与空间多尺度计算,有效地提高了多尺度数值计算效率。同时,本文将该多尺度方法应用于激光辐照下受拉铝板破坏行为的数值模拟中,模拟结果与实验结果基本一致。时空多尺度计算方法通过有限元区域采用较大时间步长来而离散元区域采用较小时间步长来实现时间多尺度,并通过增大耦合过渡区内有限元网格尺寸来优化空间多尺度。因此,对于包含大型有限元区域的复杂结构力学问题,该方法更能发挥高效率计算的优势,具有重要的工程应用前景。参考文献:
[1]ParkHS,LiuWK.Anintroductionandtutorialon
multiple2scaleanalysisinsolids[J].2004,193:173321772.
[2]AzevedoNM,LemosJV.Hybriddiscreteele2
ment/finiteelementmethodforfractureanalysis[J].ComputerMethodsinAppliedMechanicsandEngi2neering,2006,195:457924593.
Computer
MethodsinAppliedMechanicsandEngineering,
拟中质量缩放方法的运用[J].塑性工程学报,
2005,12(5):86291,100.
QianDong2sheng,HuaLin,ZuoZhi2jiang,etal.Applicationofmassscalinginsimulationofringrollingbythree2dimensionalfiniteelementmethod[J].JournalofPlasticityEngineering,2005,12(5):86291,100.
[11]孙承伟.激光辐照效应[M].北京:国防工业出版
社,2002.