粘滞材料冲击碰撞应力场的数值模拟

粘滞材料冲击碰撞应力场的数值模拟

汽车安全的数学模型是计算和建模的重要课题之一。它可以为中国快速、国际先进的汽车行业提供平台和参考。例如，弹塑材料的大变形模拟是当前科学工程计算中非常重要和紧迫的问题。它涉及飞机、火车、船只、大型建筑物等安全结构设计的概念。因此，这项研究不仅具有理论上的重要性和紧迫性，而且具有直接和巨大的经济效益。冯康、石钟慈和林群等以战略眼光对弹性结构的数学模型及其数值解法创立了系统的研究方法,从泛函分析理论高度对计算方法的高精度收敛性进行了开拓性的工作,取得了国际领先的结果.经典的Maxwell方程计算已被广泛运用于电磁场理论的数学建模,并被推广至磁流动力学计算.而近来科学家们又发现Maxwell方程可用于软固体材料(如尼龙、橡胶制品)等制作流程的流动应力场计算.PhanThien及Tanner(1987年)在Maxwell方程应用与实验的基础上建立了具有Elongation特征的应力计算,从而对粘滞弹塑性流体的正压力与切变应力作出了更全面的描述.剑桥大学科学家L.J.Gibson和M.F.Ashby(1988年)对蜂窝多孔材料力学性质作了系统的论述.EEVC、Euro-NCAP、IIHS、NHTSA等国际汽车安全法规机构对蜂窝多孔材料在汽车碰撞中的应用作了大量实验,并由此而制成了相关安全法规(1995～2005年).中国汽车安全技术研究部门,如清华大学等也开展了大量汽车碰撞试验数据采集系统的研制与公路交通FEA模拟(1990～2006年).然而,对冲击碰撞问题的信息压缩至微分方程数学模型的研究工作还尚待完善.本工作将借鉴Lagrange-Euler方法计算网格的移动及应力场变化,运用带有指数冲击项的Maxwell方程与描述大变形的Cauchy方程来模拟弹塑性大变形的流固耦合问题.通过进一步研究流变学,完善优化非线性自适应有限元计算方法.对汽车安全(EEVC安全法规)所采用的模型进行深入学习并建立具有新思想和新数学模型.在此基础上运用奇异摄动方法讨论扰动参数对自适应有限元方法的适定性问题.1类变量p-t/t的数值响应面为—弹塑性材料的冲击碰撞问题数学模型我们运用流变学理论对各向异性弹塑性P-T/T应力方程(Maxwell方程指数增值冲击项)来计算应力场分布变化律.式中,τ是应力场,D是应变,u是弹塑性材料的变形速度场,η是粘滞度,λ是弹塑性松弛量,ε是延长因子,ξ是切变因子,ρ是材料密度.其Maxwell微分方程FEA计算收敛性至少具2阶精度.另外,我们计算了由应力场分布变化˙ττ˙产生的弹塑性材料的大变形,包括受冲击边界以外Ω内u的Cauchy守衡方程:ρ˙u=[■⋅τ-ρu⋅■u]‚inΩ‚(2)初始值由静态试验结果τ(0)和试验的冲击速度u(0)决定,记为:τ(0)=τ(static),u(0)=ux0onΓ0∈Γ,其中Ω是材料体积,Γ是材料表面,Γ0是受冲击表面(沿x方向移动边界).以上应力方程(1)的离散与计算及稳定性分析已在作者2001年的论文中详细描述.有关工程实验论证及参考见以下结果(AE-MDB/IIHS/NHTSA),其中的关键问题是Maxwell方程的移动边界条件与有限元自适应算法.另外混合元等方法也非常有效.为解决方程(2)中非线性对流项引出的著名的Kelvin-Helmholtz不稳定切变流现象,作者使用分段二次试验函数即9点双二次元.此类非线性逼近的超收敛性和有效性在理论与实践上得到了验证.由加权剩余法得到的有关方程(1)、(2)的FEA逼近在自适应网格点上求解并对大变形区域中各向异性网格的加密进行智能判别与控制.我们采用Galerkin方法9点双二次元计算.网格的移动变形由Euler网格跟踪Cauchy方程的流场,然后由Lagrange网格计算Maxwell方程的流动应力场.当计算二相流时二网格存在流速差,否则二网格重合.我们在有关汽车材料部门合作下完成了数学模拟的FEA非线性分析与实验论证.从而开始构造用冲击块(impactbarrier)代替冲击试验汽车并实现微分方程理论在碰撞中的应用.我们模拟的FEA非线性结果(三维)将在第2节中显示.其x方向冲击加速度G数值解与实验吻合,由此而得到材料的大变形应力场分布.我们将在此基础上,运用自适应Lagrange-Euler方法建立高精度网格变形算法,进一步完善与有限元冲击接触问题的求解格式相容互补的自适应格式.有关Maxwell方程及P-T/T方程的收敛性与稳定性在文献中已对标准化正定格式的建立作了详细证明,并保证了方程的有限元求解至少具有2阶精度.2由弹塑料制成的影响冲突模型为结果2.1[n+y+3.]uxy+uy+y+y+y+y+y+y+y+y+y+3.2.0.FEA非线性结果确定了固体材料变形与非接触边界.由此得到的应力场分布及变形边界的速度与微分方程理论结果完全符合(见图1、2及文献).对应的P-T/T方程在frontalsolver中的预估-校正格式的后估计式为{τn+1xx}={τnxx}+{J1}n+{J1}n+1pre‚(3)式中,J1=Δtηλ∂ux∂x-Δt2λτxx-Δtε2η0(τxx+τyy)τxx-Δt2(ux∂∂x+uy∂∂y)τxx+Δt(τxx∂∂x+τxy∂∂y)ux-Δtξ2[2τxx∂ux∂x+τxy(∂ux∂y+∂uy∂x)].{τn+1xy}={τnxy}+{J2}n+{J2}n+1pre‚(4)式中,J2=Δtη2λ(∂ux∂y+∂uy∂x)-Δt2λτxy-Δtε2η0(τxx+τyy)τxy-Δt2(ux∂∂x+uy∂∂y)τxy+Δt2[(∂ux∂x+∂uy∂y)τxy+τxx∂ux∂y+τyy∂uy∂x-Δtξ212(∂ux∂y+∂uy∂x)⋅(τxx+τyy)+τxy(∂ux∂x+∂uy∂y)].{τn+1yx}={τn+1xy},(5){τn+1yy}={τnyy}+{J3}n+{J3}n+1pre,(6)式中,J3=Δtηλ∂uy∂y-Δt2λτyy-Δtε2η0(τxx+τyy)τyy-Δt2(ux∂∂x+uy∂∂y)τyy+Δt(τxy∂∂x+τyy∂∂y)uy-Δtξ2[2τyy∂uy∂y+τxy(∂ux∂y+∂uy∂x)]‚式中,{·}n+1pre是迭代{·}n的预估解已在文献中给出,初值(n=0)由实验得到.Cauchy方程在frontalsolver中采用9点双二次元(9pointbi-quadraticelement)的离散形式{Un+1x}9={Unx}9+Δtρ{(∂∂xτxx+∂∂yτxy)D-ρ[{Unx}⋅{C1}+{Unx}⋅{C2}]}9⋅{A-1}‚(7){Un+1y}9={Uny}9+Δtρ{(∂∂xτxy+∂∂yτyy)D-ρ[{Uny}⋅{C1}+{Uny}⋅{C2}]}9⋅{A-1}.(8)式中,Ux、Uy是弹塑性材料的变形速度分量,A、C、D是刚度矩阵,沿x方向移动边界初值(n=0)条件由实验得到.由此降阶处理得到有效格式,节省了大量计算.2.2有限元数值分析上述非线性显式模型在冲击碰撞问题中具有非常有效的解题功能.图1为汽车碰撞试验(test)与有限元模拟初步模型(CAE)的比较,可归纳为4种情况.在此假定,除冲击板块为可塑性材料外,冲击车、障碍墙及B柱均为不变形钢材.(1)非奇异碰撞区域.低中速冲击(u0:3～9m/s)包括表面充分光滑的钢体冲击物.无需局部网格加密,自适应方法待机,非线性代数方程组具有强收敛解.冲击特征曲线(图2(a))与冲击板块变形有很好的光滑性.(2)非奇异碰撞区域.高速冲击(u0>10m/s)包括表面充分光滑的钢体冲击物.自适应方法网格加密非常有效.非线性代数方程组具有强收敛解;冲击特征曲线(图2(b))与冲击板块变形有很好的光滑性.(3)奇异碰撞区域.低中速冲击(u0:3～9m/s)包括表面尖角的钢体冲击物.自适应方法网格加密有效但收敛速度降低,非线性代数方程组具有弱收敛解.(4)由于局部元素变形达破坏极限,自适应方法网格加密失效.冲击特征曲线(图2(c))有明显局部偏差,需要指数非线性方程求解.我们采用有限元前处理器中MAT-honeycomb对多孔材料极限元素处理方法,用实验所得参数进行断裂适应性模拟,以达到最佳材料极限特性.但是材料的微观与冲击极限的宏观特性拟合仍需进行更进一步的流变学研究,以完善此建模课题.在以上假设下,由接触表面Γ0决定的边界条件为τ(0)=τ(static),u(0)=u0onΓ0∈Γ.质量约1000kg的汽车冲击力小于F=M·a=250kN.冲击车、障碍板及B柱均为不变形钢材.这与实际情况一致(图2).若方程组(3)～(6)中ξ→0,则体积变形以正压缩为主;若ε→0,则体积变形以扭曲为主;在图2(a)～(c)中各局部应力场分布均可对应参数ξ、ε的摄动范围.若方程组(7)、(8)中各应力项→0,则体积变形近乎塑性流体,可运用自适应网格点的求解,并对大变形局部区域中各向异性网格的加密进行智能判别与控制.3自适应方法网格加密本工作采用有限元方法对Maxwell方程与Cauchy方程分别进行应力与应变计算,并对弹塑性多孔材料大变形极限问题进行了详细研究.对以流固耦合变形特征的几种自适应方法求解进行讨论.主要结果概述如下:(1)非奇异碰撞区域.低中速冲击(u0:3～9m/s).无需局部网格加密,自适应方法待机非线性代数方程组具有强收敛解.冲击特征曲线与冲击板块变形有很好的光滑性.(2)非奇异碰撞区域.高速冲击(u0:>10m/s)包括表面充分光滑的钢体冲击物.自适应方法网格加密非常有效,非线性代数方程组具有强收敛解.(3)奇异碰撞区域.低中速冲击(u0:3～9m/s)包括表面尖角的钢体冲击物.自适应方法网格加密有效但收敛速度降低,非