一种陶瓷基复合材料流固耦合响应计算方法

摘要

本发明提供一种陶瓷基复合材料流固耦合响应计算方法，方法包括：通过多尺度方法计算CMCs单胞模型加卸载应力应变迟滞曲线；通过单胞模型计算所得到的加卸载迟滞回线，插值出任意加卸载迟滞回线响应，并将其作为代理模型用于流固耦合固体域动力学计算；流固耦合界面上，流体载荷通过CFD计算，编写程序读取流体载荷并映射到固体结点上，以及读取固体结点的位移并映射到流体结点上，流体域与固体域采用相同的时间步长。本发明实现了CMCs流固耦合响应的计算，CMCs加卸载下的变刚度、迟滞行为描述简单，易于实现；动力学求解的鲁棒性好，不易发散。

权利要求

1.一种陶瓷基复合材料流固耦合响应计算方法，其特征在于，包括如下步骤：

步骤一：建立编织CMCs代表体元的有限元模型，并赋予纤维束合适的细观力学模型；

步骤二：计算编织CMCs代表体元的有限元模型的加卸载迟滞回线，插值出任意加卸载迟滞回线响应，得到任意加卸载应力应变计算方法；

步骤三：流固耦合界面上，通过CFD计算流体载荷，读取流体载荷并映射到固体结点上，并计算得到固体结点载荷；

步骤四：基于步骤二得到的任意加卸载应力应变计算方法，结合显式动力学积分和步骤三计算的固体结点载荷，求得当前时间步的CMC结构的流固耦合动力学响应；

步骤五：读取步骤四求得的流固耦合动力学响应中的固体结点位移结果并映射到流体结点上，求得耦合面流体结点位移结果，其中，流体域与固体域采用相同的时间步长；

步骤六：根据步骤五求得的流体结点位移结果对流体结点位置进行更新，并转至步骤三进行下一个时间步的CMC结构的流固耦合动力学响应计算。

2.如权利要求1所述的一种陶瓷基复合材料流固耦合响应计算方法，其特征在于：步骤二中，给予代表体元有限元模型一系列加卸载路径，加卸载过程中最大应变值逐渐增大，获得与不同最大应变值相对应的迟滞回线。

3.如权利要求1所述的一种陶瓷基复合材料流固耦合响应计算方法，其特征在于：步骤二中，利用三次多项式对迟滞回线进行拟合，拟合得到不同εi对应的多项式系数an，bn(n＝1～4)：式中，σ为应力，ε为应变，εi为第i个迟滞回线的最大应变值，+、-分别代表加载与卸载；

则任意最大应变值εt，当其处于某两个代表体元有限元模型计算所得的最大应变值之间时，即(εi＜εt＜εi+1)时，则当前迟滞回线的多项式系数插值为：当振幅由大到小时，加卸载发生在最大迟滞回线内部时，假设当前应力应变状态在P点，其位置为(εP，σP)，通过动力学数值计算，得到下一时刻的位移，则可确定相应的应变εP′，此时通过下式计算下一时刻P′点的应力值σP′：式中，A、B分别代表迟滞回线的上下顶点。

4.如权利要求1所述的一种陶瓷基复合材料流固耦合响应计算方法，其特征在于：步骤三中，通过CFD进行流体域求解，获得流固耦合面上流体单元的几何信息和载荷信息。

5.如权利要求1所述的一种陶瓷基复合材料流固耦合响应计算方法，其特征在于：步骤三中，首先进行固体单元与流体单元之间的配对，每个固体单元对应n个流体结点；然后通过下式，实现流体载荷向固体结点的映射：式中，s代表固体，f代表流体，Fsi代表作用于某固体单元第i个结点上的等效流体载荷，表示第k个流体单元作用力映射到第i个固体结点对应的等参插值系数，等参插值系数通过牛顿迭代法计算； 表示第k个流体单元对当前固体单元的作用力。

6.如权利要求1所述的一种陶瓷基复合材料流固耦合响应计算方法，其特征在于：步骤五中，首先匹配每个固体单元的某个面和面内的n个流体结点；然后通过下式，实现固体结点的位移向流体结点的映射：式中，ufj代表某流体单元上的第j个流体结点的位移， 表示固体单元第i个结点的位移映射到第j个流体结点对应的等参插值系数，等参插值系数通过牛顿迭代法计算；usi表示第i个固体结点的位移。

说明书

一种陶瓷基复合材料流固耦合响应计算方法

技术领域

[0001] 本发明属于编织陶瓷基复合材料(ceramic matrix composites,CMCs)流固耦合响应计算领域，具体涉及一种陶瓷基复合材料流固耦合响应计算方法。

背景技术

[0002] 编织CMCs在高温下，具有高比强度、高比模量、耐腐蚀等优点，是制造高温结构的理想材料。CMCs在受载时，由于陶瓷基体的脆性，材料在应变很小的情况下即会产生裂纹，因此工程应用中CMCs通常处于非线性段。在振动过程中，CMCs是处于不断加卸载情况下，会呈现复杂的变刚度、迟滞现象。此外，CMCs结构的服役环境通常伴随着流固耦合作用(fluid-structure interaction,FSI)，FSI有时甚至是导致结构失效的主要因素。因此发展一种CMCs结构的流固耦合响应的计算方法是十分必要的。

[0003] 目前，弱耦合求解方法在工程获得了广泛应用，因为可以最大程度地利用已有的固体域和流体域求解程序。然而目前在流固耦合仿真中，还少有考虑到材料加卸载下的变刚度和迟滞行为。如《一种基于交错迭代耦合技术的薄壁结构动力学热性预测方法》(中国专利CN 103177162 B)公开了一种计算薄壁流固耦合响应的方法，该专利中固体域采用Newmark法计算，未考虑材料的非线性。崔鹏、韩景龙(崔鹏,韩景龙.切尖三角翼极限环振荡的气动弹性模拟[J].航空学报,2010,31(12):2295-2301.)在切尖三角翼的流固耦合计算中，固体域考虑了材料非线性，并采用Newmark法计算动力学响应，但也未考虑材料的迟滞行为。

[0004] 前述流固耦合相关的专利和论文，固体域均未考虑材料的迟滞行为，这是因为工程应用中，设计人员通常保证材料通常工作于其线性段，就算进入的非线性段，金属材料的迟滞行为也可以忽略。然而CMCs在工程应用中，CMCs的变刚度、迟滞行为不能忽略。在振动过程中，CMCs承受的实质上是任意加卸载过程。针对任意加卸载下，CMCs的应力应变计算方法如《单向陶瓷基复合材料任意加卸载应力应变行为预测方法》(中国专利CN 104866690 B)公开了一种计算任意加卸载本构响应的计算方法，但是该方法需要反复迭代，不适合动力学计算。针对陶瓷基复合材料动力学仿真方法如《一种确定陶瓷基复合材料非线性振动响应的方法》(中国专利CN 106777595 B)公开了一种确定陶瓷基复合材料振动响应的方法，虽然该方法相较于上一个计算量大幅减少，可用于动力学计算，但该专利中在计算子迟滞回线时运算量还是较大，且在多次加卸载的情况下会积累较大的误差，从而导致动力学计算的鲁棒性较差。此外，变刚度、迟滞本构模型的引入，也给动力学方程的求解带来了极大挑战。因为加卸载下刚度会发生突变，从而引入了间断点，进而引起的动力学方程计算结果的不收敛。

[0005] 可见，目前CMCs流固耦合计算仍面临以下问题：(1)固体域缺少能够描述材料变刚度、迟滞行为的高效本构模型；(2)考虑材料变刚度、迟滞行为的固体域动力学方程求解问题。

发明内容

[0006] 本发明针对现有技术中的不足，提供一种陶瓷基复合材料流固耦合响应计算方法。

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

[0008] 一种陶瓷基复合材料流固耦合响应计算方法，其特征在于，包括如下步骤：

[0009] 步骤一：建立编织CMCs代表体元的有限元模型，并赋予纤维束合适的细观力学模型；

[0010] 步骤二：计算编织CMCs代表体元的有限元模型的加卸载迟滞回线，插值出任意加卸载迟滞回线响应，得到任意加卸载应力应变计算方法；

[0011] 步骤三：流固耦合界面上，通过CFD计算流体载荷，读取流体载荷并映射到固体结点上，并计算得到固体结点载荷；

[0012] 步骤四：基于步骤二得到的任意加卸载应力应变计算方法，结合显式动力学积分和步骤三计算的固体结点载荷，求得当前时间步的CMC结构的流固耦合动力学响应；

[0013] 步骤五：读取步骤四求得的流固耦合动力学响应中的固体结点位移结果并映射到流体结点上，求得耦合面流体结点位移结果，其中，流体域与固体域采用相同的时间步长；

[0014] 步骤六：根据步骤五求得的流体结点位移结果对流体结点位置进行更新，并转至步骤三进行下一个时间步的CMC结构的流固耦合动力学响应计算。

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

[0016] 进一步地，步骤二中，给予代表体元有限元模型一系列加卸载路径，加卸载过程中最大应变值逐渐增大，获得与不同最大应变值相对应的迟滞回线。

[0017] 进一步地，步骤二中，利用三次多项式对迟滞回线进行拟合，拟合得到不同εi对应的多项式系数an,bn(n＝1～4)：

[0018]

[0019] 式中，σ为应力，ε为应变，εi为第i个迟滞回线的最大应变值，+、-分别代表加载与卸载；

[0020] 则任意最大应变值εt，当其处于某两个代表体元有限元模型计算所得的最大应变值之间时，即(εi＜εt＜εi+1)时，则当前迟滞回线的多项式系数插值为：

[0021]

[0022] 当振幅由大到小时，加卸载发生在最大迟滞回线内部时，假设当前应力应变状态在P点，其位置为(εP,σP)，通过动力学数值计算，得到下一时刻的位移，则可确定相应的应变εP'，此时通过下式计算下一时刻P'点的应力值σP'：

[0023]

[0024] 式中，A、B分别代表迟滞回线的上下顶点。

[0025] 进一步地，步骤三中，通过CFD进行流体域求解，获得流固耦合面上流体单元的几何信息和载荷信息。

[0026] 进一步地，步骤三中，首先进行固体单元与流体单元之间的配对，每个固体单元对应n个流体结点；然后通过下式，实现流体载荷向固体结点的映射：

[0027]

[0028] 式中，s代表固体，f代表流体，Fsi代表作用于某固体单元第i个结点上的等效流体载荷， 表示第k个流体单元作用力映射到第i个固体结点对应的等参插值系数，等参插值系数通过牛顿迭代法计算； 表示第k个流体单元对当前固体单元的作用力。

[0029] 进一步地，步骤五中，首先匹配每个固体单元的某个面和面内的n个流体结点；然后通过下式，实现固体结点的位移向流体结点的映射：

[0030]

[0031] 式中，ufj代表某流体单元上的第j个流体结点的位移， 表示固体单元第i个结点的位移映射到第j个流体结点对应的等参插值系数，等参插值系数通过牛顿迭代法计算；usi表示第i个固体结点的位移。

[0032] 本发明的有益效果是：实现了CMCs流固耦合响应的计算；CMCs加卸载下的变刚度、迟滞行为的描述简单，易于实现；动力学求解的鲁棒性好，不易发散。

附图说明

[0033] 图1为本发明实施方式中插值迟滞回线应力应变曲线示意图。

[0034] 图2a为本发明实施方式中迟滞回线内卸载应力应变示意图。

[0035] 图2b为本发明实施方式中迟滞回线内加载应力应变示意图。

[0036] 图3为本发明实施方式中流体载荷映射示意图。

[0037] 图4为本发明实施方式中流固耦合位移响应示意图。

具体实施方式

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

[0039] 本发明提供了一种陶瓷基复合材料流固耦合响应计算方法，方法包括：通过多尺度方法计算CMCs单胞模型加卸载应力应变迟滞曲线；通过单胞模型计算所得到的加卸载迟滞回线，插值出任意加卸载迟滞回线响应，并将其作为代理模型用于流固耦合固体域动力学计算；流固耦合界面上，流体载荷通过CFD计算，编写程序读取流体载荷并映射到固体结点上，以及读取固体结点的位移并映射到流体结点上，流体域与固体域采用相同的时间步长。

[0040] 首先建立编织CMCs的RVE有限元模型，在RVE模型中可以加入自定义的细观本构模型，给与加卸载应变，获得加卸载本构响应曲线。在动力学计算中，直接使用RVE模型计算是极为耗时的。因此先通过RVE模型计算出一系列加卸载迟滞回线，而在动力学计算中，则在之前计算出的一系列迟滞回线之间进行插值计算。由于CMCs非线性迟滞行为与最大应变值εmax密切相关，每个迟滞回线的顶点位置和曲线形状都是由εmax决定，可通过数值拟合得到。RVE模型计算所得的第i个迟滞回线，包括加载、卸载段，可拟合成以下多项式：

[0041]

[0042] 其中，σ为应力，ε为应变，εi为第i个迟滞回线的最大应变值，+、-分别代表加载与卸载，多项式系数可表达成以下形式：

[0043]

[0044] 其中，f(εi)表示与εi相关的函数。

[0045] 在进行动力学计算过程中，首先要计算当前最大迟滞回线，若当前迟滞回线最大应变为εt，其处于RVE模型计算所得的迟滞回线对应的最大应变值为εi和εi+1之间，则当前迟滞回线的多项式系数可表达为：

[0046]

[0047] 当振幅先大后小时,则需要在迟滞回线内部插值，假设当前应力应变状态在P点，其位置为(εP,σP)。通过动力学数值计算，得到了下一时刻的位移，则可确定相应的应变εP'，此时通过下式计算下一时刻P'点的应力值σP'：

[0048]

[0049] 其中，A、B分别代表迟滞回线的上下顶点。

[0050] 通过以上过程可以计算出编织CMCs在任意加卸载下的应力应变响应，并应用于CMCs非线性动力学计算。

[0051] 在流固耦合计算中，固体域若采用Newmark等需要迭代的动力学求解方法，则CMCs的变刚度和迟滞行为，会在加卸载中引起刚度的间断，从而为动力学求解带来困难。本专利采用中心差分法等显示积分法计算振动响应，从而避免加卸载本构模型的迭代，避免了求解结果的发散。

[0052] 流固耦合的界面映射包含载荷映射和位移映射。载荷映射，在每一时间步动力学计算之前，把流体作用力映射给固体有限元结点。位移映射，在动力学计算出位移后把位移映射给流体结点。

[0053] 流体载荷映射，首先读取此时的各个耦合面单元的面心坐标以及作用力大小。由于固体域、流体域的网格是不匹配的，通常流体域的网格在耦合面上比固体域要密得多，因此需要将流体载荷映射到固体节点上。在某个耦合单元面内，有n个流体结点，则作用于固体结点的力表示为：

[0054]

[0055] 其中s代表固体，f代表流体，Fsi代表作用于某固体单元第i个结点上的等效流体载荷， 表示第k个流体单元作用力映射到第i个固体结点对应的等参插值系数； 表示第k个流体单元对当前固体单元的作用力。

[0056] 等参插值系数，需要计算流体结点在固体单元中的参数坐标。可通过数值方法如牛顿迭代法计算。位移映射与流体载荷映射类似，可表达为：

[0057]

[0058] 式中，ufj代表某耦合单元上的第j个流体结点的位移， 表示第i个固体单元的位移映射到第j个流体结点对应的等参插值系数，usi表示第i个固体结点的位移。

[0059] 接下来，结合具体实施方式对陶瓷基复合材料流固耦合响应计算方法进行具体描述，该方法包括以下步骤：

[0060] S1:建立编织CMCs代表体元的有限元模型，并赋予纤维束合适的细观力学模型。

[0061] S2:给予代表体元有限元模型一系列加卸载路径，加卸载过程中最大应变值应逐渐增大，获得与不同最大应变值相对应的迟滞回线，获得迟滞回线越多，计算结果将越精确。注意卸载时要卸载到裂纹闭合，即继续卸载，应力应变关系为线性。

[0062] S3：由于每个迟滞回线的凹凸性通常来说最多变化一次，因此利用三次多项式对这些迟滞回线进行拟合可获得足够的精度，拟合得到不同εi对应的多项式系数an,bn(n＝1～4)：

[0063]

[0064] S4:则任意最大应变值εt，当其处于某两个RVE模型计算所得的最大应变值之间时，即(εi＜εt＜εi+1)时，则当前迟滞回线的多项式系数可插值为：

[0065]

[0066] S5：当振幅由大到小时，加卸载发生在最大迟滞回线内部时，已知当前的应力应变，则下一时刻的应力可由下式求得：

[0067]

[0068] S6：进行固体单元与流体单元之间的配对，通常流体网格相较于固体网格要密得多，则每个固体单元对应n个流体结点。然后通过下式，实现流体载荷向固体结点的映射。其中的等参插值系数，已知流体面心坐标和4个固体结点坐标，可通过牛顿迭代法计算得到：

[0069]

[0070] S7：基于上述任意加卸载应力应变计算方法，结合显式动力学积分，其中，固体结点载荷为S6计算所得，即可求得当前时间步的CMC结构的流固耦合动力学响应。

[0071] S8：根据S7计算的固体结点位移结果，将固体结点的位移映射到耦合面的流体结点。固体结点位移映射，也要匹配每个固体单元的某个面和面内的n个流体结点。固体结点的位移通过下式，实现向流体结点的映射。其中的等参插值系数，与上一步类似，通过牛顿迭代法计算。

[0072]

[0073] S9：根据S8计算的耦合面流体结点位移结果对流场结点位置进行更新，进行下一个耦合步的计算，即转至S6。

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

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