CN104866690A - 单向陶瓷基复合材料任意加卸载应力应变行为预测方法 - Google Patents

Abstract

本发明涉及一种复合材料应力应变行为预测方法，特别是一种单向陶瓷基复合材料任意加卸载应力应变行为预测方法。本发明的目的在于克服现有技术的缺陷而提供一种能快速预测单向陶瓷基复合材料任意加卸载过程中应力应变行为的方法。本发明提供了一种单向陶瓷基复合材料任意加卸载应力应变行为预测方法，考虑了基体开裂、纤维断裂、界面滑移和界面磨损等失效机理。提出了正反向滑移区的产生和覆盖规律，还给出了存在任意多个正反向滑移区时的应力分布和应变。本发明给出的公式大部分都存在解析解，因此能快速预测出单向陶瓷基复合材料任意加卸载下的应力应变行为。

Description

单向陶瓷基复合材料任意加卸载应力应变行为预测方法

技术领域

本发明涉及一种复合材料应力应变行为预测方法，特别是一种单向陶瓷基复合材料任意加卸载应力应变行为预测方法。

背景技术

陶瓷基复合材料具有高比强度、高比刚度、耐高温和低密度等优点，在航空发动机燃烧室和尾喷管调节片等部件上有着广泛的应用前景。在航空发动机工作时，燃烧室和尾喷管调节片由于不稳定燃烧和高速气流的作用，会受到随机激励而产生随机振动和疲劳。在随机激励下，材料受到的是不规律的变幅循环载荷。材料任意加卸载(即变幅循环载荷)下的应力应变行为决定了部件的随机振动响应和疲劳行为。因此，研究陶瓷基复合材料在任意加卸载下的应力应变行为对其应用有着重要意义。

由于陶瓷基复合材料是一种新型结构材料，国内外还没有高效的方法预测其任意加卸载下的应力应变行为，也未见公开的发明专利。Solti等(Solti JP.Modeling of progressivedamage in fiber-reinforced ceramic matrix composites.DTIC Document；1996.)和李龙彪(李龙彪,长纤维增强陶瓷基复合材料疲劳损伤模型与寿命预测,2010,南京航空航天大学.第184页.)模拟了单向陶瓷基复合材料在等幅循环载荷下的应力应变响应。他们的计算结果与实验结果较为一致,但是无法给出任意加卸载下的应力应变响应。方光武等(方光武,高希光,宋迎东.单向纤维增强陶瓷基复合材料界面滑移规律.复合材料学报.2013；4:101-107.)和高希光等(Xiguang G,Guangwu F,Yingdong S.Hysteresis loop modelof unidirectional carbon fiber-reinforced ceramic matrix composites under anarbitrary cyclic load.Compos,B,Eng.2014；56:92-99.)基于纤维/基体间位移增量平衡原理发展了单向陶瓷基复合材料的界面摩擦模型，预测了循环加载下的应力应变行为。但是由于存在大量的数值计算，计算耗时长，这不是一种高效的方法。当前，如何快速预测单向陶瓷基复合材料任意加卸载下的应力应变行为是本技术领域重要而难以解决的技术问题。

发明内容

1、所要解决的技术问题：

本发明的目的在于克服上述现有技术的缺陷而提供一种能快速预测单向陶瓷基复合材料任意加卸载过程中应力应变行为的方法。

2、技术方案：

为了解决以上问题，本发明提供了一种单向陶瓷基复合材料任意加卸载应力应变行为预测方法，包括如下步骤：

1)判断材料是否产生基体裂纹，若未产生基体裂纹，则根据混合率公式计算应变，否则执行步骤2)；

2)判断滑移区个数；

3)基于载荷历程计算当前界面剪应力；

4)基于Weibull失效分布计算纤维断裂率D，计算新的纤维体积分数；

5)基于临界基体应变能准则计算基体裂纹间距L；

6)计算新滑移区的长度；

7)判断是否发生滑移区覆盖，若新滑移区长度大于原滑移区长度，则发生覆盖，滑移区个数n减1，并重新执行步骤4)，否则执行步骤8)；

8)计算各滑移区的应力分布；

9)计算并输出应变。

所述步骤1)中，当此时和之前的载荷应力都小于基体临界开裂应力σ_cr，则判断基体未产生裂纹。混合率公式是：

$\mathrm{\ϵ}=\frac{\mathrm{\σ}}{v_f{E}_f+v_m{E}_m}.$

所述步骤2)中，当载荷由增大变为减小时，则反向滑移区个数加1；当载荷由减小变为增大时，则正向滑移区个数加1，否则滑移区个数不变。

所述步骤3)中，τ_i计算公式为：

$({\mathrm{\τ}}_i-{\mathrm{\τ}}_{i_0})/({\mathrm{\τ}}_{i_{\min }}-{\mathrm{\τ}}_{i_0})=1-\exp (-\mathrm{\ω}{((\∫|\mathrm{d\σ}\left|\right)/{\mathrm{\σ}}_A)}^{\mathrm{\λ}}).$

所述步骤4)中，Weibull失效分布为：

D＝exp(-|σ_max/σ₀|^-m)。

新的纤维体积分数为：

$v_f={\mathrm{Dv}}_{f_0}.$

所述步骤5)中，临界基体应变能准则为：

$U_m=U_{{\mathrm{cr}}_m},$

其中，U_m＝∫_V∫_εσ_m(x)dεdV。

所述步骤6)中，当载荷变化时，只有新滑移区的长度变化，之前的各滑移区长度不变。

当只有1个正向滑移区时，新滑移区长度为：

$d_{f_1}=\frac{1}{2}{r}_f\frac{v_m{\mathrm{\σ}}_{m_0}}{v_f{\mathrm{\τ}}_i}-\frac{1}{\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i};$

当有1个正向滑移区和1个反向滑移区时，新滑移区长度为：

$d_{r_1}=\frac{1}{2}{d}_{f_1}-(\frac{v_m{r}_f}{{4v}_f{\mathrm{\τ}}_i}{\mathrm{\σ}}_{m_0}-\frac{1}{2\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i});$

当有n(n≥2)个正向滑移区和n-1个反向滑移区时，新滑移区长度为：

$d_{f_n}=\frac{1}{2}{d}_{f_1}+\frac{v_m{r}_f}{{4v}_f{\mathrm{\τ}}_i}{\mathrm{\σ}}_{m_0}-\frac{1}{2\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i}+\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}(d_{r_i}-d_{f_i}),n\≥2;$

当有n(n≥2)个正向滑移区和n个反向滑移区时，新滑移区长度为：

$d_{r_n}=\frac{1}{2}{d}_{f_1}-\frac{v_m{r}_f}{{4v}_f{\mathrm{\τ}}_i}{\mathrm{\σ}}_{m_0}+\frac{1}{2\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i}+\underset{k=2}{\overset{n}{\mathrm{\Σ}}}(d_{f_k}-d_{r_{k-1}}),n\≥2.$

所述步骤8)中，当只有1个正向滑移区时，各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}-\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x),& \frac{L}{2}-d_{f_1}\≤x\≤\frac{L}{2}\\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}-2\frac{{\mathrm{\τ}}_i}{r_f}{d}_{f_1}),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right.;$

当有1个正向滑移区和1个反向滑移区时，各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}+\frac{2{\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x),& \frac{L}{2}-d_{r_1}\≤x\≤\frac{L}{2}\\ \frac{\mathrm{\σ}}{v_f}-\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x-{2d}_{r_1}),& \frac{L}{2}-d_{f_1}\≤x\≤\frac{L}{2}-d_{r_1}\\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\×(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}-2\frac{{\mathrm{\τ}}_i}{r_f}(d_{f_1}-{2d}_{r_1})),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right.;$

当有n个(n≥2)正向滑移区和n-1个反向滑移区时，各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(x-\frac{L}{2}),& \frac{L}{2}-d_{f_n}\≤x\≤\frac{L}{2}\\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x-{2d}_{f_n}),& \frac{L}{2}-d_{r_{n-1}}\≤x\≤\frac{L}{2}-d_{f_n}\\ .& \\ .& \\ .& \\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(x-\frac{L}{2}+2\underset{k=j}{\overset{n-1}{\mathrm{\Σ}}}(d_{r_k}-d_{f_{k+1}})),& \frac{L}{2}-d_{f_j}\≤x\≤\frac{L}{2}-d_{r_j}\\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x+2\underset{k=j}{\overset{n-1}{\mathrm{\Σ}}}(d_{r_k}-d_{f_{k+1}})-{2d}_{f_j}),& \frac{L}{2}-d_{r_{j-1}}\≤x\≤\frac{L}{2}-d_{f_j}\\ .& \\ .& \\ .& \\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}& \\ +2\frac{{\mathrm{\τ}}_i}{r_f}(2\underset{k=1}{\overset{n-1}{\mathrm{\Σ}}}(d_{r_k}-d_{f_{k=1}})+d_{f_1}),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right.;$

当有n(n≥2)个正向滑移区和n个反向滑移区时，各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x),& \frac{L}{2}-d_{r_n}\≤x\≤\frac{L}{2}\\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(x-\frac{L}{2}+{2d}_{r_n}),& \frac{L}{2}-d_{f_n}\≤x\≤\frac{L}{2}-d_{r_n}\\ .& \\ .& \\ .& \\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x+2\underset{k=j+1}{\overset{n}{\mathrm{\Σ}}}(d_{r_k}-d_{f_k})),& \frac{L}{2}-d_{r_j}\≤x\≤\frac{L}{2}-d_{f_{j+1}}\\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(x-\frac{L}{2}+2\underset{k=j+1}{\overset{n}{\mathrm{\Σ}}}(d_{r_k}-d_{f_k})+{2d}_{r_j}),& \frac{L}{2}-d_{f_j}\≤x\≤\frac{L}{2}-d_{r_j}\\ .& \\ .& \\ .& \\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}& \\ \×(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}+2\frac{{\mathrm{\τ}}_i}{r_f}(2\underset{k=1}{\overset{n}{\mathrm{\Σ}}}(d_{r_k}-d_{f_k})+d_{f_1})),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right..$

所述步骤9)中应变的计算公式是：

${\stackrel{\‾}{\mathrm{\ϵ}}}_c={\stackrel{\‾}{\mathrm{\ϵ}}}_f=\frac{2}{E_{f}L}{\∫}_0^{L/2}{\mathrm{\σ}}_f\left(x\right)\mathrm{dx}+({\mathrm{\α}}_f-{\mathrm{\α}}_c)\mathrm{\ΔT}.$

3、有益效果：

本发明提供的一种单向陶瓷基复合材料任意加卸载应力应变行为预测方法，考虑了基体开裂、纤维断裂、界面滑移和界面磨损等失效机理。提出了正反向滑移区的产生和覆盖规律，还给出了存在任意多个正反向滑移区时的应力分布和应变。本发明给出的公式大部分都存在解析解，因此能快速预测出单向陶瓷基复合材料任意加卸载下的应力应变行为。

附图说明

图1是产生基体裂纹后的材料代表体元。

图2是滑移区覆盖过程。

图3是滑移区长度变化过程。

图4是代表体元含有1个粘结区、1个正向滑移区和1个反向滑移区的示意图。

图5是代表体元含有1个粘结区、2个正向滑移区和1个反向滑移区的示意图。

具体实施方式

本发明一种单向陶瓷基复合材料任意加卸载过程下应力应变行为预测方法，该方法具体步骤如下：

1)判断材料是否产生基体裂纹，若未产生基体裂纹，则根据混合率公式计算应变，否则执行步骤2)，产生基体裂纹后的材料代表体元，如图1所示；

2)判断滑移区个数；

3)基于载荷历程计算当前界面剪应力τ_i；

4)基于Weibull失效分布计算当前纤维断裂率D，并计算新的纤维体积分数v_f；

5)基于临界基体应变能准则计算基体裂纹间距L；

6)计算新滑移区的长度；

7)判断是否发生滑移区覆盖，若新滑移区长度大于原滑移区长度，则发生覆盖，覆盖过程如图2所示，滑移区个数n减1，并重新执行步骤4)，否则执行步骤8)；

8)计算各滑移区的应力分布；

9)计算并输出应变。

所述步骤1)中当此时和之前的载荷应力都小于基体临界开裂应力σ_cr，则判断基体未产生裂纹。混合率公式是：

$\mathrm{\ϵ}=\frac{\mathrm{\σ}}{v_f{E}_f+v_m{E}_m}.$

所述步骤2)中，当载荷由增大变为减小时，则反向滑移区个数加1，当载荷由减小变为增大时，则正向滑移区个数加1，否则滑移区个数不变。

所述步骤3)中，τ_i计算公式为：

$({\mathrm{\τ}}_i-{\mathrm{\τ}}_{i_0})/({\mathrm{\τ}}_{i_{\min }}-{\mathrm{\τ}}_{i_0})=1-\exp (-\mathrm{\ω}{((\∫|\mathrm{d\σ}\left|\right)/{\mathrm{\σ}}_A)}^{\mathrm{\λ}})$

所述步骤4)中，Weibull失效分布为：

D＝exp(-|σ_max/σ₀|^-m)

新的纤维体积分数为：

$v_f={\mathrm{Dv}}_{f_0}$

所述步骤5)中，临界基体应变能准则为：

$U_m=U_{{\mathrm{cr}}_m}$

其中，U_m＝∫_V∫_εσ_m(x)dεdV。

所述步骤6)中，当载荷变化时，只有新滑移区的长度变化，之前的各滑移区长度不变，滑移区长度变化过程如图3所示。

当只有1个正向滑移区时，新滑移区长度为

$d_{f_1}=\frac{1}{2}{r}_f\frac{v_m{\mathrm{\σ}}_{m_0}}{v_f{\mathrm{\τ}}_i}-\frac{1}{\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i};$

当有1个正向滑移区和1个反向滑移区时，新滑移区长度为：

$d_{r_1}=\frac{1}{2}{d}_{f_1}-(\frac{v_m{r}_f}{{4v}_f{\mathrm{\τ}}_i}{\mathrm{\σ}}_{m_0}-\frac{1}{2\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i});$

当有n(n≥2)个正向滑移区和n-1个反向滑移区时，新滑移区长度为：

$d_{f_n}=\frac{1}{2}{d}_{f_1}+\frac{v_m{r}_f}{{4v}_f{\mathrm{\τ}}_i}{\mathrm{\σ}}_{m_0}-\frac{1}{2\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i}+\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}(d_{r_i}-d_{f_i}),n\≥2;$

当有n(n≥2)个正向滑移区和n个反向滑移区时，新滑移区长度为：

$d_{r_n}=\frac{1}{2}{d}_{f_1}-\frac{v_m{r}_f}{{4v}_f{\mathrm{\τ}}_i}{\mathrm{\σ}}_{m_0}+\frac{1}{2\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i}+\underset{k=2}{\overset{n}{\mathrm{\Σ}}}(d_{f_k}-d_{r_{k-1}}),n\≥2.$

所述步骤8)中，当只有1个正向滑移区时，各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}-\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x),& \frac{L}{2}-d_{f_1}\≤x\≤\frac{L}{2}\\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}-2\frac{{\mathrm{\τ}}_i}{r_f}{d}_{f_1}),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right.;$

当有1个正向滑移区和1个反向滑移区时，各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}+\frac{2{\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x),& \frac{L}{2}-d_{r_1}\≤x\≤\frac{L}{2}\\ \frac{\mathrm{\σ}}{v_f}-\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x-{2d}_{r_1}),& \frac{L}{2}-d_{f_1}\≤x\≤\frac{L}{2}-d_{r_1}\\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\×(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}-2\frac{{\mathrm{\τ}}_i}{r_f}(d_{f_1}-{2d}_{r_1})),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right.;$

当有n个(n≥2)正向滑移区和n-1个反向滑移区时，各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(x-\frac{L}{2}),& \frac{L}{2}-d_{f_n}\≤x\≤\frac{L}{2}\\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x-{2d}_{f_n}),& \frac{L}{2}-d_{r_{n-1}}\≤x\≤\frac{L}{2}-d_{f_n}\\ .& \\ .& \\ .& \\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(x-\frac{L}{2}+2\underset{k=j}{\overset{n-1}{\mathrm{\Σ}}}(d_{r_k}-d_{f_{k+1}})),& \frac{L}{2}-d_{f_j}\≤x\≤\frac{L}{2}-d_{r_j}\\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x+2\underset{k=j}{\overset{n-1}{\mathrm{\Σ}}}(d_{r_k}-d_{f_{k+1}})-{2d}_{f_j}),& \frac{L}{2}-d_{r_{j-1}}\≤x\≤\frac{L}{2}-d_{f_j}\\ .& \\ .& \\ .& \\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}& \\ +2\frac{{\mathrm{\τ}}_i}{r_f}(2\underset{k=1}{\overset{n-1}{\mathrm{\Σ}}}(d_{r_k}-d_{f_{k=1}})+d_{f_1}),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right.;$

当有n(n≥2)个正向滑移区和n个反向滑移区时，各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x),& \frac{L}{2}-d_{r_n}\≤x\≤\frac{L}{2}\\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(x-\frac{L}{2}+{2d}_{r_n}),& \frac{L}{2}-d_{f_n}\≤x\≤\frac{L}{2}-d_{r_n}\\ .& \\ .& \\ .& \\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x+2\underset{k=j+1}{\overset{n}{\mathrm{\Σ}}}(d_{r_k}-d_{f_k})),& \frac{L}{2}-d_{r_j}\≤x\≤\frac{L}{2}-d_{f_{j+1}}\\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(x-\frac{L}{2}+2\underset{k=j+1}{\overset{n}{\mathrm{\Σ}}}(d_{r_k}-d_{f_k})+{2d}_{r_j}),& \frac{L}{2}-d_{f_j}\≤x\≤\frac{L}{2}-d_{r_j}\\ .& \\ .& \\ .& \\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}& \\ \×(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}+2\frac{{\mathrm{\τ}}_i}{r_f}(2\underset{k=1}{\overset{n}{\mathrm{\Σ}}}(d_{r_k}-d_{f_k})+d_{f_1})),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right..$

所述步骤9)中应变的计算公式是：

${\stackrel{\‾}{\mathrm{\ϵ}}}_c={\stackrel{\‾}{\mathrm{\ϵ}}}_f=\frac{2}{E_{f}L}{\∫}_0^{L/2}{\mathrm{\σ}}_f\left(x\right)\mathrm{dx}+({\mathrm{\α}}_f-{\mathrm{\α}}_c)\mathrm{\ΔT}.$

下面结合附图和具体实施例对本发明进行详细说明。

实施例：经过N个载荷循环，第N个循环卸载结束后材料代表体元有1个粘结区、1个正向滑移区和1个反向滑移区，如图4所示。当前开始第N+1个循环的加载。

1)判断材料是否产生基体裂纹。当前状态下，材料已发生损伤，产生基体裂纹。

2)判断滑移区个数。当前状态是载荷由减小变为增大，所以正向滑移区个数加1，有2个正向滑移区，1个反向滑移区，如图5所示。

3)基于载荷历程计算当前界面剪应力τ_i。本实施例中，当前状态下经过的载荷历程为：

$\∫\left|\mathrm{d\σ}\right|={\mathrm{\σ}}_1^{\max }+|{\mathrm{\σ}}_1^{\min }-{\mathrm{\σ}}_1^{\max }|+\underset{n=2}{\overset{N}{\mathrm{\Σ}}}\left(\right|{\mathrm{\σ}}_n^{\max }-{\mathrm{\σ}}_{n-1}^{\min }|+|{\mathrm{\σ}}_n^{\min }-{\mathrm{\σ}}_n^{\max }\left|\right).$

其中和分别是第n个载荷循环的峰值和谷值。则当前状态下界面剪应力τ_i可由下式计算得到：

$\begin{array}{c}({\mathrm{\τ}}_i-{\mathrm{\τ}}_{i_0})/({\mathrm{\τ}}_{i_{\min }}-{\mathrm{\τ}}_{i_0})\\ =1-\exp \{-\mathrm{\ω}{[({\mathrm{\σ}}_1^{\max }+|{\mathrm{\σ}}_1^{\min }-{\mathrm{\σ}}_1^{\max }|+\underset{n=2}{\overset{N}{\mathrm{\Σ}}}\left(\right|{\mathrm{\σ}}_n^{\max }-{\mathrm{\σ}}_{n-1}^{\min }|+|{\mathrm{\σ}}_n^{\min }-{\mathrm{\σ}}_n^{\max }\left|\right))/{\mathrm{\σ}}_A]}^{\mathrm{\λ}}\}\end{array}.$

4)基于Weibull失效分布计算当前纤维断裂率D，计算纤维断裂后的新体积分数v_f。当前纤维断裂率为：

D＝exp(-|σ_max/σ₀|^-m)。

新的纤维体积分数为：

$v_f={\mathrm{Dv}}_{f_0}=\exp (-|{\mathrm{\σ}}_{\max }/{\mathrm{\σ}}_0{|}^{-m})v_{f_0}.$

5)基于临界基体应变能准则计算基体裂纹间距L。当前状态下的基体应变能为：

$\begin{array}{c}{U}_m=\frac{A_m}{E_m}{\∫}_0^{L/2}{\mathrm{\σ}}_m^2\left(x\right)\mathrm{dx}\\ =\frac{A_m}{E_m}{\∫}_0^{L/2}{(\frac{\mathrm{\σ}}{v_m}-\frac{v_f}{v_m}{\mathrm{\σ}}_f)}^2\mathrm{dx}\\ =\frac{A_m{\mathrm{L\σ}}^2}{{2E}_m{v}_m^2}-\frac{{2A}_m{v}_f\mathrm{\σ}}{E_m{v}_m^2}{\∫}_0^{L/2}{\mathrm{\σ}}_f\left(x\right)\mathrm{dx}+\frac{A_m{v}_f^2}{E_m{v}_m^2}{\∫}_0^{L/2}\left(x\right)\mathrm{dx}\\ =\frac{A_m{\mathrm{L\σ}}^2}{2E_m{v}_m^2}-\frac{{2A}_m{v}_f\mathrm{\σ}}{E_m{v}_m^2}\{\frac{\mathrm{\σ}}{v_f}{d}_{f_1}+{\mathrm{\σ}}_{f_0}(\frac{L}{2}-d_{f_1})+\frac{{\mathrm{\τ}}_i}{r_f}(-d_{f_1}^2-{2d}_{f_1}^2+{4d}_{f_1}{d}_{r_1}+{2d}_{f_2}^2-{4d}_{f_1}{d}_{f_2})\\ +\frac{1}{\mathrm{\β}}\tanh \left[\mathrm{\β}(\frac{L}{2}-d_{f_1})\right][\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}+\frac{{2\mathrm{\τ}}_i}{r_f}(-d_{f_1}+{2d}_{r_1}-{2d}_{f_2})]\}\\ +\frac{A_m{v}_f^2}{E_m{v}_m^2}\{{\mathrm{\σ}}_{f_0}^2(\frac{L}{2}-d_{f_1})+\frac{A_2^2}{4\mathrm{\β}}\sinh \mathrm{\β}(L-{2d}_{f_1})+\frac{A_2^2}{2}(\frac{L}{2}-d_{f_1})+\frac{{2\mathrm{\σ}}_{f_0}{A}_2}{\mathrm{\β}}\sinh \mathrm{\β}(\frac{L}{2}-d_{f_1})\\ +\frac{r_f}{{6\mathrm{\τ}}_i}[2{(\frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(d_{r_1}-{2d}_{f_2}))}^3+{\left(\frac{\mathrm{\σ}}{v_f}\right)}^3-{(\frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(-d_{f_1}+{2d}_{f_1}-{2d}_{f_2}))}^3-2{(\frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(-d_{f_2}))}^3]\},\end{array}$

L可由计算得出。

其中：

$A_2=\frac{(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}+2\frac{{\mathrm{\τ}}_i}{r_f}(2(d_{r_1}-d_{f_2})+d_{f_1}))}{\cosh \left[(L/2-d_{f_1})\right]}$

${\mathrm{\σ}}_{f_0}=\frac{E_f}{E_1}\mathrm{\σ}+E_f({\mathrm{\α}}_1-{\mathrm{\α}}_f)\mathrm{\ΔT}$

${\mathrm{\σ}}_{m_0}=\frac{E_m}{E_1}\mathrm{\σ}+E_m({\mathrm{\α}}_1-{\mathrm{\α}}_m)\mathrm{\ΔT}$

${\mathrm{\β}}^2=\frac{{8E}_1}{v_m{E}_f{E}_m{r}_f^2}{[\frac{1}{G_f}+\frac{1}{G_m}(\frac{2}{v_m^2}\ln \frac{1}{v_f}-3-2\frac{v_f}{v_m})]}^{-1}.$

6)计算新滑移区的长度。只有新滑移区的长度会变化，原有的1个正向滑移区长度和1个反向滑移区长度保持不变。新滑移区长度为：

$d_{f_2}=\frac{1}{2}{d}_{f_1}+\frac{v_m{r}_f}{{4v}_f{\mathrm{\τ}}_i}{\mathrm{\σ}}_{m_0}-\frac{1}{2\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i}+d_{r_1}-d_{f_1}.$

7)判断是否发生滑移区覆盖。当前状态下，载荷由减小变为增大，所以载荷较小，，不会发生滑移区覆盖。

8)计算各滑移区的应力分布。各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(x-\frac{L}{2}),& \frac{L}{2}-d_{f_2}\≤x\≤\frac{L}{2}\\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x-{2d}_{f_2}),& \frac{L}{2}-d_{r_1}\≤x\≤\frac{L}{2}-d_{f_2}\\ \frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(x-\frac{L}{2}+2(d_{r_1}-d_{f_2})),& \frac{L}{2}-d_{f_1}\≤x\≤\frac{L}{2}-d_{r_1}\\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}+2\frac{{\mathrm{\τ}}_i}{r_f}(2(d_{r_1}-d_{f_2})+d_{f_1})),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right..$

其中表示粘结区，表示第1个正向滑移区，表示第1个反向滑移区，表示第2个正向滑移区。

9)计算并输出应变。当前状态下应变为：

$\begin{array}{c}{\stackrel{\‾}{\mathrm{\ϵ}}}_c={\stackrel{\‾}{\mathrm{\ϵ}}}_f=\frac{2}{E_{f}L}{\∫}_0^{L/2}{\mathrm{\σ}}_f\left(x\right)\mathrm{dx}+({\mathrm{\α}}_f-{\mathrm{\α}}_c)\mathrm{\ΔT}\\ =\frac{2}{E_{f}L}\{\frac{\mathrm{\σ}}{v_f}{d}_{f_1}+{\mathrm{\σ}}_{f_0}(\frac{L}{2}-d_{f_1})+\frac{{\mathrm{\τ}}_i}{r_f}(-d_{f_1}^2-{2d}_{r_1}^2+{4d}_{f_1}{d}_{r_1}+{2d}_{f_2}^2-{4d}_{f_1}{d}_{f_2})\\ +\frac{1}{\mathrm{\β}}\tanh \left[\mathrm{\β}(\frac{L}{2}-d_{f_1})\right][\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}+\frac{{2\mathrm{\τ}}_i}{r_f}(-d_{f_1}+{2d}_{r_1}-{2d}_{f_2})]\}+({\mathrm{\α}}_f-{\mathrm{\α}}_c)\mathrm{\ΔT}\end{array}.$

虽然本发明已以较佳实施例公开如上，但它们并不是用来限定本发明的，任何熟习此技艺者，在不脱离本发明之精神和范围内，自当可作各种变化或润饰，因此本发明的保护范围应当以本申请的权利要求保护范围所界定的为准。

Claims (9)

1.一种单向陶瓷基复合材料任意加卸载应力应变行为预测方法，其特征在于：包括如下步骤：

1)判断材料是否产生基体裂纹，若未产生基体裂纹，则根据混合率公式计算应变，否则执行步骤2)；

2)判断滑移区个数；

3)基于载荷历程计算当前界面剪应力；

4)基于Weibull失效分布计算纤维断裂率D，计算新的纤维体积分数；

5)基于临界基体应变能准则计算基体裂纹间距L；

6)计算新滑移区的长度；

7)判断是否发生滑移区覆盖，若新滑移区长度大于原滑移区长度，则发生覆盖，滑移区个数n减1，并重新执行步骤4)，否则执行步骤8)；

8)计算各滑移区的应力分布；

9)计算并输出应变。

2.如权利要求1所述的单向陶瓷基复合材料任意加卸载应力应变行为预测方法，其特征在于：所述步骤1)中，当此时和之前的载荷应力都小于基体临界开裂应力，则判断基体未产生裂纹，混合率公式是：

3.如权利要求1所述的单向陶瓷基复合材料任意加卸载应力应变行为预测方法，其特征在于：所述步骤2)中，当载荷由增大变为减小时，则反向滑移区个数加1；当载荷由减小变为增大时，则正向滑移区个数加1，否则滑移区个数不变。

4.如权利要求1所述的单向陶瓷基复合材料任意加卸载应力应变行为预测方法，其特征在于：所述步骤3)中，计算公式为： $({\mathrm{\τ}}_i-{\mathrm{\τ}}_{i_0})/({\mathrm{\τ}}_{i_{\min }}-{\mathrm{\τ}}_{i_0})=1-\exp (-\mathrm{\ω}{((\∫|\mathrm{d\σ}\left|\right)/{\mathrm{\σ}}_A)}^{\mathrm{\λ}}).$

5.如权利要求1所述的单向陶瓷基复合材料任意加卸载应力应变行为预测方法，其特征在于：所述步骤4)中，Weibull失效分布为：D＝exp(-|σ_max/σ₀|^-m)，新的纤维体积分数为： $v_f={\mathrm{Dv}}_{f_0}.$

6.如权利要求1所述的单向陶瓷基复合材料任意加卸载应力应变行为预测方法，其特征在于：所述步骤5)中，临界基体应变能准则为：其中，U_m＝∫_V∫_εσ_m(x)dεdV。

7.如权利要求1所述的单向陶瓷基复合材料任意加卸载应力应变行为预测方法，其特征在于：所述步骤6)中，当载荷变化时，只有新滑移区的长度变化，之前的各滑移区长度不变。

当只有1个正向滑移区时，新滑移区长度为：

$d_{f_1}=\frac{1}{2}{r}_f\frac{v_m{\mathrm{\σ}}_{m_0}}{v_f{\mathrm{\τ}}_i}-\frac{1}{\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i};$

当有1个正向滑移区和1个反向滑移区时，新滑移区长度为：

$d_{r_1}=\frac{1}{2}{d}_{f_1}-(\frac{v_m{r}_f}{4v_f{\mathrm{\τ}}_i}{\mathrm{\σ}}_{m_0}-\frac{1}{2\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i})\text{;}$

当有n(n≥2)个正向滑移区和n-1个反向滑移区时，新滑移区长度为：

$d_{f_n}=\frac{1}{2}{d}_{f_1}+\frac{v_m{r}_f}{4v_f{\mathrm{\τ}}_i}{\mathrm{\σ}}_{m_0}-\frac{1}{2\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i}+\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}(d_{r_i}-d_{f_i}),n\≥2;$

当有n(n≥2)个正向滑移区和n个反向滑移区时，新滑移区长度为：

$d_{r_n}=\frac{1}{2}{d}_{f_1}-\frac{v_m{r}_f}{4v_f{\mathrm{\τ}}_i}{\mathrm{\σ}}_{m_0}+\frac{1}{2\mathrm{\β}\tanh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\frac{{\mathrm{\τ}}_{\mathrm{ult}}}{{\mathrm{\τ}}_i}+\underset{k=2}{\overset{n-1}{\mathrm{\Σ}}}(d_{f_k}-d_{r_{k-1}}),n\≥2.$

8.权利要求1所述的单向陶瓷基复合材料任意加卸载应力应变行为预测方法，其特征在于：所述步骤8)中，当只有1个正向滑移区时，各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}-\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x),& \frac{L}{2}-d_{f_1}\≤x\≤\frac{L}{2}\\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}-2\frac{{\mathrm{\τ}}_i}{r_f}{d}_{f_1}),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right.;$

当有1个正向滑移区和1个反向滑移区时，各区域应力分布为：

${\mathrm{\σ}}_f\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{\σ}}{v_f}+\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x),& \frac{L}{2}-d_{r_1}\≤x\≤\frac{L}{2}\\ \frac{\mathrm{\σ}}{v_f}-\frac{{2\mathrm{\τ}}_i}{r_f}(\frac{L}{2}-x-{2d}_{r_1}),& \frac{L}{2}-d_{f_1}\≤x\≤\frac{L}{2}-d_{r_1}\\ {\mathrm{\σ}}_{f_0}+\frac{\cosh \left(\mathrm{\βx}\right)}{\cosh \left[\mathrm{\β}(L/2-d_{f_1})\right]}\×(\frac{v_m}{v_f}{\mathrm{\σ}}_{m_0}-2\frac{{\mathrm{\τ}}_i}{r_f}(d_{f_1}-{2d}_{r_1})),& 0\≤x\≤\frac{L}{2}-d_{f_1}\end{array}\right.;$

当有n个(n≥2)正向滑移区和n-1个反向滑移区时，各区域应力分布为：

当有n(n≥2)个正向滑移区和n个反向滑移区时，各区域应力分布为：

9.权利要求1所述的单向陶瓷基复合材料任意加卸载应力应变行为预测方法，其特征在于：所述步骤9)中应变的计算公式是： ${\stackrel{\‾}{\mathrm{\ϵ}}}_c={\stackrel{\‾}{\mathrm{\ϵ}}}_f=\frac{2}{E_{f}L}{\∫}_0^{L/2}{\mathrm{\σ}}_f\left(x\right)\mathrm{dx}+({\mathrm{\α}}_f-{\mathrm{\α}}_c)\mathrm{\ΔT}.$