切向气流下激光烧蚀全编织CFRP的仿真方法

切向气流下激光烧蚀全编织cfrp的仿真方法
技术领域
1.本发明属于激光加工和毁伤领域，尤其涉及一种切向气流下激光烧蚀全编织cfrp的仿真方法。


背景技术：

2.碳纤维增强树脂基复合材料(cfrp)是以碳纤维为增强材料、树脂为基体的层状复合材料。cfrp具有比强度高、比模量高、硬度高、耐高温、耐腐蚀、低膨胀性、抗疲劳性好、比重小和优异的力学性能等特点。全编织cfrp其由于拥有比单向cfrp和正交cfrp更高的强度和更优越的力学性能，因而被更广泛地应用于航空航天和军事武器制造，尤其在多用途无人机领域。现有的导弹拦截成本高，精度有限；而高能激光武器具有持续追踪打击、光速伤害、成本低廉等优势，因此展开激光对高速气流下全编织cfrp的烧蚀特性研究至关重要。
3.激光烧蚀cfrp是极为复杂的物化过程，其包括树脂的热解、热解产物的反应和流动、内部压力导致的不可逆膨胀和力学剥蚀现象、残碳和碳纤维的氧化以及氧化导致的材料强度的变化、热解羽烟颗粒对激光的屏蔽效应和对氧气传质的阻碍；若材料表面存在高速流动的气体，其烧蚀过程还包括外部的力学剥蚀效应，氧气传质的动态演化、材料表面对流换热导致的温度演化等；当温度达到碳的气化温度时，材料表面会发生升华。因此，造成切向气流下激光烧蚀cfrp材料的毁伤过程可描述为：热解、氧化、内部力学剥蚀、外部力学剥蚀、升华。激光作用cfrp材料过程非常剧烈，通过实验研究想量化烧蚀过程中的各类贡献是非常困难的，进而数值仿真已成为越来越重要的分析手段。
4.现有的数值仿真主要针对低功率下激光切割cfrp研究、低流速下的氧化烧蚀研究、高流速下的外部力学剥蚀研究等，鲜有充分考虑各个物化阶段性贡献分析的研究，对切向气流下激光烧蚀全编织cfrp的研究更为有限。全编织材料的特殊结构在烧蚀过程中表现得较单向和正交结构更为复杂，热物性参数更难以定量分析。在研究氧化烧蚀时和外部力学剥蚀过程中，很多研究直接将流体域氧气浓度和边界流速作为恒定常数计算，未准确考虑随烧蚀坑演化而产生得变量分布。目前而言，针对cfrp材料的内部力学剥蚀单纯考虑温度演化，并未结合实际去考虑内部力学剥蚀的阶段性特点，即随着热解的完全，内部力学剥蚀速率逐渐趋于零。
5.因此，必须发展切向气流下激光烧蚀全编织cfrp仿真方法，定性定量地分析热解、氧化、升华、外部力学剥蚀、内部力学剥蚀对烧蚀过程中的阶段性贡献，从而能够精确化激光对高速靶目标的毁伤机理过程，实现高效激光打击策略的制定。


技术实现要素：

6.本发明的目的在于提供一种切向气流下激光烧蚀全编织cfrp的仿真方法。
7.实现本发明目的的技术解决方案为：一种切向气流下激光烧蚀全编织cfrp仿真方法，包括如下步骤：
8.第一步：采用体积平均的方法对材料的密度、恒压热容、热导率进行均一化处理；
9.第二步：建立流体域和多孔介质域基本控制方程，以动量守恒定理、质量守恒定理和能量守恒定理为基础，并采用移动网格方法追踪界面退移和烧蚀过程；
10.第三步：构造致密结构物理模型，追踪烧蚀坑内材料表面的流体速度变化，通过边界相似性原理映射到材料表面，并替代原方程的流速项进行修正；
11.第四步：以致密结构为理论依据推导并设置材料各向异性的热导率；
12.第五步：考虑热解气体的压力变化做功和粘性耗散，并添加到多孔介质能量守恒方程进行修正；
13.第六步：追踪烧蚀坑内材料表面的流体的温度和压力的演化，并结合理想气体方程计算氧气浓度的分布，并替代原方程的恒定浓度项进行修正；
14.第七步：追踪材料表面的达西速度场演化，对内部力学剥蚀方程进行修正；
15.第八步：设定边界条件和初始值；
16.第九步：结合上述推导和计算的速率方程，设置网格移动条件和移动速率，通过有限元迭代计算，分析烧蚀过程中各个剥蚀贡献的演化。
17.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其所述处理器执行所述计算机程序时实现上述切向气流下激光烧蚀全编织cfrp仿真方法的步骤。
18.与现有技术相比，本发明的显著优点为：本发明提出了高速切向气流下和高功率激光密度作用下编织cfrp材料的致密结构外部力学剥蚀模型，在此模型的基础上推导了全编织cfrp材料的各向异性热导率。同时考虑了烧蚀坑涡流演化对氧气传质和外部力学剥蚀的影响，并结合材料表面达西速度场的演化对内部力学剥蚀效应进行修正，考虑热解气体压力变化做功和粘性耗散，进一步对材料内部能量守恒方程进行修正。该仿真方法通过建立流固耦合模型，同时分析了碳纤维增强环氧树脂复合材料在激光烧蚀过程中的树脂热解反应、碳纤维和残碳的氧化和升华剥蚀、外部切向气流导致的力学剥蚀以及内部热解气体产生的力学剥蚀问题，解决了碳纤维增强环氧树脂复合材料在激光烧蚀不同阶段中不同剥蚀的贡献分析问题。该数值仿真可以重现一定流速范围内的任意形式激光对cfrp烧蚀过程，为高流速切向气流下高功率激光毁伤复合材料提供理论和模拟依据。
附图说明
19.图1是切向气流下激光烧蚀全编织cfrp材料的仿真方法的二维建模示意图。
20.图2是致密结构物理模型示意图。
21.图3是烧蚀坑内不规则涡流温度和速度分布示意图。
22.图4是材料中心点的达西速度演化曲线。
23.图5是网格剖分示意图。
24.图6是不同剥蚀的贡献演化曲线。
具体实施方式
25.为使本领域技术人员更加清楚和明确本发明的技术方案，下面结合实施例及附图对本发明作近一步详细的描述，但本发明的实施方式不限于此。
26.本发明提供了一种切向气流下激光烧蚀全编织cfrp材料的仿真方法，建立数值模
型的过程中，引入高速切向气流下和高功率激光密度作用下编织cfrp材料的致密结构外部力学剥蚀模型，在此模型的基础上推导了全编织cfrp材料的各向异性热导率。充分考虑了烧蚀坑涡流演化对氧气传质和外部力学剥蚀的影响，并结合材料表面达西速度场的演化对内部力学剥蚀效应进行修正，考虑热解气体压力变化做功和粘性耗散，进一步对材料内部能量守恒方程进行修正，以此提高流固耦合模型的计算精度。
27.本发明的二维数值模型的示意图如图1所示。在保证数值计算结果可靠性的基础上，对模型做出以下假设和简化：
28.(1)假设热解、氧化、升华均不会影响流体域气体组分；
29.(2)忽略热解气体之间以及与材料之间复杂的化学反应；
30.(3)忽略材料的热膨胀和收缩效应；
31.(4)忽略材料不规则烧蚀拔出导致的流场重新分布；
32.该方法具体步骤如下：
33.第一步：采用体积平均的方法对材料的密度、恒压热容、热导率进行均一化处理。
[0034]vi
＝mi/ρii＝f,p,ch,g
[0035][0036][0037][0038][0039]
上式中各个组份的密度是假设恒定的，f,p,ch,g依次对应碳纤维、树脂、残碳、热解气体，表示各组份的体积分数，ρi表示各组份的恒定密度，c
pi
表示各组份额恒压热容，ξ为热解气体质量分数。
[0040]
第二步：建立控制方程，以动量守恒定理、质量守恒定理和能量守恒定理为基础，并采用移动网格方法追踪界面退移和烧蚀过程。
[0041]
流体域计算采用可压缩流体—守恒形式纳维斯托克斯方程：
[0042][0043][0044][0045]
[0046][0047][0048]
上式中，ρ
air
为空气密度，为空气速度矢量，p
air
为空气压强，μ
air
为空气运动粘度，f为力源项，c
p,air
为空气恒压热容，λ
air
为空气热导率，i为单位矩阵张量，q
vd
为粘性耗散项，q
p
为压力功项，为哈密顿算子。
[0049]
多孔介质域的基本控制方程结合了修正的傅里叶方程、单步阿伦纽斯方程、单相达西定律公式、康尼尔-卡曼渗透率模型理论，具体描述为：
[0050][0051][0052][0053][0054][0055]
上式中ρ
eff
为材料整体等效密度，c
p,eff
为材料整体等效比热容，为哈密顿算子，l
p
为热解潜热，ρg、ρ
p
分别为热解气体和树脂的密度，分别为孔隙率和树脂体积分数，为热解气体达西速度矢量，为热解反应指前因子，e
p
为热解反应活化能，为反应末端树脂体积分数，为初始树脂体积分数，κ为渗透率，pg为热解气体导致的达西压力，μg为热解气体动力粘度，d
particle
为多孔介质的粒径尺寸，r为普适气体常量，t为基体温度，t为计算时间，ξ为热解气体质量分数，n为反应级数。
[0056]
第三步：构造外部力学剥蚀致密结构物理模型，如图2所示，并追踪烧蚀坑内材料表面的流体速度变化，通过边界相似性原理映射到材料表面，并替代原方程的流速项进行修。
[0057]
在高功率激光辐照和高速切向气流作用下，碳纤维的力学剥蚀更为剧烈，一方面高温迅速降低了碳纤维的强度，另一方面树脂迅速热解较低功率下的缓慢热解更能保持材料的致密结构，因此基于以上的事实，本发明构造了外部力学剥蚀致密结构物理模型。
[0058]
流固耦合模型中，材料表面存在速度流体域，因此理论上紧贴材料表面的流速应是0。但是事实并非如此，材料表面不规则的烧蚀形貌会破坏理想的速度流体域，因此为表征切向气流对材料的力学剥蚀效应，通过追踪材料表面上方一定距离的流体速度，通过边界相似映射原理可以在不规则的烧蚀坑中准确地追踪变化的速度场，依次来修正力学剥蚀速率，进而更为准确地计算力学剥蚀的贡献。
[0059][0060][0061][0062][0063][0064][0065][0066][0067][0068][0069]
上式中v
out
为外部力学等效剥蚀速率，d
p
为树脂基体外部力学剥蚀等速率，du为外部径向碳纤维力学剥蚀速率，df为轴向碳纤维力学剥蚀速率，δ为纤维层与树脂层厚度的比值，ρf、ρ
p
为碳纤维和树脂的恒定密度，为碳纤维和树脂气化反应的指前因子，e
af
、e
ap
为碳纤维和树脂气化反应的活化能，λf、λ
p
为碳纤维和树脂的热导率，cf、c
p
为碳纤维树脂的恒压热容，分别为碳纤维轴向和径向的初始强度以及树脂的基体的初始强度，σf、σf⊥
、σ
pt
分别为轴向、径向和树脂基体的实时强度，分别为碳纤维和树脂基体强度变化因子，为修正强度因子的温度相关变量，分别为碳纤维、树脂、残碳的实时体积分数，为树脂的初始体积分数，tw为材料表面实时温度，t0为材料表面初始温度，t为激光作用时间，r为普适气体常量，p
∑
为材料表面压头。
[0070]
追踪烧蚀坑内材料表面的流体速度变化，如图1所示，并替代原方程的流速项进行修正，其中材料表面压头用材料表面流速修正描述为：
[0071]
p
∑
＝ρ
air
(u
airsurface
(t)
·
cos(θ))2[0072]
其中ρ
air
为空气密度，u
airsurface
(t)为随烧蚀坑演化的材料表面气流速度大小，θ为速度法向与水平面的夹角。
[0073]
第四步：以致密结构为理论依据推导并设置材料各向异性的热导率，描述为：
[0074][0075][0076][0077][0078][0079][0080][0081][0082]
上式中λ
solid
为固体等效热导率，λ1,λ2为材料表面方向热导率，λ3为材料厚度方向热导率；λ
ch
、λ
p
、λf分别为残碳、树脂、碳纤维的恒定热导率；分别为碳纤维、残碳、树脂的实时体积分数，为碳纤维的初始体积分数，a
*
、n1、n2、b
*
、sb、s
p
为中间变量，η为残碳热导率和基体热导率的比值。
[0083]
第五步：考虑热解气体的压力变化做功和粘性耗散，并添加到多孔介质能量守恒方程进行修正。
[0084]
热解气体的产生会导致材料内部压力的变化，与此同时由于内部与边界的压力梯度，热解气体会不断地往边界溢出，在热解气体逃逸的过程中压力在不断变化，并且热解气体与孔隙之间存在粘性，因此考虑压力变化做功和粘性耗散能更为精确地计算材料内部温度的演化规律。
[0085][0086][0087][0088]
[0089]
上式中ρ
eff
为描述多孔介质的等效密度，c
p,eff
为等效恒压热容，λ
eff
为等效热导矩阵，ρ
p
为树脂恒定密度，为树脂实时体积分数，l
p
为树脂热解潜热，c
pg
为热解气体恒压热容，ρg为热解气体密度，为热解气体达西速度矢量，pg为热解气体压力，q
vd
′
为热解气体流动产生的粘性耗散，q
p
′
为内部热解气体压力变化做功项，t为材料温度，t表示计算时间，为哈密顿算子。
[0090]
第六步：追踪烧蚀坑内近材料表面的流体温度和压力演化，并结合理想气体方程计算氧气浓度的分布，进而替代原方程的恒定浓度项进行修正。
[0091]
在流固耦合模型中，随着烧蚀坑的加深，烧蚀坑内会形成从下风去流往上风区的热流，并且随着流入速度的增大，这种涡旋热流愈发的不稳定，且往无序地方向发展，如图3所示。因此在烧蚀坑内会产生温度剧烈扰动的气流区域，气流温度大幅度地无规则演化将会导致烧蚀坑内氧气浓度分布产生差异，结合烧蚀坑内变化发展的压力，基于理想气体方程，可以计算处烧蚀坑内氧气浓度的演化分布，并以此修正氧化烧蚀速率方程，使得计算结果更为精确。
[0092][0093][0094]
上式中v
oxi
为氧化反应速率，a
oxi
为氧化反应指前因子，e
oxi
为氧化反应活化能，mi为参与氧化反应的表面材料平均摩尔质量，为参与氧化反应的等效表面密度，p
air,surface
为近烧蚀坑材料表面氧气分压，t
air,surface
为近烧蚀坑材料表面气体温度，t
surface
为材料表面温度，ζ为修正系数，r为普适气体常量。
[0095]
第七步：追踪材料表面的达西速度场演化，根据材料表面达西场的分布特征，对内部力学剥蚀方程进行修正。
[0096]
在激光烧蚀材料过程中，树脂率先在较低的温度发生热解反应，由于温升速率极大，热解气体短时间大量聚集，导致材料近表面内部压力迅速增大。在此阶段中，依据实验观察到的结果，发现内部力学剥蚀效应最为剧烈，而随着树脂热解的完全，没有更多的热解气体聚集并从表面逃逸，因此内部力学剥蚀效应逐渐降低并趋于稳定，如图4所示。基于以上事实，本发明通过追踪材料表面达西速度场的演化并以此修正内部力学剥蚀速率，表述为:
[0097][0098][0099]
上式中v
int
为内部力学剥蚀速率，为树脂热解反应的指前因子，λ
p
为树脂热导
率，ρ
p
为树脂密度，c
p
为树脂恒压热容，k为渗透率，l0为孔隙尺寸，σ
pt
为树脂强度，e
p
为树脂热解反应活化能，r为普适气体常数，t
surface
为材料表面温度，σ
compsite,t
为复合材料强度，φ为修正项，χ为修正系数，u为材料表面达西速，u
max
为材料表面达西速度峰值常数，ω为常数。
[0100]
第八步：依据赫兹-朗缪尔气化方程、克劳修斯-卡佩龙饱和蒸汽压方程计算材料表面的气化速率。
[0101][0102][0103][0104]
δh
sub
＝mcl
sub
[0105]
上式中，v
sub
为升华速率，m
sub
为碳耗散的质量流密度，为材料表面参与气化的等效密度，p
th
为饱和蒸汽压，β为逆扩散系数，δh
sub
为碳的气化焓，mc为碳的摩尔质量，l
sub
为碳的气化潜热，t
sub
为碳的气化温度，t
surface
为材料表面温度，p0为标准大气压，r为普适气体常量。
[0106]
第九步：设定边界条件和初始条件。
[0107]
流固边界上的能量守恒：
[0108][0109]
上式中λ
eff
为材料整体等效热导率，为温度梯度，α为材料对激光的吸收系数，i
lase
为功率密度分布，ε为辐射系数，σ为斯特藩-玻尔兹曼常量，t
surface
为材料表面温度，te为环境温度，为材料表面参与氧化和升华的等效密度，l
oxi
为氧化潜热，l
sub
为升华潜热，v
oxi
为氧化速率，v
sub
为升华速率。
[0110]
流固传热边界条件和初始条件：
[0111]
流体域上方边界：
[0112]
流体域入口处：
[0113]
流体域出口处：
[0114]
在固体下放边界上：
[0115]
初始条件：t
priamry
＝273.15[kl
[0116]
上式中为垂直于边界的法向矢量，λ
air
为空气热导率，为温度梯度，ρ
air
为空气密度，为入口焓差，h
al
为自然对流换热系数，t
enviroment
为环境温度，t为材料温
度。
[0117]
流体流动的边界条件和初始条件：
[0118]
流体域入口处设置法向流动：u
in
＝u
air
[0119]
流体域出口设置静压：p
out
＝1[atm]
[0120]
流体域上方边界设置对称：
[0121]
流体域下方延申边界设置滑移条件(具体如图所示)：u
air
＝u
in
[0122]
固体下方边界设置无滑移条件：u
air
＝0
[0123]
上式中为垂直于边界的法向矢量，u
air
为来流速度，p
out
为出口压力。
[0124]
达西定律边界条件和初始条件：
[0125]
多孔介质边界设置出口压力：p
boundary
＝1[atm]
[0126]
多孔介质内部初条件：p
priamry
＝1[atm]
[0127]
上式中p
boundary
为边界压力，p
priamry
为多孔介质内部初始压力。
[0128]
表征树脂热解的域微分方程初始条件：
[0129]
表征碳纤维体烧蚀的域微分方程初始条件：
[0130]
上式中为树脂初始体积分数，为碳纤维初始体积分数。
[0131]
第十步：对网格实行四边形划分，如图5所示，并设置网格移动条件和速率，通过有限元迭代计算，分析烧蚀过程中各个剥蚀贡献的演化，如图6所示。
[0132]vtotal
＝v
out
v
int
v
oxi
v
sub
[0133]
图6为材料中心点的烧蚀速率演化曲线，从图中可以看出在最初的烧蚀阶段，内部力学剥蚀和外部力学剥蚀速率增长趋势明显，这吻合热解气体导致的材料表面膨胀的物理事实，随着热解气体的完全，内部力学剥蚀在2.5s后基本对剥蚀不产生明显贡献，而伴随着凹槽深度的增加，近壁面的流体速度也逐渐降低，在2.3s后基本对剥蚀不产生明显贡献。随着温度的上升，材料表面的氧化速率逐渐上升至稳定值0.2mm/s左右，而在案例模型中由于温度未达到碳的气化温度，因此在图中未显示明显的气化剥蚀速率。在整个烧蚀过程中，从0至1.5s，整体烧蚀速率迅速上升，在1.5s左右达到峰值0.85mm/s，此后至2.3左右降低，并基本维持在0.23mm/s左右，材料也因此进入稳态烧蚀阶段。
[0134]
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。