CN110348165B

CN110348165B - 基于结构网格的陶瓷基复合材料细观建模与力学计算方法

Info

Publication number: CN110348165B
Application number: CN201910653580.7A
Authority: CN
Inventors: 高希光; 董洪年; 宋迎东; 张盛; 于国强; 陈正扬
Original assignee: Nanjing University of Aeronautics and Astronautics
Current assignee: Nanjing University of Aeronautics and Astronautics
Priority date: 2019-07-18
Filing date: 2019-07-18
Publication date: 2020-06-30
Anticipated expiration: 2039-07-18
Also published as: CN110348165A

Abstract

本发明属于复合材料建模与力学计算技术领域，具体涉及基于结构网格的陶瓷基复合材料细观建模与力学计算方法。本发明基于结构网格，由给定的独立的几个细观几何参数，区分像素单元所属材料类别，从而实现陶瓷基复合材料细观建模。采用均匀化理论，依据六组不同周期性边界条件下的应力应变关系，获得柔度矩阵，最终实现弹性参数计算。本发明参数化地实现陶瓷基复合材料细观几何模型参数化建模及六面体网格自动划分，并实现力学计算。与传统方法相比，大大降低了难度，减少了对人的依赖，缩减了建模与计算的时间。

Description

基于结构网格的陶瓷基复合材料细观建模与力学计算方法

技术领域

本发明属于复合材料建模与力学计算技术领域，具体涉及基于结构网格的陶瓷基复合材料细观建模与力学计算方法。

背景技术

陶瓷基复合材料作为一种先进的新型材料已经成为航空航天等特种工程技术行业中的理想备选材料。为了实现陶瓷基复合材料的结构设计及工程应用，必须要对其进行细观建模及力学计算。目前，常见的细观建模方法可以分为非结构网格法和结构网格法。

基于非结构网格，孔春元等(孔春元，孙志刚，高希光，宋迎东：2.5维C/SiC复合材料单胞模型及刚度预测，26，2459-2467，(2011))实现了2.5维编织陶瓷基复合材料的细观建模与刚度预测。然而，在非结构化网格中，空间拓扑关系将随着几何参数的变化而发生变化，从而极可能引起建模失败。此外，在网格划分方面，尚无法实现六面体网格的自动划分。这使得细观建模成为一件繁琐而耗时的工作。

参数化建模一直是工程技术人员不断追求的目标。结构化网格区域内所有的内部点都具有相同的毗邻单元，可以很容易地实现区域的边界拟合，具有网格生成速度快，数据结构简单等优点，可以容易地实现复杂结构的参数化建模及自动网格划分。中国专利CN103871059B《纤维增强复合材料等效弹性参数的计算方法》提供了一种纤维增强复合材料等效弹性参数的计算方法，但是在进行细观建模时采用了XCT技术，而不是参数化的建模方法。即其只能对事先制备完成的材料进行细观建模与力学计算，而无法通过给定任意参数而实现材料的高通量细观建模。目前尚没有一种有效的基于结构网格的陶瓷基复合材料细观建模与力学计算方法。

发明内容

本发明针对现有技术中的不足，提供一种基于结构网格的陶瓷基复合材料细观建模与力学计算方法，可实现陶瓷基复合材料细观几何模型参数化建模及六面体网格自动划分，并实现力学计算。

为实现上述目的，本发明采用以下技术方案：

基于结构网格的陶瓷基复合材料细观建模与力学计算方法，其特征在于，包括如下步骤：

步骤一：计算建模所需几何参数；

步骤二：根据几何参数生成单胞，并划分结构网格；

步骤三：区分单胞中的像素单元材料；

步骤四：设置像素单元材料属性；

步骤五：对单胞施加周期性边界条件并计算柔度矩阵；

步骤六：根据柔度矩阵计算材料的弹性参数。

为优化上述技术方案，采取的具体措施还包括：

进一步地，所述步骤一具体如下：

1)给定独立的几何参数，包括经纱走向波长LC、经纱走向幅值AC和经纱截面宽WC；

2)计算出建模所需的所有几何参数：假设纱线截面由两条形状相同的正弦曲线组成，其形状由正弦曲线波长λ_c＝LC和幅值A_c＝AC，以及短轴间距L_c＝AC决定；几何参数还包括单胞长L₁＝LC，单胞宽L₂＝2·WC，单胞高H＝2·AC。

进一步地，所述步骤二中，划分结构网格具体如下：

分别在长、宽、高三个维度，将单胞进行N1、N2、N3等分；由此，整个单胞被划分成N1*N2*N3个像素单元，同时，在单个长、宽、高的边上，分别出现N1+1、N2+1、N3+1个节点，整个单胞将含有(N1+1)*(N2+1)*(N3+1)个节点。

进一步地，所述步骤三中，采用线性规划的方法区分经纱及纬纱区域，具体如下：

1)二维纱线截面中的区域划分

通过比较各像素单元中心坐标与边界曲线之间的关系，判断该像素单元是否位于闭合区域内；

对于经纱截面而言，二维经纱截面中的区域划分为

其中(y，z)为像素单元中心坐标，f₁(y)为上边界曲线函数，f₂(y)为下边界曲线函数，S(y，z)为坐标的状态函数，S(y，z)＝1表示给定坐标在闭合区域内，属于经纱区域，S(y，z)＝0表示给定坐标在闭合区域之外，不属于经纱区域；

当截面中心坐标为(y＝0，z＝0)时

f₁(y)＝L_c/2

f₂(y)＝-L_c/2

对于纬纱截面而言，二维纬纱截面中的区域划分为

其中(x，z)为像素单元中心坐标，f₁(x)为上边界曲线函数，f₂(x)为下边界曲线函数，S(x，z)为坐标的状态函数，S(x，z)＝1表示给定坐标在闭合区域内，属于纬纱区域，S(x，z)＝0表示给定坐标在闭合区域之外，不属于纬纱区域；

当截面中心坐标为(x＝0，z＝0)时

f₁(x)＝A_c·cos((2π/λ_c)·x)-(A_c-L_c/2)

f₂(x)＝A_c·cos((2π/λ_c)·x+π)+(A_c-L_c/2)

2)考虑纱线走向的纱线高度修正

对于经纱走向而言，当x坐标发生变化时，即沿经纱走向位置发生变化时，经纱截面位置以正弦曲线的形式上下波动，则

其中，(x，y，z)为考虑经纱走向的像素单元中心坐标，f₃(x)为考虑经纱走向的正弦曲线，函数S(x，y，z)为坐标的状态函数；

f₃(x)＝AC·(sin((2π/LC)·x+π/2)-1)

对于纬纱走向而言，当x坐标发生变化时，即沿纬纱走向位置发生变化时，纬纱截面位置上下保持不变，则

其中，(x，y，z)为考虑经纱走向的像素单元中心坐标，f₃(y)为考虑经纱走向的正弦曲线函数；

f₃(y)＝0。

进一步地，所述步骤四中，基于像素单元所在位置信息，不同的像素单元被赋予不同的材料属性，具体如下：

经纱单元被赋予经纱材料属性；纬纱单元被赋予纬纱材料属性；孔隙单元被赋予基体材料属性；经纱被视为横观各向同性材料；纬纱被视为横观各向同性材料，其材料参数等同于经纱的材料参数，区别仅在于材料主方向的不同；孔隙被视为各向同性材料。

进一步地，所述步骤五具体如下：

对单胞施加周期性位移边界条件：

其牛

是单胞的平均应变，j+和j-表示垂直于x_j轴的一对边界面，u表示位移，Δx表示相对点的坐标差；

给定六组不同的周期性边界条件，然后分别进行有限元分析；根据均匀化理论，定义材料的宏观应力

为材料细观应力的体积平均，材料的宏观应变

由周期性边界条件直接得到：

其中V是单胞体积，σⁱ是第i个像素单元的应力；计算出六组宏观应力和宏观应变，六组宏观应力和宏观应变具有如下的关系：

其中上标i(i＝1，…，6)表示第i组周期性边界条件下的结果，下标ij(i，j＝1，…，6)表示材料的主方向，S_ij(i，j＝1，…，6)为柔度矩阵系数；

由上式可得柔度矩阵表达式为

本发明的有益效果是：可以由给定的独立的几个细观几何参数，参数化地实现陶瓷基复合材料细观几何模型参数化建模及六面体网格自动划分，并实现力学计算。本发明与传统方法相比，大大降低了难度，减少了对人的依赖，缩减了建模与计算的时间。

附图说明

图1是结构网格划分结果示意图。

图2是二维经纱截面中的区域划分示意图。

图3是二维纬纱截面中的区域划分示意图。

图4是经纱走向示意图。

图5是纬纱走向示意图。

具体实施方式

现在结合附图对本发明作进一步详细的说明。

实施案例：基于结构网格的平纹编织SiC/SiC陶瓷基复合材料细观建模与刚度计算。

步骤一：计算建模所需几何参数

(1)给定独立的几何参数，包括经纱走向波长LC＝10mm、经纱走向幅值AC＝0.24mm和经纱截面宽WC＝1.74mm。

(2)计算出建模所需的所有几何参数

假设纱线截面由两条形状相同的正弦曲线组成，其形状由正弦曲线波长λ_c＝LC＝10mm和幅值A_c＝AC＝0.24mm，以及短轴间距L_c＝AC＝1.74mm决定。除此之外，还需要的几何参数有单胞长L₁＝LC＝10mm，单胞宽L₂＝2·WC＝3.48mm，单胞高H＝2·AC＝0.48mm。

步骤二：划分结构网格

如图1所示，分别在长、宽、高三个维度，将单胞进行48、8、24等分。由此，整个单胞将被划分成48*8*24＝9216个像素单元。同时，在单个长、宽、高的边上，将分别出现48+1、8+1、24+1个节点。由此，整个单胞将含有(48+1)*(8+1)*(24+1)＝11025个节点。

步骤三：区分像素单元材料

为了区分经纱及纬纱区域，采用了线性规划的方法。

(1)二维纱线截面中的区域划分

通过比较各像素单元中心坐标与边界曲线之间的关系，判断该像素单元是否位于闭合区域内。

对于经纱截面而言，如图2所示，给出了二维经纱截面中的区域划分示意图。

其中(y，z)为像素单元中心坐标，y指该点在y坐标上的值，z指该点在z坐标上的值，f₁(y)为上边界曲线函数，f₂(y)为下边界曲线函数，S(y，z)为坐标的状态函数，S(y，z)＝1表示给定坐标在闭合区域内，属于经纱区域，S(y，z)＝0表示给定坐标在闭合区域之外，不属于经纱区域。

当截面中心坐标为(y＝0，z＝0)时

f₁(y)＝L_c/2＝0.87mm

f₂(y)＝-L_c/2＝-0.87mm

对于纬纱截面而言，如图3所示，给出了二维纬纱截面中的区域划分示意图。

其中(x，z)为像素单元中心坐标，x指该点在x坐标上的值，z指该点在z坐标上的值，f₁(x)为上边界曲线函数，f₂(x)为下边界曲线函数，S(x，z)为坐标的状态函数，S(x，z)＝1表示给定坐标在闭合区域内，属于纬纱区域，S(x，z)＝0表示给定坐标在闭合区域之外，不属于纬纱区域。

当截面中心坐标为(x＝0，z＝0)时

f₁(x)＝A_c·cos((2π/λ_c)·x)-(A_c-L_c/2)＝0.24·cos((2π/10)·x)-(0.24-0.87)mm

f₂(x)＝A_c·cos((2π/λ_c)·x+π)+(A_c-L_c/2)＝0.24·cos((2π/10)·x+π)+(0.24-0.87)mm

(2)考虑纱线走向的纱线高度修正

由于纱线具有起伏波动的特征，纱线截面并不总是处在同一高度水平上。

对于经纱走向而言，如图4所示，当x坐标发生变化时，即沿经纱走向位置发生变化时，经纱截面位置以正眩曲线的形式上下波动。

其中，(x，y，z)为考虑经纱走向的像素单元中心坐标，f₃(x)为考虑经纱走向的正弦曲线函数。

f₃(x)＝AC·(sin((2π/LC)·x+π/2)-1)＝0.24·(sin((2π/10)·x+π/2)-1)mm

对于纬纱走向而言，如图5所示，当x坐标发生变化时，即沿纬纱走向位置发生变化时，纬纱截面位置上下保持不变。

其中，(x，y，z)为考虑经纱走向的像素单元中心坐标，f₃(y)为考虑经纱走向的正弦曲线函数。

f₃(y)＝0

步骤四：设置像素单元材料属性

基于像素单元所在位置信息，不同的单元被赋予了不同的材料属性及编号。其中，经纱单元被赋予经纱材料属性，材料编号为1；纬纱单元被赋予纬纱材料属性，材料编号为2；孔隙单元被赋予基体材料属性，材料编号为3。

经纱被视为横观各向同性材料，其材料参数为：E₁＝241GPa，E₂＝210GPa，μ₁₂＝0.2，μ₂₃＝0.3，G₁₂＝60.5GPa。

纬纱被视为横观各向同性材料，其材料参数等同于经纱的材料参数，区别仅在于材料主方向的不同。

孔隙被视为各向同性材料，其材料参数为：E＝1GPa，μ＝0.2。

步骤五：施加周期性边界条件并计算柔度矩阵

对单胞施加周期性位移边界条件：

其中

是单胞的平均应变，j+和j-表示垂直于x_j轴的一对边界面，u表示位移，Δx表示相对点的坐标差。

如表1所示，给定六组不同的周期性边界条件。

表1六组不同的周期性边界条件

然后分别进行有限元分析。根据均匀化理论，定义材料的宏观应力

为材料细观应力的体积平均，材料的宏观应变

可由周期性边界条件直接得到：

其中V是单胞体积，σⁱ是第i个单元的应力。计算出六组宏观应力和宏观应变：

其中上标i(i＝1，…，6)表示第i组周期性边界条件下的结果，下标ij(i，j＝1，…，6)表示材料的主方向，S_ij(i，j＝1，…，6)为柔度矩阵系数。

由上式可得柔度矩阵表达式为

步骤六：计算弹性参数

由柔度矩阵系数可得材料的弹性参数：

μ₁₂＝-E₁·S₁₂＝-152.1×10⁹·(-0.000000000001075)＝0.164，

μ₁₃＝-E₁·S₁₃＝-152.1×10⁹·(-0.000000000001239)＝0.188，

μ₂₃＝-E₂·S₂₃＝-158.6×10⁹·(-0.000000000001675)＝0.266

至此，材料的弹性参数便已计算得到，其中E₁、E₂、E₃分别为材料在1、2和3方向的弹性模量，G₁₂、G₁₃、G₂₃分别为材料在12、13和23方向的剪切模量，μ₁₂、μ₁₃、μ₂₃分别为材料在12、13和23方向的泊松比。上述各方向均为业内规定的方向。

需要注意的是，发明中所引用的如“上”、“下”、“左”、“右”、“前”、“后”等的用语，亦仅为便于叙述的明了，而非用以限定本发明可实施的范围，其相对关系的改变或调整，在无实质变更技术内容下，当亦视为本发明可实施的范畴。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims

1.基于结构网格的陶瓷基复合材料细观建模与力学计算方法，其特征在于，包括如下步骤：

步骤一：计算建模所需几何参数；所述步骤一具体如下：

2)计算出建模所需的所有几何参数：假设纱线截面由两条形状相同的正弦曲线组成，其形状由正弦曲线波长λ_c＝LC和幅值A_c＝AC，以及短轴间距L_c＝AC决定；几何参数还包括单胞长L₁＝LC，单胞宽L₂＝2·WC，单胞高H＝2·AC；

步骤二：根据几何参数生成单胞，并划分结构网格；

步骤三：区分单胞中的像素单元材料；所述步骤三中，采用线性规划的方法区分经纱及纬纱区域，具体如下：

1)二维纱线截面中的区域划分

对于经纱截面而言，二维经纱截面中的区域划分为

其中(y,z)为像素单元中心坐标，f₁(y)为上边界曲线函数，f₂(y)为下边界曲线函数，S(y,z)为坐标的状态函数，S(y,z)＝1表示给定坐标在闭合区域内，属于经纱区域，S(y,z)＝0表示给定坐标在闭合区域之外，不属于经纱区域；

当截面中心坐标为(y＝0,z＝0)时

f₁(y)＝L_c/2

f₂(y)＝-L_c/2

对于纬纱截面而言，二维纬纱截面中的区域划分为

其中(x,z)为像素单元中心坐标，f₁(x)为上边界曲线函数，f₂(x)为下边界曲线函数，S(x,z)为坐标的状态函数，S(x,z)＝1表示给定坐标在闭合区域内，属于纬纱区域，S(x,z)＝0表示给定坐标在闭合区域之外，不属于纬纱区域；

当截面中心坐标为(x＝0,z＝0)时

f₁(x)＝A_c·cos((2π/λ_c)·x)-(A_c-L_c/2)

f₂(x)＝A_c·cos((2π/λ_c)·x+π)+(A_c-L_c/2)

2)考虑纱线走向的纱线高度修正

其中，(x,y,z)为考虑经纱走向的像素单元中心坐标，f₃(x)为考虑经纱走向的正弦曲线，函数S(x,y,z)为坐标的状态函数；

f₃(x)＝AC·(sin((2π/LC)·x+π/2)-1)

其中，(x,y,z)为考虑经纱走向的像素单元中心坐标，f₃(y)为考虑经纱走向的正弦曲线函数；

f₃(y)＝0

步骤四：设置像素单元材料属性；

步骤五：对单胞施加周期性边界条件并计算柔度矩阵；

步骤六：根据柔度矩阵计算材料的弹性参数。

2.如权利要求1所述的基于结构网格的陶瓷基复合材料细观建模与力学计算方法，其特征在于：所述步骤二中，划分结构网格具体如下：

3.如权利要求1所述的基于结构网格的陶瓷基复合材料细观建模与力学计算方法，其特征在于：所述步骤四中，基于像素单元所在位置信息，不同的像素单元被赋予不同的材料属性，具体如下：

4.如权利要求1所述的基于结构网格的陶瓷基复合材料细观建模与力学计算方法，其特征在于：所述步骤五具体如下：

对单胞施加周期性位移边界条件：

其中

为材料细观应力的体积平均，材料的宏观应变

由周期性边界条件直接得到：

其中上标i，i＝1,…,6表示第i组周期性边界条件下的结果，下标ij，i,j＝1,…,6表示材料的主方向，S_ij，i,j＝1,…,6为柔度矩阵系数；

由上式可得柔度矩阵表达式为

5.如权利要求4所述的基于结构网格的陶瓷基复合材料细观建模与力学计算方法，其特征在于：所述步骤六中，由柔度矩阵系数得到材料的弹性参数：

μ₁₂＝-E₁·S₁₂,μ₁₃＝-E₁·S₁₃,μ₂₃＝-E₂·S₂₃

其中E₁、E₂、E₃分别为材料在1、2、3方向的弹性模量，G₁₂、G₁₃、G₂₃分别为材料在12、13、23方向的剪切模量，μ₁₂、μ₁₃、μ₂₃分别为材料在12、13、23方向的泊松比。