基于修正任意拉格朗日欧拉算法的结构体入水仿真

曹力飞,李炳彦,刘 艳

(1.西安机电信息技术研究所,陕西 西安 710065;2.中国航发成都发动机有限公司,四川 成都 610503)

0 引言

入水冲击过程是一种典型的流固耦合问题,针对流固耦合的研究从20世纪30年代就已经开始,1929年Von Karman对楔形体和平头体的入水载荷进行初步计算,之后Wagner对Von Karman的算法进行改进,用于完成分析水面抬升现象,提出了销斜面升角理论。发展至今对于流固耦合的力学分析计算主要有两种方案:CFD仿真方法和ALE流固耦合方法。CFD仿真方案主要是通过动网格方法[1-2],实现作用过程中流体与固体的实时耦合计算[3],但该方案存在以下问题:1) 固体运动幅度太大会导致计算的域过大,进而致使网格数量太大,大幅增加计算周期与计算压力;2) 流体域的网格划分需要对固体部分进行掏空与绑定,对于复杂固体的网格划分难度较大[4]。例如:机翼在流场中的摆动、桥梁在风场中的震动等。ALE流固耦合方法的特点是流体与固体单独进行网格划分,计算过程中网格流体网格不发生变化,只在耦合时在耦合界面进行物理量的交换[5]。例如:流体在箱体中的运动、高超声速入水、流水冲击固体等案例[7]。

针对结构体入水的研究主要有CFD仿真方法和ALE流固耦合方法:任意拉格朗日-欧拉(ALE)方法分析结构体对于流体的流动特性,该方法可以分析高速入水与冲击过程,可以确定结构体入水的整个过程的力学特性,但是计算量过大。CFD算法关于流体对于结构体的粘性摩擦分布特性的研究,使用范围主要集中在连续流体领域和低速流体领域,对于流-固高速冲击问题无法计算,并且对于计算的网格质量要求较高,否则容易产生发散。

随着工程系统的复杂程度越来越高,网格的数量越来越大。一方面,更多的网格数量可以得到更加逼近真实的仿真环境;另一方面,更多的网格数量又势必会降低仿真效率。为了提高计算效率,对仿真服务器和计算机性能要求越来越高,影响产品研制周期,通过网格模型优化的方法已经越来越无法满足工程计算的需求。对采用定步长ALE算法的大变形、大数据量交换的结构体高速入水耦合计算影响更加明显。针对上述的不足,本文提出基于冲击修正的任意拉格朗日-欧拉(C-ALE)算法的结构体入水计算方法。

1 任意拉格朗日欧拉(ALE)算法

ALE(arbitrary Lagrange-Euler)连续介质力学有两种经典的运动描述方法,分别是Lagrange描述和Euler描述。Lagrange描述:网格节点固定在固定物质点并随之运动,因此在描述运动边界或者运动界面时非常方便;但当物质发生大变形时常常产生网格纠缠,轻则影响单元近似精度,重则使坐标中的Jocobian行列式的值近似于零或负数,从而使计算终止,或者引起局部严重误差。Euler描述:网格节点固定在空间始终不动,因此在描述大变形时,不会产生纠缠问题;但是,存在网格和物质的相对处理对流效应更加困难和无法精确确定运动边界的问题。为了克服Lagrange描述和Euler描述各自存在的缺点,引入相对坐标任意移动的方法即ALE方法。

结构体的入水过程是一个非线性过程,流体域包含了气液多相流流动,ALE算法主要包含质量、动量和能量守恒方程。

质量守恒方程:

动量守恒方程:

式中,X为ALE坐标;x为空间坐标;ci为ALE描述的对流速度,ci=vi-wi;vi为流体指点的物质速度;wi为网格速度;ρw为流体速度;bi为流体体积;σij为应力涨量。

能量守恒方程:

σij=-psδij+μd(vi,j+vj,i),

(3)

式(3)中,ps为水的静压;μd为动力粘性系数;δij为克罗内常数。

2 C-ALE算法的结构体入水计算方法

本文基于C-ALE算法实现结构体入水计算,其原理是在任意拉格朗日-欧拉(ALE)方法的基础上引入贯穿量C。求解过程如图1所示,采用交替计算的方法求解流固耦合问题,将流体和固体划分为两个独立的求解域,在耦合界面进行物理参数交换求解之后,分别对流体域和固体域进行单独计算。既保证流体域和固体域的计算可靠性,又保证耦合区间的物理量传递。利用C-ALE方法单独求解流体域,耦合界面上的结点力传递给结构,作为固体域的力边界条件,单独求解动力学方程。当流体域和固体域同时达到精度要求时结束迭代,根据耦合界面的运动情况更新单元网格进入下一时间步的求解。

图1 C-ALE流固耦合计算模型Fig.1 C-ALE fluid-solid calculation model

具体计算过程中,检查每个计算步长的网格节点是否有物质穿透,没有则对该网格不做处理。如果发生物质穿透,则通过界面力F参与计算。计算过程中界面力与贯穿量C成正比:

F=kt·C,

(4)

式(4)中,kt物质模型刚度系数,C贯穿量。

(5)

式(5)中,ΔEm为动能变化量,Fd为冲击阻力,s水流冲击距离。

通过贯穿量C的引入,加快了在耦合界面的能量传递参数的计算速度,尤其在高速流体冲击过程中,对于结构体的冲击应力的计算。但是,对于低速入水的结构体入水过程影响不大。

2.1 算例模型

入水模型由入水结构体、弹体、空气和水组成,如图2所示。其中弹体和入水结构体作为固体域参与计算,水和空气作为流体域参与计算。耦合界面包含水-空气流体耦合、弹体与入水结构体固体耦合、入水结构体和水组成流固耦合。

图2 结构体入水仿真模型Fig.2 Structure water entry simulation model

流体域主要采用欧拉网格进行描述,而固体域则采用拉格朗日网格描述。材料网格之间可以进行材料交换,而不需要产生网格运行与生成[8]。对固体和流体域在建立耦合边界条件之后,分别进行网格划分。网格划分如图3、图4所示。

图3 流体域网格划分Fig.3 Fluid domain meshing

通过粘性应力本构关系,建立水的Gruneisen状态方程材料模型,模型参数如表1所示。

(6)

式(6)中,μ为比体积;p为压力;E为单位体积内能;C,s1,s2,s3为材料常数;γ0为水的密度;a为对γ0的一阶体积修正。

图4 固体域网格划分Fig.4 Solid domain meshing

表1 水和空气的Gruneisen状态方程材料模型参数Tab.1 Gruneisen equation of state material model parameters of water and air

入水结构体选用铝合金2A12,材料模型选取Johnson-Cook模型,模型参数如表2所示。Johnson-Cook模型本质是将应变、应变率和温度三个变量进行分离,利用乘积关系处理三者对动态屈服应力的影响,具有形式简单、各项物理意义明确的优点。

表2 2A12材料Johnson-Cook模型参数Tab.2 Johnson-Cook model parameters of 2A12

弹体部分处理为质量单元,在不影响计算准确性的前提下以增加计算效率。

2.2 流体域计算准确性

通过空泡理论对仿真的准确性进行验证。入水速度取350 m/s。对于结构体入水难以准确测量,只能通过空泡的外形进行对比验证。基于空泡理论可以相对准确的计算入水过程的空泡外形[9]。

式(7)中,db为弹体直径;Cd为入水阻力系数,该系数与入水的速度有关。

根据采用C-ALE算法进行对比可以得到仿真空泡轮廓,弹丸入水过程体积分数变化如图5所示。

图5 弹丸入水过程体积分数变化Fig.5 Variation of volume fraction of projectile during water

C-ALE算法与空泡理论计算对比曲线如图6所示。可知基于C-ALE算法入水计算曲线与空泡理论曲线特征相近,空泡深度空泡和开口尺寸基本一致,因此,基于C-ALE算法的结构体入水计算结果准确,算法合理可行。

图6 理论空泡曲线与计算空泡曲线对比Fig.6 Comparison of theoretical and calculated cavitation curves

2.3 固体域计算准确性

以350 m/s的速度进行基于冲击修正的任意拉格朗日-欧拉(C-ALE)方法的结构体入水数值分析。入水过程结构体速度和加速度变化曲线如图7、图8所示。

仿真结果显示:弹丸入水过程结构体速度稳定下降,原因是由于空泡的产生只有入水结构体的部分接触水面并产生阻力;至6 ms后下降速度加快,分析原因为锥形段部分开始浸水导致摩擦阻力增加。

图7 入水过程结构体速度变化Fig.7 Variation of structure velocity during water entry

入水结构体的变形主要发生在入水结构体接触水面的瞬间,此时过载虽然没有形成最大过载,但是此时接触面积较小。

图8 入水过程弹丸加速度变化Fig.8 Variation of projectile acceleration during water entry

在入水结构体中形成超过极限屈服载荷的压力,入水结构体产生形变。最大形变为3.520 6 mm,时间发生在0.122 ms,与入水过载峰值时间基本保持一致。随后形变发生一定回弹,并保持稳定,最终形变为3.231 1 mm。入水结构体的变形阶段主要发生在入水结构体与水接触的初期。当入水结构体侵入水面之后,随着锥形弹丸与水的接触面积的增加,变形不再变化。

根据图9入水前后结构体形变可知:仿真结果显示形变主要发生在结构体的头部,在锥形段部分没有发生相关形变。

外场试验时,测试弹丸以350 m/s的速度进行入水冲击试验,试验前后入水结构体变形情况如图10所示。

试验结果表明:1) 结构体入水的形变主要发生在结构体的头部,这与仿真结果保持一致;2) 外场试验实测产品产生3.10 mm的变形,仿真计算入水结构体产生3.231 1 mm变形,仿真计算与外场试验结果基本一致。

图9 入水后结构体形变云图Fig.9 Deformation of structure after water entry

图10 入水前后结构体变形(试验结果)Fig.10 Variation of projectile acceleration during water entry with structure

验证了冲击修正的任意拉格朗日-欧拉(C-ALE)方法对锥形壳体入水冲击计算准确,可以用于指导结构体高速入水的仿真。

3 适应性与数据分析

C-ALE算法与ALE算法采用同样的网格模型,因此可以通过ALE算法进一步验证C-ALE算法对结构体入水过程仿真的适用性。

3.1 C-ALE算法与结构体形变收敛特性

由图11和表3可知,采用同一网格模型,在相同输入条件下,通过ALE算法进行计算。ALE算法计算收敛时间约为8 976 s,收敛后,入水结构体的结构形变为3.76 mm。C-ALE算法计算收敛时间约为3 342 s,入水结构体变形3.230 9 mm。相对误差从21.2%减小到4.2%,绝对误差由0.66 mm减小到0.13 mm,收敛时间8 976 s减小到3 342 s。可以得出:针对结构体入水计算过程,相同输入条件下,C-ALE算法比ALE算法具有更小的计算误差和更快的收敛速度。

图11 不同计算方法的收敛过程Fig.11 Convergence process of different calculation methods

表3 不同计算方法的结构体入水计算结果Tab.3 Calculation results of water inflow of structure with different calculation

3.2 入水深度与计算时间的关系

为了进一步验证C-ALE算法对于结构体入水深度的计算时间,在350 m/s入水速度条件下,采用同样的网格和材料模型,对比ALE算法对结构体入水50 mm深度的仿真过程进行监测,获得C-ALE和ALE算法入水深度与计算时间关系如图12所示。

图12 C-ALE和ALE算法入水深度与计算时间关系Fig.12 C-ALE和ALE algorithm relationship between water depth and calculation time

根据图12可知,在高速入水过程中,随着入水深度的增加,C-ALE算法相比ALE算仿真时间大幅减少。

3.3 入水速度与计算时间的关系

为确定不同入水速度对于算法仿真速度的影响,入水速度分别取50、100、150、200、250、300、350 m/s,采用同样的网格和材料模型,对比ALE算法对结构体入水50 mm深度的仿真过程进行监测,获得不同入水速度C-ALE和ALE下入水深度与计算时间如图13所示。

图13 不同入水速度C-ALE和ALE入水深度与计算时间Fig.13 Different entry velocity C-ALE and ALE entry depth and calculation time

由图13可知:1) 在相同网格和材料模型下,不同入水速度的C-ALE仿真过程均比ALE算法的仿真时间少;2) 在相同网格和材料模型下,当入水速度大于150 m/s,C-ALE算法与ALE算法的仿真时间差距增大。50 mm深度不同入水速度下结构体形变与计算时间如表4所示。

表4 50 mm深度不同入水速度下结构体形变与计算时间Tab.4 Structure deformation and calculation time with different water entry velocity at 50mm depth

根据表4可知,入水速度在150 m/s时候,结构体开始产生形变,随着入水速度的增加,结构体的形变增加,C-ALE算法中引入贯穿量C对于仿真的影响增大,对比ALE算法,进一步缩短了仿真时间。

3.4 网格数量与收敛性

进一步研究不同网格数量对C-ALE算法在入水仿真过程的收敛特性的影响,分别对如下网格特性进行仿真计算。A网格特性:流体域网格数量675,固体域网格数量373;B网格特性:流体域网格数量2 325,固体域网格数量1 562;C网格特性:流体域网格数量6 566,固体域网格数量2 987。不同网格特性的收敛过程如图14所示。仿真条件为350 m/s入水冲击过程,实际形变为冲击试验形变。C-ALE算法下不同网格特性下结构体形变收敛时间如表5所示。

图14 不同网格特性的收敛过程Fig.14 Convergence process of different mesh methods

表5 不同网格特性下结构体形变与收敛时间Tab.5 Structure deformation and calculation time with different mesh methods

由图14和表5可知,不同流体域与固体域的网格数量仅影响仿真的收敛速度,对仿真的准确度没有影响。

4 结论

本文提出了基于冲击修正的任意拉格朗日-欧拉(C-ALE)算法的结构体入水计算方法,该方法引入贯穿量C对任意拉格朗日-欧拉方法进行修正,采用流体域与固体域定步长交替计算的方法进行求解。通过空泡理论和外场试验结果验证了C-ALE计算结果的准确性;同时,对比采用相同输入条件和网格模型的ALE算法可以得出,对于结构体入水过程的求解,C-ALE算法比ALE算法具有更快的收敛速度和更小的计算误差。在相同网格模型和材料模型下,不同的入水速度,随着结构体形变的增加,C-ALE算法可以进一步缩减仿真时间。可以得出如下结论:通过贯穿量C的引入,加快了耦合界面的参数传递效率,提高了仿真计算速度。