基于多学科耦合的内燃机热应力分析

基于多学科耦合的内燃机热应力分析

热应力耦合模拟研究随着颅内工程越来越向高度发展，挖掘机的热量问题越来越突出。内燃机的热负荷通常要从两个方面来考虑:一是受热部件因为温度过高而失效;二是受热部件由于温度梯度过大而产生很大的热应力,以至发生热疲劳而损坏。由于计算机辅助工程的数值模拟理论和计算机硬件条件的限制,早期的内燃机设计主要以经验设计和反复多次的试验检验为主,受热零部件热应力的耦合数值模拟研究很少涉及。近10年来,计算机软硬件技术的飞跃发展使运用CAE技术分析受热零部件的热应力成为可能。在大幅度降低内燃机的设计成本和缩短研发周期的同时,内燃机的性能也得到了更大改进。内燃机系统是一个极其复杂的牵涉到多个学科的耦合系统,其中许多现象是多个物理场相互作用而形成的。为了得到更高精度的数值模拟结果,综合考虑多学科物理现象的耦合变得越来越重要。内燃机热流体-热应力的耦合分析也是基于这种背景下发展起来的。1弱耦合求解方法多学科耦合分析是指考虑了属于不同学科研究范围内两个或多个物理场之间相互作用后对研究对象进行的数值模拟。如考虑流体和固体之间互相作用,求解流场运动与结构变形互相影响的流-固耦合分析;考虑流场速度场与温度场变化引起的热应力变化就是典型的流-热-固(热流体-热应力)耦合分析等等。耦合分析通常有两种解法:一种是强耦合,即需要对包括所有物理方程的矩阵直接求解或者迭代求解;另一种是对多个物理场分别列出其物理方程,按照一定的顺序逐个对各物理场的物理方程进行迭代求解,通过将前一个物理场的计算结果作为载荷施加到下一个物理场的物理方程中进行耦合。这种按物理场分区迭代求解耦合问题的方法称为弱耦合。采用强耦合方式需要对各学科有较深刻的理解,程序的开发成本较高,对硬件的要求也较高。而弱耦合方式可以利用各学科现有的程序,考虑边界上的相互影响。而且,由于对每个物理场的分析都是相对独立的,可以比较容易地采用计算机并行技术解决大型的复杂工程实际问题。鉴于弱耦合求解方式灵活、方便的特点,国外的许多学者和研究机构都在这方面进行了卓有成效的研究,这其中以德国科学计算与算法研究所开发的MPCCI多学科耦合分析工具最为有名,它根据弱耦合原理,开创性地在多学科的商用软件之间建立了功能强大的链接。而国内在这方面的研究还鲜有报导。在此背景下,结合作者多年进行弱耦合分析时积累的经验,以内燃机热流体-热应力分析为例,介绍弱耦合分析时应该注意的一些关键问题,对解决热流体-热应力问题有所帮助。2后一阶段结构分析热流体-热应力弱耦合分析通常包含两个分析阶段,一是热流体分析(CFD)阶段,二是结构分析(FEM)阶段。这里涉及到一个很重要的问题,就是如何将前一阶段分析所得到的温度等作为载荷传递到后一阶段结构分析中去。文中按两个阶段的固体域网格的异同分别进行论述。2.1直接转换方法与热应力分析如果网格完全一致,那么在热应力分析时,经过热流体计算后存贮在每个节点的温度数据可以作为输入数据直接调入结构分析所需的节点上,不需要进行插值,文中称这种方法为直接转换法。如图1所示,CFD分析阶段的计算领域包括了流体和固体领域,而热应力分析阶段只包括固体领域。采用直接转换方法,可以直接将在热流体分析阶段获得的固体部分的温度分布作为负荷输入到热应力分析中去。为了表述简明,以下各图中均略去了流体部分的网格。使用直接转换法的限制是,固体域在热流体分析和结构分析阶段必须采用相同的网格类型。例如,使用以节点为中心的有限体积法建立的模型如果只采用四面体网格的话,就可以使用直接转换法,因为结构分析软件一般都能支持四面体网格。如果结构分析软件能和前一个阶段的热流体计算一样能支持多种网格类型(四面体、六面体、金字塔体和棱柱体),直接转换法也十分方便。2.2网格节点的插值假如两个阶段采用的网格不一致,那么就需要将前一阶段网格节点上温度等作为热负荷插值到后一阶段的网格节点上去。在插值过程中,就有可能产生误差。2.2.1固体域的温度场计算将CFD分析中所得到的温度信息作为负荷传递到热应力分析中通常有两种方法。第一种是考虑得到固体部分的热边界条件后,先求出固体域的温度场,再进一步求出热应力。如图2所示,通过插值求得与流体域接触的固壁表面(粗线所代表的面)的对流换热系数及表面流体温度,然后综合其他热边界条件计算得到固体域的温度场。利用该温度场进行热应力分析,文中称这种方法为面映射法。文献论述了应用该方法计算排气歧管热应力问题的过程。这里需要强调的是,为了得到较准确的对流换热系数及表面流体温度,在划分网格时应注意在流体一侧插入边界层单元。边界层单元的第一层厚度应该尽量保持常数且必须足够精细。图3为边界层网格的处理方法示意图。面映射方法的优点是固体域的温度场可以不在热流体分析阶段求出。因此,可以在热流体分析阶段只划分流体域网格,而在热应力分析阶段只划分固体域网格。这样,在两个阶段分析时都能用较少网格得到结果。2.2.2体映射分析第二种方法就是将流体域和固体域网格一起进行热流体分析,把同时得到的固体域的节点温度值直接插值到热应力分析的有限元网格的节点上去。文中称这种方法为体映射法。图4为采用体映射法示意图,该方法的特点是固体部分网格疏密与类型等均不受前一阶段分析的限制。由于在第二阶段不需要进行温度场的计算,所以比面映射法少一个温度场分析过程。相对直接转换法,体映射法在第二阶段分析时网格单元类型不受前一阶段热流体分析中固体网格类型的限制,显得更灵活。另一方面,目前CFD分析的网格大都是一次线性单元,而结构分析往往为了提高计算精度要求采用二次非线性单元。采用体映射方法完全能够满足不同分析系统对于网格类型的要求。2.3热应力分析中“映射目标”的应用在热流体-热应力分析中,根据映射原理称热流体分析计算结果文件为“映射源”,对应地将热应力分析的计算文件称为“映射目标”。在映射过程中,如果“映射目标”的一个点(如图5中黑色圆点所示)的位置在“映射源”领域内,映射数据由该点周围的点的值内插得到;相反,“映射目标”的一个点的位置在“映射源”领域外,映射数据由离该点最近的点的值外插得到。3中小型结构单元的选用在热流体分析阶段使用SC/Tetra软件,结构热应力分析采用ANSYS软件。这里要强调的是,文中引进的概念和方法适用于任何商业化软件或自编的CFD和有限元分析软件。SC/Tetra是以节点为中心的有限体积法的CFD软件。它支持混合网格,即可同时使用六面体网格、棱柱体网格、金字塔型网格和四面体网格。SC/Tetra与结构分析软件有良好的接口,它的结果可以作为载荷提供给结构分析软件。ANSYS是基于有限元法(FEM)的结构分析软件,拥有丰富的单元库、材料库和出色的求解器,能够求解各类结构问题。在ANSYS里进行热应力分析时,通常有4种单元可以选用,其单元构成如图6～图9所示,括号内注明各单元在ANSYS中的名称。目前的弱耦合解析往往涉及两个独立的软件体系,它们之间单元类型的匹配并不是很容易。比如在ANSYS里,45号单元不支持金字塔型网格;72、92号单元只支持四面体网格。如果在热流体分析中固体部分包含有金字塔型网格,那么进行结构分析时45、72、92号单元都不能使用。换句话说,如果在热流体分析中固体域由混合网格构成的话,95号单元是唯一可以使用的。这些问题在具体分析时经常遇到,尤其要注意。4排气歧管断裂故障的发生机制排气歧管是内燃机的主要受热部件,工作热负荷大,热应力高,工作环境极其恶劣。内燃机在进行耐久性试验时往往会发生排气歧管断裂的故障。由于断裂现象产生的原因较复杂,如单纯热应力引起,排气管与内燃机共振引起,或是两者综合作用引起等。为了深入分析断裂产生的原因,以便排除故障并提出改进措施,很有必要进行排气歧管的热流体-热应力的耦合分析。4.1计算公式的确定采用某内燃机排气歧管作为研究对象,如图10所示。在热流体分析阶段,控制方程采用可压缩N-S方程,湍流模型采用标准k-ε模型,壁面按标准壁面函数计算。计算的曲轴转角范围为720～1440°CA,进出口边界条件由AVLBOOST计算给定,进口按照一个循环内平均流量和平均温度给定,出口给定压力。在SC/Tetra软件里生成由混合网格构成的流体域网格,共230217个单元,74577节点。在热应力分析阶段,不考虑所连接缸盖的热应力,在排气歧管螺孔表面选取一些节点施加3个方向的位移约束,压板处同样选取一些节点施加Z方向位移约束。其他边界条件限于篇幅,不再详述。4.2结果与分析4.2.1采用体映射方法求解结构分析的网格对于直接转换方法,固体域和流体域网格均在SC/Tetra前处理模块中一起划分,因此很难同时兼顾到两个分析阶段对网格的要求,只考虑了流体分析的要求,固体域管道径向仅有1层单元,不能满足结构分析的多层网格的要求;而在面映射法和体映射法中固体域允许采用不同的网格,如可以对固体域的管厚方向划分2层单元,这样既可以保证结构分析的网格要求,又不至于增加热流体分析的负担,如图11所示。图12为采用体映射方法得到的排气歧管的位移分布。经计算可知,直接转换法和面映射方法求解得到的排气歧管的位移分布与图12所示的分布相似;采用体映射方法得到的位移量最大,面映射方法次之,直接转换法最小,最大误差为7%左右。由此可以证明求解热应力问题时,以上3种方法都是可行的,对于后两种方式,只要插值函数选择合适就能够保证精度。4.2.2d不同单元的热应力变化图13和图14分别为体映射方法下,采用45号和92号单元得到排气歧管的VonMises应力分布。经过计算可知,92号和95号单元的VonMises应力分布近似;而45号单元的应力分布与前两者有较大的差别,且最大应力点数值比前两者约低35%。这是因为ANSYS中的45号4面体单元是退化1次单元,是一种典型的3维等应变单元,这种单元的内部应变保持不变,应力随之也保持不变。故在用ANSYS计算热应力或者机械应力时须慎用45号退化单元。表1列出了体映射法下3种不同单元计算所需的CPU时间和内存。采用2阶的92号和95号单元的CPU时间约是1阶的45号单元的9倍多。采用20节点的95单元要比10节点的92单元花费更多的时间和内存。所以在硬件条件允许的情况下,采用92单元既能节省时间又保证计算的精度;在受到硬件条件限制的情况下,应该采用1阶4节点的72单元。但由于72号单元是线性单元,因此比92号2次单元精度低。由于在高温下排气歧管的材料屈服极限比常温下降很多,从图14中可以看出,计算排气歧管管身部位的热应力约为98MPa,是热应力危险区域。通过进一步的分析,如热疲劳分析、模态分析等就能找出排气歧管产生断裂的根本原因,进而提出改进措施,解决排气歧管的实际问题。5分析单元类型(1)以内燃机热流体-热应力分析为例,重点阐述应用