NUMERICAL SIMULATION OF MICRO-ABLATION BEHAVIOR FOR CARBON FIBER REINFORCED COMPOSITES1)
LiWei中图分类号:O346.1
文献标识码:A
通讯作者:
收稿日期:2018-10-24
网络出版日期:2019-05-18
版权声明:2019力学学报期刊社 所有
基金资助:
展开
摘要
关键词:
Abstract
Keywords:
-->0
PDF (22787KB)元数据多维度评价相关文章收藏文章
本文引用格式导出EndNoteRisBibtex收藏本文-->
引 言
高超声速飞行器日益受到广泛重视,而其所涉及的气动加热问题亦成为学术界研究的热点[1-3]. 碳纤维增强复合材料如碳/碳复合材料[4]、碳/酚醛复合材料[5]以及新一代的酚醛浸渍碳纤维复合材料[6]等,作为烧蚀型防热材料在再入式飞行器如飞船返回舱等热防护系统上得到了广泛的应用[7-9]. 自20世纪60年代以来,宏观烧蚀理论[10-15]逐渐建立起来并已经成功应用于诸多热防护系统设计. 宏观烧蚀理论认为材料表面是均匀的碳材料,并且材料表面各处以相同的反应速率烧蚀. 但实际上烧蚀复合材料的基体相和纤维相的烧蚀速率差异较大,这种差异使得材料在烧蚀过程中形成复杂的表面微细观形貌. 同时形成粗糙的烧蚀表面又会引起近壁面附近的气流及热流发生改变,进而引起表面烧蚀状态的改变. 已有研究[16-17]发现:即使宏观组分相同而微细观结构不同的碳纤维增强复合材料,其烧蚀性能也存在显著差异,这种差异不能通过宏观烧蚀模型来描述. 因此需要从材料微观角度研究碳材料的微观烧蚀机制,从而为烧蚀防热的碳纤维增强复合材料设计和制造奠定理论基础.已有****试图从微细观尺度上探究碳纤维增强复合材料的烧蚀问题. Georges等[18]通过建立表面传热、传质以及界面后退模型,从理论上分析了碳基材料表面粗糙度的产生及其特征,发现了不同尺度粗糙度产生的可能原因,从而解释再入体表面流态从层流向湍流转捩的原因. Aspa等[19]对固体火箭发动机喉道处的碳/碳复合材料细观烧蚀现象进行了研究,假设烧蚀表面附近气体输运以纯扩散方式为主,以有限体积法求解气体扩散,以体积分数法(VOF)捕捉气$\!$-$\!$-$\!$固界面,对纤维和基体在水蒸气和二氧化碳氛围下的细观烧蚀进行数值模拟,获得了细观形貌演化以及等效烧蚀速率. Vignoles等[20]在文献[19]的基础上,进一步研究了编织碳/碳复合材料在固体火箭发动机喷管中的细观烧蚀问题,研究发现:尽管流动起重要作用,但烧蚀形貌细观演化的产生首先是由材料结构引起的;并且随着温度的升高,材料的氧化深度逐渐减小,氧化反应集中于烧蚀表面附近,材料的几何结构对等效烧蚀速率的影响逐渐减弱;材料的等效反应速率不满足混合率,而是由最弱相决定的. Lachaud等[21-23]分别从实验、理论和数值模拟方面,对碳纤维复合材料微细观烧蚀行为研究,获得了温度达到898,K的碳纤维与碳基体的本征氧化速率,由精细化实验观测到碳纤维复合材料在898,K条件下的微细观形貌特征,并通过布朗运动算法模拟了微细观烧蚀形貌的产生过程;通过理论分析,得到了单根纤维在烧蚀过程中形貌的解析解. 国内亦有****对微观尺度研究了烧蚀问题,黄飞等[24]应用直接蒙特卡罗模拟(DSMC)对简化的碳/碳复合材料烧蚀形貌进行模拟,探讨了氧分子与材料表面碰撞的作用规律,分析了纤维烧蚀结构的演变特性. Deng等[25]对微观烧蚀条件下产生的"笋尖"状形貌建立侵蚀模型,由材料力学理论给出了单根纤维在微观烧蚀条件下的断裂准则. 综上所述,目前国内外****主要针对纤维垂直于烧蚀表面的单根碳纤维烧蚀过程进行了研究,并且取得了一些重要结果.
由于碳纤维增强复合材料制造工艺不同,其碳纤维方向和烧蚀表面之间存在倾斜角度,探讨倾斜角度对微观烧蚀行为的影响有助于改进烧蚀复合材料的编织工艺. 另一方面,对于基体抗氧化改性的复合材料,纤维因为烧蚀速率大于基体而形成深坑,纤维快速烧蚀形成的深坑将显著影响材料表面的微观氧通道的生成和发展,因此有必要讨论倾斜角度的影响以指导抗氧化基体改性. 本文主要针对纤维存在倾斜角的碳纤维增强复合材料微观烧蚀问题进行研究. 首先通过材料烧蚀过程分析,建立碳材料微观烧蚀的简化模型,然后对模型进行求解并验证,最后分析了不同碳纤维倾斜角对材料微观烧蚀行为的影响规律.
1 微细观烧蚀模型
对于碳/碳复合材料来说,在烧蚀材料表面仅含有碳元素,表面碳材料的质量损失包含两部分因素:热化学烧蚀和机械剥蚀. 机械剥蚀是指材料表面在高温气流下发生的颗粒剥落或由热应力引起的片状崩落[26-27]. 对于常见的烧蚀型复合材料的使用条件,机械剥蚀并不严重,因此可忽略机械剥蚀.[15]. 热化学烧蚀指材料在高温下发生的氧化、氮化和升华反应. 对于大部分应用场合,碳纤维增强复合材料将主要面对氧化作用产生的烧蚀,因此本文只研究氧化引起的微观烧蚀现象.烧蚀表面是高速的外流场边界层,在边界层内发生复杂的传热、传质过程[14-15,26]. 但是已有研究指出在微观尺度下边界层的氧输运现象可以认为是由纯扩散现象主控[20],氧气将由边界层外缘通过纯扩散过程进入边界层内部和材料表面的碳元素反应. 在材料表面发生碳的氧化,主要包含以下几个反应
C + O$_2$ $\to$ CO$_2$ (1)
C + $\dfrac 12 $ O$_2$ $\to$ CO (2)
2CO $\rightleftharpoons$ CO$_2$ + C (3)
高温条件下的碳氧反应是复杂的非均相反应.
由于碳氧反应的速率受到扩散控制、杂质催化等现象的干扰[21],长
期以来很难获得碳氧反应的本征反应速率,且已有的碳氧反应测量值相差可达数千倍[28].
Lachaud通过精细化的实验设计[21],在898,K温度下测量了碳纤维复合材料的本征反应速率,并观测了微观烧蚀
形貌,发现微观尺度条件下碳纤维的烧蚀形貌呈"针尖状".
在898,K温度条件下,将只发生反应(1),这时控制方程[22-23]将如下所示
$S\left( {x,y,z,t} \right) = 0 (4)$
$\dfrac{\partial S}{\partial t} + {\pmb v}_e \cdot \nabla S = 0 (5)$
$ {\pmb v}_e = \varOmega _e J_e {\pmb n} (6)$
${\pmb n} = \nabla S / \left\| {\nabla S} \right\| (7)$
$J_e = k_e C (8)$
$\dfrac{\partial C}{\partial t} + \nabla \cdot \left( { - D\nabla C} \right) = 0 (9)$
其中,$S$是材料表面函数,在气相$\!$-$\!$-$\!$固相交界面有$S =0$, $C$指反应气体的摩尔浓度,$ {\pmb v}_e $ 是界面后退速度, $C_0 $是边界层外缘浓度,$k$是非均相反应速率,下标$e$代表不同相的碳元素(基体、纤维), $\varOmega _e $是固体相$e$的摩尔体积,${\pmb n}$ 是固体相界面的外法向.
Lachaud[23]采用布朗运动求解方程组(4)$\sim$(9),所得到的形貌与其在898,K 下的精细实验[21]观测结果十分接近,证实了这种方法的可行性. 并且Lachaud认为碳纤维微观烧蚀行为由两个无量纲参数$Sh$和$A$决定,其中$Sh = {R_{\rm f} k_i }/ D$ 体现氧化反应和气体扩散速率之比,$R_{\rm f} $ 是纤维半径,$k_i $ 是基体氧化反应速率,通常认为$Sh > 100$ 是扩散控制,而$Sh<0.01$ 是速率控制;$A = {k_i \varOmega _i }/ {k_{\rm f} \varOmega _{\rm f} }$ 体现基体相和纤维相在相同的氧摩尔浓度下反应体积之比,其中$\varOmega _i $ 是基体相的摩尔体积,$\varOmega _{\rm f} $是纤维相的摩尔体积, $k_{\rm f} $是碳纤维的氧化反应速率.
在实际的烧蚀条件下,氧化烧蚀受诸多因素(温度、氧分压、非连续流动等)控制,且材料表面温度远高于898,K时表面的碳氧反应会包含反应(1)$\sim$(3),而且在更高温度时的碳氧本征反应速率尚无可信数据,因此要想完全精确的模拟微观烧蚀过程是十分困难的. 但是我们依然可以用变化的$A$和${Sh}$体现出温度、氧分压以及反应速率的相对关系,给出各种烧蚀因素对烧蚀行为的影响趋势.
2 数值模拟及求解
2.1 数值模拟方法
在微观烧蚀模型中,固体界面在氧化作用下会发生烧蚀后退,且气相扩散和界面后退现象相互耦合. 一方面,材料表面氧化速率与表面附近的氧浓度有关,同时表面后退又影响浓度分布;另一方面,材料的结构在氧化过程中,会发生界面几何的融合、消失,这些都给该问题的数值模拟带来极大挑战. 本问题有以下难点需要解决:(1) 如何处理固体相的演化?
(2) 如何在不断演化的气$\!$-$\!$-$\!$固界面上计算氧化速率?
(3) 如何将互相耦合的界面运动和气体扩散相耦合?
上述问题类似于计算流体力学中的气$\!$-$\!$-$\!$液两相流,可以借鉴两相流的方法. 数值求解此类问题的方法通常分为两大类:基于拉格朗日描述的界面追踪法和基于欧拉坐标系的界面捕捉法[29]. 界面追踪法需要用网格追踪界面的移动,这类方法计算量大、且难以处理几何拓扑变化,因此本研究采用界面捕捉法.
VOF (volume-of-fluid)方法是解决两相流的常用方法,类似地引入固体体积分数法[19-20],定义固体体积分数$\varepsilon _{\rm s} = \dfrac{V_{{\rm solid}} }{V}$,结合有限体积法对扩散和烧蚀问题进行了求解,采用这种方法,得到
$ \dfrac{\partial C}{\partial t} + \nabla \cdot \left( J \right) = 0 \,, \hskip 6mm \ \hbox{in fluid} (10)$
$\dfrac{\partial \varepsilon _{\rm s} }{\partial t} + v_{\rm s} \nabla \cdot \left( { - J} \right) = 0\,, \ \hbox{in solid} (11)$
$1=\left\{\begin{matrix} - 1( 1 - \varepsilon _{\rm s} ) D\nabla C , \hbox{in fluid and solid} \\ - kC|| {\nabla \varepsilon _{\rm s} } ||^{ - 1}\nabla \varepsilon _{\rm s} , \hbox{fluid-solid interface} \\\end{matrix}\right.$
上述方程组,引入固体体积分数实现了气相和固相的处理.但是还需要解决的问题是识别气$\!$-$\!$-$\!$固界面,这就需要对计算过程中的固体域进行重构,并识别出各相间的接触面积. 本文采用PLIC (piecewise linear interface calculation)算法[30]进行几何重构,并识别界面,原理如图1所示.
显示原图|下载原图ZIP|生成PPT
图1计算各相间接触面积示意图
-->Fig. 1Illustration of phase contact area
-->
由于上述方程组中,摩尔浓度$C$和$\varepsilon _{\rm s}$ 相互耦合,在存在界面移动的条件下难以解耦.该问题中为了提高计算效率,且保持较高精度,我们采用显式算法求解固体相分布,对相分布进行光滑处理并采用PLIC算法[30]进行重构. 由重构所得的界面识别出$\partial \varOmega _{\rm f / f} $ 以及$\partial \varOmega _{\rm f / s} $,然后采用隐式算法求解浓度分布.从而采用松耦合算法进行解耦,耦合算法顺序如图2所示.
显示原图|下载原图ZIP|生成PPT
图2松耦合算法示意图
-->Fig. 2Illustration of loosely coupled algorithm
-->
采用自编程序对上述问题进行求解,程序由C语言与MATLAB语言混合编写.基于有限体积法计算,空间离散采用中心差分格式, 由于松耦合算法的时间精度最高仅为一阶[31-32],因此计算$C^{t + \Delta t}$ 时间采用一阶隐式格式推进.
2.2 松耦合算法验证
对于一维有界面后退的问题,不涉及复杂的几何重构,可利用本程序与商业软件进行对比.商业软件基于COMSOL多物理场软件, 采用有限元法、动网格法求解此问题. 为此构建一维烧蚀问题.如图3所示的计算域,在$t= 0$时刻$A$-$B$是气相域,$B$-$C$是固体域,$B$是氧化界面. 在一定的温度下$B$发生氧化,界面向右移动,如果气相域的浓度分布、界面位置吻合,则说明程序的松耦合算法精度足够. 其中,$L_{ AC} = 1.0\times 10^{ - 3}$,m,$L_{ AB} = 0.53\times 10^{ - 3}$,m,$v_{\rm s} = 6.67\times 10^{ - 2}$,m$^3$/mol,$k_{\rm c}= 2\times 10^{ - 2}$,m/s, $D = 2.0\times 10^{-4}$\,m$^2$/s,$C_{x = A} = 0.23$,mol/m$^3$.
显示原图|下载原图ZIP|生成PPT
图3松耦合算法验证示意图
-->Fig. 3Sketch map of the validation of loose coupling algorithm
-->
界面$B$在氧化条件下由初始位置运动到$C$位置,总共约1.7,s, 选取$t =0.5$\,s和$t =1.6 $,s两个时刻对比两种求解方式的结果.
由图4和表1可以看出,程序计算所得的气相分布、界面位置与商业软件吻合很好,说明程序的耦合算法及离散方法是可靠的.
Table 1
表1
表1界面位置
Table 1Interface position
新窗口打开
显示原图|下载原图ZIP|生成PPT
图4氧浓度分布
-->Fig. 4Oxygen concentration distribution
-->
2.3 三维烧蚀形貌的验证
Lachaud等[22]给出纤维单丝在微观烧蚀条件下的解析解,我们将根据解析解的形貌和计算所得的形貌进行对比验证.对基体包裹的纤维单丝在氧气扩散条件下的微观烧蚀行为进行模拟,计算模型如图5所示.因为纤维和基体之间存在界面薄层,而本文经过研究发现未考虑界面时所得的烧蚀形貌与文献中一致,为简化模型并减少计算量,在本文的计算模型均不考虑纤维和基体的界面.
显示原图|下载原图ZIP|生成PPT
图5计算域示意图
-->Fig. 5Schematic diagram of mode
-->
Lachaud[22]指出微观尺度的烧蚀形貌,与边界层厚度以及浓度$C_0 $ 无关,为模拟高空中的氧浓度本文取$C_0 =0.23,mol\/m^3$. 由于精确的碳氧反应速率数据难以获得,在计算中我们采用以下数据[22]:$D = 2\times 10^{ -4}$,m$^2/s, R_{\rm f}= 5,\mu$m, $\varOmega _i = \varOmega _{\rm f} = 6\times 10^{ - 6}$,m$^3$/mol, $A= 5$. 在具体计算时根据$Sh$ 和$A$推出所需的反应速率,例如$Sh = 0.01$ 时有$k_i = 0.4$\,m/s, $k_{\rm f} =0.08$\,m/s. 实际在烧蚀过程中,碳氧化的本征反应速率将随着温度的升高呈指数倍增大,因此在不同温度下的烧蚀行为将由$Sh$ 体现,在特定的$Sh$ 下认为材料处于等温状态.
由图6可见本研究所用程序计算结果与解析解吻合良好,在$A =5$条件下基体烧蚀速率大而纤维烧蚀速率低,使烧蚀过程中纤维高于基体平面. 在$Sh = 0.01$时稳定的烧蚀形态中纤维呈"笋尖状";随着$Sh$增大($Sh = 0.1$和$Sh = 1$),"笋尖"高度降低;当$Sh =5$时纤维形态呈"半球形";当$Sh = 10$时纤维呈"平台状"而略微高于基体表面;最后当$Sh =100$时纤维和基体后退速率几乎相同,纤维将不再高于基体平面.
显示原图|下载原图ZIP|生成PPT
图6数值模拟形貌与解析解形貌对照
-->Fig. 6Comparison of the morphologies from simulation results and analytical results
-->
在$A =5$条件下的材料稳定烧蚀形貌如图7所示. 在不同$Sh$下材料产生不同形貌的原因是氧化反应速率和氧扩散速率之间的竞争作用主导的.当$Sh$接近0.01时气体扩散速率远远高于氧化反应速率,这时氧气供应充足,在充分的氧化作用下靠近气相扩散边界处的纤维将比靠近底端的纤维更快的氧化,因此纤维直径将从上到下依次递减从而形成"笋尖状"形貌. 而当$Sh$接近100时,反应速率远远大于氧气扩散速率,这时在气相$\!$-$\!$-$\!$固相交界面的氧气消耗量远远高于供应量,在缺氧条件下纤维和基体的后退速率都将被相对缓慢的气体扩散速率制约,此时烧蚀进入所谓的"扩散控制"状态,纤维和基体将同步后退.
显示原图|下载原图ZIP|生成PPT
图 7单根纤维在不同Sh数下的烧蚀形貌($A=5$)
-->Fig. 7Ablation morphologies of single fiber at different Sh numbers ($A=5$)
-->
3 倾斜角度对微观烧蚀行为的影响
真实的碳纤维增强复合材料中,纤维与烧蚀表面会存在夹角,且纤维与基体的抗氧化性能不同,和相同摩尔量的氧气反应后产生的碳体积损失也不同. 根据前文2.1节定义的无量纲参数$A$,可以将纤维和基体的抗氧化能力分为两大类:当$A >1$时,纤维抗氧化性强于基体(以下简称"纤维强基体弱",如图8(a)所示[33]),可以观察到笋尖状纤维烧蚀形貌.当$A < 1$时,纤维抗氧化性弱于基体(以下简称"纤维强基体弱",如图8(b)所示),则可以观测到纤维快速烧蚀留下的深坑.显示原图|下载原图ZIP|生成PPT
图8实验观测的典型微观烧蚀形貌
-->Fig. 8Typical micro-ablation morphologies
-->
为研究碳纤维角度对材料微观烧蚀行为的影响,对不同的两种碳纤维和基体抗氧化性能关系($A =5$和$A=0.2$),在不同$Sh$条件 下的微观烧蚀行为进行模拟. 如图9所示,其中倾斜角度 $\theta$,指纤维取向与烧蚀表面的夹角,取0$^\circ$$\sim $ 90$^\circ$. 不同倾斜角度对应的纤维体积分数介于0.22$\sim$0.26间.
显示原图|下载原图ZIP|生成PPT
图9有倾斜角度的材料示意图
-->Fig. 9Schematic diagram of material with different fiber inclination angle
-->
由于高温条件下的碳氧反应十分复杂,因此缺乏精确的碳氧反应速率数据.但可以通过本模型定性地分析碳纤维复合材料在不同反应 机理下的烧蚀形貌特征,并寻求基本规律. 在参考文献中所提供数据的基础上,保持其他参数不变,通过改变 $k_i $, $k_{\rm f}$的参数获得不同的$A$和${Sh}$,分别模拟不同碳氧反应段的典型形貌规律.
3.1 纤维抗氧化性能比基体强的情况
首先针对碳纤维抗氧化性强而基体抗氧化性弱的情况进行模拟. 如图 10给出了$Sh = 0.01 $时不同倾斜角材料在3个典型时刻的形貌,可以看出在所有角度和时刻,纤维均高于基体表面,形成 笋尖状形貌.并且对于不同角度和时刻,纤维的顶点处于几乎相同的高度. 图 10是$Sh = 0.01 $时不同倾斜角材料的质量损失过程,可以看出不同的倾斜角失重速率有 差异,但是差异不大.由于计算模型在水平方向不可能无限大,造成不同倾斜角的模型本身存在纤维基体含量的差异,因此从模拟图中体现出的质量损失速率微小差异,并不是单纯由于材料结构的不同而引起的.Lachaud等[22]同样通过解析解分析认为倾斜角对材料的等效烧蚀速率没有影响. 因此可以得出结论,在$Sh=0.01$条件 下不同倾斜角对材料的烧蚀行为影响不显著,即使有影响也很微小.
另一方面,可以看出在大倾斜角度条件下存在部分纤维完全失去基体保护而"悬空"的现象,例如图 10中$\theta = 60^\circ$, $90^ \circ $的形貌,产生这种形貌是由于基体烧蚀速率大于纤维烧蚀速率,使纤维露出一定高度$h_{f}$,并且高度在$Sh\mbox{ = 0.01}$时与$A$有关. 在大倾斜角度时,高出基体$h_f$的纤维将完全失去基体保护而"悬空". 但是可以看出图10中$\theta = 45^ \circ$, $60^ \circ $和$75^ \circ$时的"悬空"是由于模型尺寸造成的,当模型足够大时则不会出现该现象. 而对于$\theta = 90^ \circ $则是完全是因为基体的后退造成的,在现实条件下这种"悬空"将使得纤维整体剥蚀.
显示原图|下载原图ZIP|生成PPT
图10不同时刻微观烧蚀形貌 ($A=5,Sh=0.01$)
-->Fig. 10Micro-ablation morphologies at different times ($A=5,Sh=0.01$)
-->
显示原图|下载原图ZIP|生成PPT
图11不同纤维倾斜角度的材料质量损失速率($A=5, Sh=0.01$)
-->Fig. 11Mass loss rate of materials at different fiber inclination angles ($A=5, Sh=0.01$)
-->
对于$Sh = 1 $条件下,各倾斜角的材料烧蚀行为与$Sh = 0.01$时类似(图12),只是纤维和基体之间烧蚀速率差异减小了,并且$\theta = 90^ \circ $时依然存在纤维"悬空"现象.结合形貌和失重速率(图13)依然可以得知$Sh = 1$时倾斜角对烧蚀行为没有显著影响.
显示原图|下载原图ZIP|生成PPT
图12不同时刻微观烧蚀形貌($A=5, Sh=1$)
-->Fig. 12Micro-ablation morphologies at different times ($A=5, Sh=1$)
-->
显示原图|下载原图ZIP|生成PPT
图13不同纤维倾斜角度的材料质量损失速率($A=5, Sh=1$)
-->Fig. 13Mass loss rate of materials at different fiber inclination angles ($A=5$, $Sh=1$)
-->
随着$Sh$进一步增大时,纤维和基体之间烧蚀差异逐渐消失,倾斜角对烧蚀行为将完全没影响,纤维和基体将以相同速率烧蚀,如图14所示.
显示原图|下载原图ZIP|生成PPT
图14不同时刻微观烧蚀形貌 ($A=5$, $Sh=100$)
-->Fig. 14Micro-ablation morphologies at different times ($A=5$, $Sh=100$)
-->
3.2 纤维抗氧化性能比基体弱的情况
对于纤维抗氧化性弱而基体强的情况,采用$A$=0.2进行碳纤维在不同倾斜角度下的微观烧蚀行为模拟. 图15给出$Sh =0.01$时形貌剖视图,可以发现不同碳纤维倾斜角的材料微观烧蚀行为有很大差异:随着碳纤维倾斜角的增大,烧蚀速率显著降低.这是因为碳纤维烧蚀速率大于基体烧蚀速率,碳纤维首先发生氧化形成氧通道,而在相同碳纤维烧蚀量的条件下,碳纤维倾斜角越小则氧通道进入材料内部越深,使得材料内部更多的暴露于氧环境下,造成更大的烧蚀速率.随着碳纤维倾斜角的增大,氧通道更多的在水平方向产生,而在沿深度的方向碳纤维由抗氧化强的基体分层包裹,使得材料烧蚀速率逐渐降低. 图16可知,随着倾斜角度的增大,材料质量损失速率逐渐降低,当$\theta = 0^ \circ $的材料完全氧化时$\theta = 90^ \circ $的材料尚有约60%质量的材料.显示原图|下载原图ZIP|生成PPT
图15不同时刻微观烧蚀形貌 ($A=0.2, Sh=0.01$)
-->Fig. 15Micro-ablation morphologies at different times ($A=0.2, Sh=0.01$)
-->
显示原图|下载原图ZIP|生成PPT
图16不同纤维倾斜角度的材料质量损失速率($A=0.2$, $Sh =0.01$)
-->Fig. 16Mass loss rate of materials at different fiber inclination angles ($A=0.2$, $Sh =0.01$)
-->
当$Sh = 1$时,由图17可知不同碳纤维倾斜角的烧蚀形态依然存在差异,但是差异已经非常小. 这是因为当$Sh =1$时氧气的 扩散速率和材料表面的碳氧反应速率处于同等量级.这时在烧蚀表面的氧浓度已经不足以供应充足的氧向材料内部大量的扩散,造成基体烧蚀表面和碳纤维烧蚀表面的高度差异不大.图18的材料质量损失速率进一步说明:相比$Sh = 0.01$时,在$Sh =1$条件下的碳纤维倾斜角对材料的烧蚀行为影响确实在减弱.
显示原图|下载原图ZIP|生成PPT
图17不同时刻微观烧蚀形貌($A=0.2$, $Sh =1$)
-->Fig. 17Micro-ablation morphologies at different times ($A=0.2$, $Sh =1$)
-->
显示原图|下载原图ZIP|生成PPT
图18不同纤维倾斜角度的材料质量损失速率($A=0.2$, $Sh =1$)
-->Fig. 18Mass loss rate of materials at different fiber inclination angles ($A=0.2$, $Sh =1$)
-->
4 结 论
对微观尺度下的碳纤维复合材料烧蚀行为进行模拟. 对单根纤维的烧蚀行为进行模拟并验证了计算程序的正确性.对有倾斜角 的多根碳纤维在不同条件下的烧蚀行为进行模拟. 主要有以下结论:(1) 当纤维抗氧化性能比基体强时,基体烧蚀速率大于纤维烧蚀速率,使纤维高于基体表面.在小$Sh$时形成"笋尖"状烧蚀形貌,并随着$Sh$逐步增大,纤维露出的高度逐渐降低;在大$Sh$时纤维和基体将同步后退.纤维倾斜角度对材料的微观烧蚀行为没有显著影响.
(2) 当纤维抗氧化性能比基体弱时,碳纤维倾斜角度对材料的微观烧蚀有显著影响. 在小$Sh$时,倾斜角度越大,材料的烧蚀速率越小,表面后退速度也降低;随着$Sh$的增大,倾斜角对烧蚀行为的影响逐渐减弱.
本文对碳纤维增强复合材料的微观烧蚀行为进行了初步的探索,对于不同温度、不同气体介质以及复杂的材料组分等情况还需要进一步研究.
The authors have declared that no competing interests exist.
参考文献 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子
[1] | . , 在30-70km空域机动飞行的高超声速飞行器的优点是可以耦合利用所处空域的空气产生的升力和高速飞行的离心力进行远距离机动滑翔飞行,具有重要的实用价值尽管过去数十年在高超声速流动研究方面取得显著进展,但在设计研究近空间远程滑翔的高超声速飞行器方面仍然存在许多挑战,特别是对特定飞行条件下的流动机理了解不清楚本文介绍了作者研究团队在开展近空间高超声速飞行器有关的关键气动问题方面的研究进展,主要包括:建立了近空间高超声速飞行的流动模型,发展了系统的相关计算空气动力学方法,针对高空高速飞行条件下稀薄气体效应和真实气体效应的耦合作用影响研究了合适的滑移边界条件,考虑了不同组分存在条件下的温度、速度和压力的滑移效应影响;提出了飞行器气动外形的动态优化方法,获得了可工程实用化的高升阻比飞行器气动外形;建立了高速飞行器动稳定性理论,在实现高超声速飞行器动态稳定飞行方面取得重大进展;最后讨论了高超声速飞行器设计中进一步需要关注的若干关键技术和科学问题、可能解决的途径及其所涉及的学科发展方向. . , 在30-70km空域机动飞行的高超声速飞行器的优点是可以耦合利用所处空域的空气产生的升力和高速飞行的离心力进行远距离机动滑翔飞行,具有重要的实用价值尽管过去数十年在高超声速流动研究方面取得显著进展,但在设计研究近空间远程滑翔的高超声速飞行器方面仍然存在许多挑战,特别是对特定飞行条件下的流动机理了解不清楚本文介绍了作者研究团队在开展近空间高超声速飞行器有关的关键气动问题方面的研究进展,主要包括:建立了近空间高超声速飞行的流动模型,发展了系统的相关计算空气动力学方法,针对高空高速飞行条件下稀薄气体效应和真实气体效应的耦合作用影响研究了合适的滑移边界条件,考虑了不同组分存在条件下的温度、速度和压力的滑移效应影响;提出了飞行器气动外形的动态优化方法,获得了可工程实用化的高升阻比飞行器气动外形;建立了高速飞行器动稳定性理论,在实现高超声速飞行器动态稳定飞行方面取得重大进展;最后讨论了高超声速飞行器设计中进一步需要关注的若干关键技术和科学问题、可能解决的途径及其所涉及的学科发展方向. |
[2] | . , 气动热风洞实验是地面研究和预测飞行器气动热环境的重要手段之一, 但由于风洞实验模拟能力的限制, 风洞实验的流场参数和模型的几何尺度都会与实际飞行情况存在一定的差别, 导致地面风洞实验中得到的模型表面气动加热率数据无法直接用于飞行条件下的热环境预测和热防护设计. 以往通过针对具体飞行器的试验结果进行数据拟合后外插的气动热关联换算方法指向性较强, 没有考虑到气动热的具体影响参数, 存在一定局限性, 难以外推应用于其他外形的飞行器. 为解决通过气动热风洞实验数据外推预测飞行条件下气动热的技术难题, 基于无量纲NS方程和边界层理论分析研究了影响气动热的主要参数, 并通过推导化简边界层近似解热流公式, 针对层流流态建立了气动热关联换算方法, 可以考虑当地边界层外缘参数的影响, 具有一定通用性. 在此基础上, 利用建立的方法将Reentry-F飞行器缩比模型的风洞实验数据换算到该飞行器飞行条件下的典型工况, 并与飞行测量结果进行了比较, 外推预测结果与飞行测量结果符合较好, 表明建立的关联方法可以用于气动热风洞实验数据的外推换算. . , 气动热风洞实验是地面研究和预测飞行器气动热环境的重要手段之一, 但由于风洞实验模拟能力的限制, 风洞实验的流场参数和模型的几何尺度都会与实际飞行情况存在一定的差别, 导致地面风洞实验中得到的模型表面气动加热率数据无法直接用于飞行条件下的热环境预测和热防护设计. 以往通过针对具体飞行器的试验结果进行数据拟合后外插的气动热关联换算方法指向性较强, 没有考虑到气动热的具体影响参数, 存在一定局限性, 难以外推应用于其他外形的飞行器. 为解决通过气动热风洞实验数据外推预测飞行条件下气动热的技术难题, 基于无量纲NS方程和边界层理论分析研究了影响气动热的主要参数, 并通过推导化简边界层近似解热流公式, 针对层流流态建立了气动热关联换算方法, 可以考虑当地边界层外缘参数的影响, 具有一定通用性. 在此基础上, 利用建立的方法将Reentry-F飞行器缩比模型的风洞实验数据换算到该飞行器飞行条件下的典型工况, 并与飞行测量结果进行了比较, 外推预测结果与飞行测量结果符合较好, 表明建立的关联方法可以用于气动热风洞实验数据的外推换算. |
[3] | . , . , |
[4] | . , ABSTRACT The reactivity of fine-weave pierced 3D carbon-carbon composites in air at temperatures up to 3000 °C was studied. The corrosion morphology and microstructure of oxidized samples were investigated by XPS, SEM, and XRD techniques, and the non-equilibrium nature of the oxidation process was pointed out. A thermochemical ablation model of CC composites controlled by gas phase diffusion and reaction kinetics was developed. |
[5] | . , A lightweight carbon-phenolic ablator, with a density of 0.5聽g/cm 3 , designed to be used as a thermal protection system for a re-entry space vehicle, was manufactured by infiltration of a carbon felt with a phenolic resin. A sample of this ablative material was tested in a Plasma Wind Tunnel (PWT) facility, simulating erosion and heat flux conditions consistent with an orbital reentry. The surface temperature of the test article was monitored during the PWT test. Microstructural and microtomographic analyses were carried out on the tested sample to investigate the effect of the high heat flux exposure on the composite material, by measuring the amount of ablation and the depth of pyrolyzation. Moreover a finite element model was implemented in order to rebuild the PWT test. Very encouraging results were obtained in terms of surface insulation capacity and surface recession. The pyrolysis and erosion of the ablator was simulated by implementing a complex finite element model, with results in very good agreement with experimental evidences. |
[6] | , Publication » Phenolic Impregnated Carbon Ablators (PICA) for Discovery class missions. |
[7] | . , This pioneering book presents new models for the thermomechanical behavior of composite materials and structures taking into account internal physico-chemical transformations such as thermodecomposition, sublimation and melting at high temperatures (up to 3000 K). It is of great importance for the design of new thermostable materials and for the investigation of reliability and fire safety of composite structures. It also supports the investigation of interaction of composites with laser irradiation and the design of heat-shield systems. Structural methods are presented for calculating the effective mechanical and thermal properties of matrices, fibres and unidirectional, reinforced by dispersed particles and textile composites, in terms of properties of their constituent phases. Useful calculation methods are developed for characteristics such as the rate of thermomechanical erosion of composites under high-speed flow and the heat deformation of composites with account of chemical shrinkage. The author expansively compares modeling results with experimental data, and readers will find unique experimental results on mechanical and thermal properties of composites under temperatures up to 3000 K. Chapters show how the behavior of composite shells under high temperatures is simulated by the finite-element method and so cylindrical and axisymmetric composite shells and composite plates are investigated under local high-temperature heating. The book will be of interest to researchers and to engineers designing composite structures, and invaluable to materials scientists developing advanced performance thermostable materials.Dimitrienko, Yu I |
[8] | , The Mars Science Laboratory heat shield was designed to withstand a fully turbulent heat pulse using information from ground testing and computational analysis on a preflight design trajectory. Instrumentation on the flight heat shield measured in-depth temperatures to permit reconstruction of the surface heating. The data indicate that boundary-layer transition occurred at five of seven measurement locations before peak heating. Data oscillations at three pressure measurement locations may also indicate transition. This paper presents the heat shield temperature and pressure data, possible explanations for the timing of boundary-layer transition, and a comparison of reconstructed and computational heating on the actual trajectory. A smooth-wall boundary-layer Reynolds number that was used to predict transition is compared with observed transition at various heat shield locations. A single transition Reynolds number criterion does not uniformly explain the timing of boundary-layer transition observed duri... |
[9] | |
[10] | . , |
[11] | . , http://arc.aiaa.org/doi/abs/10.2514/2.3469 |
[12] | . , Abstract The Porous-material Analysis Toolbox based on OpenFOAM-extend (PATO) is a modular analysis platform specifically implemented to test innovative physics-based models for porous materials submitted to high-temperature environments or other unusual conditions. The governing equations implemented in the different modules are volume-averaged forms of the mass-, momentum-, and energy-conservation equations for porous media. Current developments are focussed on ablative materials. PATO is currently composed of two types of modules: (1) global analysis modules, that may be used to run a full ablative material response, with an applied/macroscopic scale point of view; (2) elementary analysis modules, that may be used to study specific fundamental aspects, with a detailed/microscopic scale point of view. Examples of global and elementary analysis applications are presented: volume-averaged simulation of the oxidation of the fibers in a carbon-fiber preform, multidimensional pyrolysis-gas flow in a cylinder facing an arc-jet, comparison of equilibrium and finite-rate chemistry in a carbon/phenolic ablative material. |
[13] | . , In-depth thermal responses and surface recession of charring materials are two main performances of the thermal protection system in reentry vehicles subjected to aerodynamic heating. To explore the thermal behavior of the charring ablator, we developed a new one-dimensional pyrolysis layer model with the heat and mass transfer through charring materials undergoing decomposition and surface ablation. Taking AVCOAT composites in an Apollo capsule during reentry as an example, the calculation codes are written to simulate its nonlinear thermal behavior caused by temperature dependent thermal properties, two moving interfaces and one moving boundary. Numerical results indicate that the pyrolysis layer model is suitable for solving the problem of the thermal behavior of the charring ablator. The nonlinear analysis is essential in the calculation process although it brings oscillations, which have little influence on surface ablation and decomposition. The thermal blockage effect of the pyrolysis gas from decomposition is significant in TPS during earth entry. This study will be helpful for the design of the thermal protection system in reentry vehicles. |
[14] | |
[15] | |
[16] | . , 为了更好地研究三维整体编织碳 /酚醛复合材料的成型工艺对材料烧蚀性能的影响 ,选用不同纤维体积分数和不同编织结构的碳 /酚醛复合材料试样 ,进行烧蚀试验 ,利用TalyScan 15 0型表面粗糙度测试仪对试样烧蚀后的表面进行测试 ,并采用多种分析方法对测试结果进行分析。从分析结果可以看出三维四向结构的碳 /酚醛复合材料随着纤维体积分数的增加 ,烧蚀性能变好 ,较低纤维体积分数 (5 0 % )的三维五向结构碳 /酚醛复合材料具有较好的烧蚀性能。 . , 为了更好地研究三维整体编织碳 /酚醛复合材料的成型工艺对材料烧蚀性能的影响 ,选用不同纤维体积分数和不同编织结构的碳 /酚醛复合材料试样 ,进行烧蚀试验 ,利用TalyScan 15 0型表面粗糙度测试仪对试样烧蚀后的表面进行测试 ,并采用多种分析方法对测试结果进行分析。从分析结果可以看出三维四向结构的碳 /酚醛复合材料随着纤维体积分数的增加 ,烧蚀性能变好 ,较低纤维体积分数 (5 0 % )的三维五向结构碳 /酚醛复合材料具有较好的烧蚀性能。 |
[17] | . , 本文对宏观热化学烧蚀理论进行总结和归纳,应用这种理论分别推导了碳-碳复合材科和碳化钨-碳复合材料烧蚀性能预测公式,并与电弧加热器的试验结果进行了比较,两者的差别在工程允许范围内。 . , 本文对宏观热化学烧蚀理论进行总结和归纳,应用这种理论分别推导了碳-碳复合材科和碳化钨-碳复合材料烧蚀性能预测公式,并与电弧加热器的试验结果进行了比较,两者的差别在工程允许范围内。 |
[18] | . , Ablation of carbon-based materials is a key issue in atmospheric reentry; it displays a strong coupling between mass, momentum and heat transfers, the importance of which relies on the surface roughness. A new possible physical cause for roughness set-up is investigated, based on the coupling between diffusive transfer in the surrounding fluid on one hand, and heterogeneous reaction or sublimation on the other. Considering mass transfer in a 2D, isothermal, vertical-flux approximation, the surface is proved to be able to acquire, among others, a stable stationary morphology made of circle arcs connected by symmetrical singular points. Such a morphology has indeed been observed in the case of graphite ablation, and the computed roughness length scale, arising from the diffusion-to-reaction ratio, is compatible with observed data. A similar model based on the presence of a thermal gradient yields similar results, but with a larger length scale, also compatible with other observations. |
[19] | , ABSTRACT In this paper, we are interested in the ablation of carbon/carbon composites, as found, for instance, in rocket nozzle applications, where aggressive gases cause the composite surface recession. Global recession velocity and the appearing roughness are function of the subscale description. We model ablation in terms of surface oxidation of a heterogeneous material by a binary gas mixture. We assume the gas transport to be purely diffusive and temperature almost uniform.A finite volume code has been developed based on a VOF method using complex representation of interface and elementary surfaces, and sequential solving of the coupled equations for fluid transport and surface recession. The code was used to perform many numerical experiments, which allowed us to identify different recession regimes as a function of the characteristic dimensionless numbers.A one-dimensional model was built using the concept of effective surface reaction. This approach replaces the non-uniform reactivity of a geometrically complex surface by an effective reactivity of a plane homogeneous interface that gives macro-scale fluxes similar to the average micro-scale fluxes. |
[20] | . , The modelling of ablation of carbon/carbon (/) composites utilized as rocket engine hot parts is addressed under the angle of the competition between bulk transport of reactants and heterogeneous mass transfer, associated to reactivity contrasts between constituent phases. A numerical solver based on a simple model and built on a VOF technique allows direct simulation at two scales. Its application to actual complex materials is performed; the results are consistent with experimental data and help understanding the origin of the material behaviour, either in terms of acquired surface morphology or in terms of effective recession rate. |
[21] | . , The oxidation reactivities of two C/C composites and of their available components (fibers, bulk matrices) are determined by measurement of mass loss rate in a cylindrical oxidation reactor under dry air and at atmospheric pressure. In order to identify reactional and diffusional regimes, and to provide a safe method for the identification of the intrinsic heterogeneous reaction rates at high temperatures, a modeling approach has been developed. Diffusion of the oxidant is considered throughout the reactor for all the samples (global-scale modeling) in combination with convection and reactions. Fibers have been arranged in unidirectional bundles in which diffusion is also considered (local-scale modeling). The importance of the reaction rate relatively to global-scale and local-scale diffusion has been evaluated. When reaction is slow enough, diffusion effects can be neglected; in the converse case, the intrinsic reaction rates are extracted from the experimental data using the models. Incidentally, the comparison of the fiber and matrix intrinsic reactivities supports the idea that the composites feature a complex oxidation behavior mainly based on a weakest link process. This result is illustrated and discussed using SEM and TEM investigations. [All rights reserved Elsevier]. |
[22] | . , Following an analysis of surface roughness features that develop on a 3D C/C composite during ablation, i.e. wall recession by oxidation and/or sublimation, a modeling strategy is set up in order to predict the composite behavior from that of its components. It relies on two changes of scale: (i) microscopic scale (fiber, matrix) to mesoscopic scale (bundle) and (ii) mesoscopic scale (bundle, matrix) to macroscopic scale (composite). The physical basis is a general model for receding surfaces under a gasification process coupled to mass transfer. At each scale, the 3D surface equation is analytically solved in steady state considering a 1-D mass transfer perpendicular to the overall surface. The models are validated by comparison to experimental data. |
[23] | . , Ablation of carbon鈥揷arbon composites (C/C) results in a heterogeneous surface recession mainly due to some gasification processes (oxidation, sublimation) possibly coupled to bulk mass transfer. In order to simulate and analyse the material/environment interactions during ablation, a Brownian motion simulation method featuring special Random Walk rules close to the wall has been implemented to efficiently simulate mass transfer in the low P茅clet number regime. A sticking probability law adapted to this kind of Random Walk has been obtained for first-order heterogeneous reactions. In order to simulate the onset of surface roughness, the interface recession is simultaneously handled in 3D using a Simplified Marching Cube discretization. This tool is validated by comparison to analytical models. Then, its ability to provide reliable and accurate solutions of ablation phenomena in 3D is illustrated. |
[24] | . , 将DSMC方法运用在材料细观烧蚀机理分析中,对碳/碳材料在典型烧蚀过程中的简化模型进行了微米量级下的氧扩散特性分析,旨在研究氧分子在材料表面缝隙中的扩散特性及其与壁面的作用规律,计算分析结果验证了纤维烧蚀结构演变特性。结果表明:沿缝隙深度方向,O2与壁面的碰撞频率不断降低;在缝隙入口处烧蚀最快,缝隙深处最慢,碳纤维随着烧蚀过程的进行不断尖化直至达到强度极限产生折断剥离;随着烧蚀过程的不断进行,O2与壁面的碰撞频率增加,材料将烧蚀得更快;在同一烧蚀条件下,O2与壁面发生烧蚀反应的概率越大,将消耗O2越快,致使O2与壁面的碰撞频率降低。 . , 将DSMC方法运用在材料细观烧蚀机理分析中,对碳/碳材料在典型烧蚀过程中的简化模型进行了微米量级下的氧扩散特性分析,旨在研究氧分子在材料表面缝隙中的扩散特性及其与壁面的作用规律,计算分析结果验证了纤维烧蚀结构演变特性。结果表明:沿缝隙深度方向,O2与壁面的碰撞频率不断降低;在缝隙入口处烧蚀最快,缝隙深处最慢,碳纤维随着烧蚀过程的进行不断尖化直至达到强度极限产生折断剥离;随着烧蚀过程的不断进行,O2与壁面的碰撞频率增加,材料将烧蚀得更快;在同一烧蚀条件下,O2与壁面发生烧蚀反应的概率越大,将消耗O2越快,致使O2与壁面的碰撞频率降低。 |
[25] | . , |
[26] | |
[27] | . , 烧蚀材料在高温环境中因物理、化学和力学因素造成质量损失,影响其性能。建立相应的热化学和 热力学烧蚀模型,分析其传热特性和力学性能的演化规律,建立材料的热毁损判据,是复合材料防热结构设计所面临的重要问题。本文综述了近年来文献中对炭基复 合材料烧蚀机理的研究,并对几种烧蚀模型进行了评价。 . , 烧蚀材料在高温环境中因物理、化学和力学因素造成质量损失,影响其性能。建立相应的热化学和 热力学烧蚀模型,分析其传热特性和力学性能的演化规律,建立材料的热毁损判据,是复合材料防热结构设计所面临的重要问题。本文综述了近年来文献中对炭基复 合材料烧蚀机理的研究,并对几种烧蚀模型进行了评价。 |
[28] | . , The analysis of carbon oxidation data presented in a previous Fuel paper is shown to contain an error, as a result of which the intrinsic reactivity of carbon to oxygen is under-estimated by a factor of between one and four. |
[29] | . , A method is presented to solve two-phase problems involving a material quantity on an interface. The interface can be advected, stretched, and change topology, and material can be adsorbed to or desorbed from it. The method is based on the use of a diffuse interface framework, which allows a simple implementation using standard finite-difference or finite-element tech niques. Here, finite-difference methods on a block-structured adaptive grid are used, and the resulting equations are solved using an on-linear multigrid method. Interfacial flow with soluble surfactants is used as an example of the application of the method, and several test cases are presented demonstrating its accuracy and convergence. |
[30] | . , Motivated by the need for three-dimensional methods for interface calculations that can deal with topology changes, we describe a numerical scheme, built from a volume-of-fluid interface tracking technique that uses a piecewise-linear interface calculation in each cell. Momentum balance is computed using explicit finite volume/finite differences on a regular cubic grid. Surface tension is implemented by the continuous surface stress or continuous surface force method. Examples and verifications of the method are given by comparing simulations to analytical results and experiments, for sedimenting droplet arrays and capillary waves at finite Reynolds number. In the case of a pinching pendant drop, both three-dimensional and axisymmetric simulations are compared to experiments. Agreement is found both before and after the reconnections. |
[31] | . , 采用亚迭代耦合求解流体动力学方程和刚体动力学方程数值研究方法,建立流体与刚体运动耦合的无时间滞后的数值手段来分析研究飞行器非定常的耦合运动特征,预测失稳运动的发生与发展.研究表明:在流动结构失稳后机体运动逐步形成极限环振荡的周期性自维持运动,滚转力矩系数滞回曲线呈现典型的"双8"稳定形态,前缘涡的非对称飘起与再附是典型三角翼摇滚的驱动形式;三角翼构型的自由滚振幅随攻角变化呈现阶跃式的变化;自激振荡的幅度及频率与滚动惯量的大小存在对应的规律性关系.非定常自激振动的翼摇滚具有多控制参数、非稳定、非线性的运动特 . , 采用亚迭代耦合求解流体动力学方程和刚体动力学方程数值研究方法,建立流体与刚体运动耦合的无时间滞后的数值手段来分析研究飞行器非定常的耦合运动特征,预测失稳运动的发生与发展.研究表明:在流动结构失稳后机体运动逐步形成极限环振荡的周期性自维持运动,滚转力矩系数滞回曲线呈现典型的"双8"稳定形态,前缘涡的非对称飘起与再附是典型三角翼摇滚的驱动形式;三角翼构型的自由滚振幅随攻角变化呈现阶跃式的变化;自激振荡的幅度及频率与滚动惯量的大小存在对应的规律性关系.非定常自激振动的翼摇滚具有多控制参数、非稳定、非线性的运动特 |
[32] | . , A methodology for designing formally second-order time-accurate and yet loosely-coupled partitioned procedures for the solution of nonlinear fluid鈥搒tructure interaction (FSI) problems on moving grids is presented. Its key components are a fluid time-integrator that is provably second-order time-accurate on moving grids, the midpoint rule for advancing in time the solution of the structural dynamics equations of motion, a second-order structure predictor for bypassing the inner-iterations encountered in strongly-coupled solution procedures, and a carefully designed algorithm for time-integrating the motion of the fluid-mesh. Following this methodology, two different loosely-coupled schemes are constructed for the solution of transient nonlinear FSI problems and proved to be second-order time-accurate. Three-dimensional numerical results pertaining to the simulation of the aeroelastic response to a gravity excitation of a complete F-16 configuration are also presented. In addition to confirming the theoretical results discussed in this paper, these numerical results highlight a very stable behavior of the designed loosely-coupled partitioned procedures. |
[33] | . , . , |