一种复合材料层合结构的非概率动力可靠性评估方法

摘要

本发明公开了一种面向全寿命周期的复合材料层合结构非概率动力可靠性评估方法。该方法综合考虑复材层合结构载荷的动态波动效应以及材料性能的累积退化，基于有限元思想和区间数学方法，构建了该结构动特性的非概率区间过程模型；进而结合首次穿越理论和轻质层合结构的失效判定准则，定义了结构非概率动力可靠性指标，并探索了高效稳健的求解策略。本发明在进行可靠度计算过程中合理表征了不确定性对结构全生命周期内动力安全的综合影响，为确保对其开展精细化设计提供必要的数据参考。

说明书

技术领域

本发明涉及轻质复材层合结构的安全性评估技术领域，特别涉及考虑不确定性、载荷交变及材料退化共同作用下结构安全性能的有效认知与定量表征，为合理制定复合材料层合结构的选材及铺层方案设计提供重要的理论支持。

背景技术

由于质量轻、强度高、多功能，使复合材料层合板成为具有优良特性的工程结构件。此外，鉴于其自身具有结构设计裕度大、性能匹配度高等天然优势，复材层合结构不仅被大量应用于航空航天、船舶、兵器等军工领域，还作为最基本和最主要的构件频繁出现在民用客机、汽车、建筑等民用结构系统中。因此，针对复合材料层合结构的力学特性分析与安全态势评估技术研究具有重要的理论意义与工程实用价值。

然而，复合材料由于成型工艺以及加工批次的差异，材料缺陷导致的分散性不可避免。此外，结构的服役环境复杂多样，外部激励的未确知性同样客观存在。加之随着服役时间的累积，材料性能的退化与载荷的交变效应严重影响着复合材料层合结构的实际使用性能，如若无法有效准确地预估结构的真实服役状态，将导致重大的安全事故发生。综合上述情况，针对轻质复材层合板结构开展不确定性传播分析与动力可靠性评估方法研究已受到学术界和工程界的高度重视。

当前，国内外学者与工程技术人员对复合材料结构的可靠性分析方法研究主要集中在两个方面：(1)基于概率统计理论的结构不确定性传播影响预测技术；(2)基于准静态假设的可靠度计算方法研究。上述工作一定程度上丰富了复合材料结构的分析与强度理论，但是忽略了随机方法对样本信息的依赖性以及时间累积效应下结构失效事件的相关性，大大限制了其理论的工程实用化进程。

由于实际工程中贫信息、少数据的情况时有发生，建立以非概率理论框架为基础的不确定性表征技术、时变可靠度建模与求解技术具有显著的现实意义。目前，相关研究工作尚不成熟，针对复合材料层合结构的方案设计经常无法严格满足所需的应用要求，亦或是安全冗余度过大，造成严重的资源浪费与时间成本损耗。鉴于此，本发明将重点探究静动力不确定性作用下轻质层合结构安全性能的退化历程，为其轻量化、多功能设计提供理论保障。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供一种针对轻质复合材料层合板结构的安全性评价方法，充分考虑实际工程问题中普遍存在的静动力不确定性因素，以提出的非概率动力可靠性度量评判结构安全与否的量化指标，所得到的强度校核结果更加符合真实情况，工程适用性更强。

本发明采用的技术方案实现步骤如下：

第一步：根据复材层合结构的本构特征及载荷边界条件，基于最小势能原理构建有限元列式，推导出如下动力方程：其中，δ分别表示节点加速度、速度和位移，M,C,K分别为总体质量矩阵、总体阻尼矩阵和总体刚度矩阵，F为载荷列向量，N为节点个数。这里，复合材料层合板结构总体刚度矩阵的计算如下积分式表达：

$\left[K\right]={\int }_V{\left[B\right]}^T\left[D\right]\left[B\right]dV=\underset{k=1}{\overset{n}{\Sigma }}{\int }_{V_k}{\left[B\right]}^T\left[D\right]\left[B\right]{\mathrm{dV}}_k=\underset{k=1}{\overset{n}{\Sigma }}{\int }_{z_{k-1}}^{z_k}{\int }_A{\left[B\right]}^T\left[D\right]\left[B\right]dAdz$

其中，V代表全域体积，V_k表示第k层板的体积，n为层板总数，A为层板面积，B为应变矩阵。

第二步：引入强度比R，结合Tsai-Wu强度失效判定准则，代入到步骤一中的动力方程中，可得到临界失效载荷的显式表达式：

F_iσ_imax+F_ijσ_imaxσ_jmax＝1i,j＝1,2,…6

$[F_1{\sigma }_1+F_2{\sigma }_2]R+[F_{11}{\sigma }_1^2+F_{22}{\sigma }_2^2+F_{66}{\tau }_{12}^2+2F_{12}{\sigma }_1{\sigma }_2]R^2=1$

其中，F_i和F_ij表示强度特征量，σ＝[σ₁,σ₂,τ₁₂]^T为应力向量，σ_imax对应最危险点处的应力值，角标i和j分别表示单元对应的自由度计数指标。这里所述的强度比R其物理意义为极限载荷与真实载荷的比值，具体公式为：

$R=\frac{-Coe\_B+\sqrt{Coe\_B^2+4Coe\_A}}{2Coe\_A}$

其中，分别定义中间参量Coe_A＝F_ijσ_iσ_j和Coe_B＝F_iσ_i，σ_i和σ_j分别对应计数指标i和j处的单元应力。强度特征量F_i和F_ij取决于材料工程强度参数{X_t,X_c,Y_t,Y_c,S₂}，分别表示材料纵向拉伸强度、纵向压缩强度、横向拉伸强度、横向压缩强度以及剪切强度，即：

$F_1=\frac{1}{X_t}-\frac{1}{X_c},F_2=\frac{1}{Y_t}-\frac{1}{Y_c},F_{11}=\frac{1}{X_t{X}_c}$

$F_{22}=\frac{1}{Y_t{Y}_c},F_{66}=\frac{1}{S_2},F_{12}=-\frac{1}{2\sqrt{X_t{X}_c{Y}_t{Y}_c}}$

第三步：综合考虑存在于第一步和第二步中结构参数中的不确定性效应，不确定性效应包括材料分散性引起的M,C,K的变动、载荷F的未确知性引起响应δ的变化、强度判定准则的模糊性导致强度特征量F_i和F_ij的不一致。引入区间向量x∈x^I＝(x₁,x₂,...,x_m)以及区间过程模型X(t)∈X(t)^I＝(X₁(t),X₂(t),...,X_n(t))，得到有限样本信息条件下结构静动力不确定性参数的数学表达：

X(t)^U/L＝(X₁(t)^U/L,X₂(t)^U/L,...,X_n(t)^U/L)＝(X₁(t)^c±X₁(t)^r,X₂(t)^c±X₂(t)^r,...,X_m(t)^c±X_m(t)^r)其中，上标U代表参量的取值上界，上标L代表参量的取值下界，上标c代表中心值，上标r代表半径。静动力不确定性参数包括：载荷边界参数、材料特征参数、测量精度参数以及设计准则(设计许用值)参数等，为便于计算，通常表示为标准化形式，具体如下：

x＝[x^L,x^U]＝[x^c-x^r,x^c+x^r]

＝x^c+x^r[-1,1]

＝x^c+x^r×ξ₁

和

X(t)＝[X^L(t),X^U(t)]＝[X^c(t)-X^r(t),X^c(t)+X^r(t)]

＝X^c(t)+X^r(t)[-1,1]

＝X^c(t)+X^r(t)×ξ₂

其中，x^c＝(x₁^c,x₂^c,...,x_m^c)和x^r＝(x₁^r,x₂^r,...,x_m^r)表示区间向量x的均值和半径，X^c(t)＝(X₁^c(t),X₂^c(t),...,x_n^c(t))和X^c(t)＝(X₁^r(t),X₂^r(t),...,x_n^r(t))表示区间过程X(t)的均值和半径；ξ₁,ξ₂∈Ξ⁴，Ξ⁴定义为所有元素包含在[-1,1]内的4维向量集合，符号“×”定义为两个向量各对应元素相乘的算子，乘积仍为维数为4的向量。

第四步：将第三步中表征的静动力不确定参数代入到第一步和第二步的方程中，形成考虑时间效应的区间集格式，即给定任意时刻有：

$\begin{array}{c}{\left[M^I(x,X(t\left)\right)\right]}_{5N\times 5N}{\left\{{\stackrel{··}{\delta }}^I(x,X(t\left)\right)\right\}}_{5N}+{\left[C^I(x,X(t\left)\right)\right]}_{5N\times 5N}{\left\{{\stackrel{·}{\delta }}^I(x,X(t\left)\right)\right\}}_{5N}\\ +{\left[K^I(x,X(t\left)\right)\right]}_{5N\times 5N}{\left\{{\delta }^I(x,X(t\left)\right)\right\}}_{5N}={\left\{F^I(x,X(t\left)\right)\right\}}_{5N}\end{array}.$

其中，δ^I(x,X(t))分别为考虑不确定性影响下的加速度、速度和位移历程，M^I(x,X(t)),C^I(x,X(t)),K^I(x,X(t))分别表示不确定性影响下的总体质量、总体阻尼和总体刚度矩阵历程。

对上述连续时刻下的平衡方程进行时间离散化处理，结合Newmark迭代求解算法和区间差分格式，实现任意离散点处临界许用载荷的上下界计算，即：

进而，获得完整寿命期内复材层合结构临界许用载荷的区间历程，即上界和下界上式中，δ^I(x,X(t))的计算采用节点位移的差分迭代格式，针对确定性问题，δ^I(x,X(t))采用退化的δ(t_k)表示，计算如下：

$\left(\begin{array}{c}\delta \left(t_k\right)=\underset{j=1}{\overset{k}{\Sigma }}\frac{1}{2}\Delta t·\left(\stackrel{·}{\delta }\right(t_{j-1})+\stackrel{·}{\delta }(t_j\left)\right)+\delta \left(0\right)=\underset{j=1}{\overset{k}{\Sigma }}\frac{1}{2}\Delta t·\left(\underset{i=1}{\overset{j-1}{\Sigma }}\right(\stackrel{··}{\delta }\left(t_i\right)·\Delta t)+\underset{i=1}{\overset{j}{\Sigma }}(\stackrel{··}{\delta }\left(t_i\right)·\Delta t\left)\right)+\stackrel{·}{\delta }\left(0\right)·\Delta t+\delta \left(0\right)\\ =\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j-1}{\Sigma }}\stackrel{··}{\delta }\left(t_i\right)·{\mathrm{\Delta t}}^2+\frac{1}{2}\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j-1}{\Sigma }}\stackrel{··}{\delta }\left(t_i\right)·{\mathrm{\Delta t}}^2+\stackrel{·}{\delta }\left(0\right)·\Delta t+\delta \left(0\right)\end{array}\right)$

其中，t_k表示第k个离散时刻，即t_k＝kΔt，δ(t_k)表示t_k时刻下的节点位移向量，δ(0)，表示初始位移向量、速度向量和加速度向量，j和k分别表示迭代离散过程的计数指标。并且：

${\delta }^L\left(t_k\right)=\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j-1}{\Sigma }}{\stackrel{··}{\delta }}^L\left(t_i\right)·{\mathrm{\Delta t}}^2+\frac{1}{2}\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j}{\Sigma }}{\stackrel{··}{\delta }}^L\left(t_i\right)·{\mathrm{\Delta t}}^2+\stackrel{·}{\delta }\left(0\right)·\Delta t+\delta \left(0\right)$

和

${\delta }^U\left(t_k\right)=\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j-1}{\Sigma }}{\stackrel{··}{\delta }}^U\left(t_i\right)·{\mathrm{\Delta t}}^2+\frac{1}{2}\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j}{\Sigma }}{\stackrel{··}{\delta }}^U\left(t_i\right)·{\mathrm{\Delta t}}^2+\stackrel{·}{\delta }\left(0\right)·\Delta t+\delta \left(0\right)$

第五步：根据实际加载历程P^I(t)，构建与步骤四获得的临界许用载荷历程间的应力-强度区间过程干涉模型，并定义时变极限状态函数如下：

引入非概率区间过程理论，实现任意离散时刻极限状态中心值G^c、半径值G^r以及任意微小时间增量[iΔt,(i+1)Δt]内协方差Cov_G(iΔt,(i+1)Δt)和相关系数ρ_G(iΔt,(i+1)Δt)的显式表达。这里，时变极限状态函数G(t,x,X(t))的定义及相关性特征的求解过程如下：

极限状态函数G(t,x,X(t))同样被定义为一个区间过程，即G(t,x,X(t))∈G^I(t,x,X(t))。对于任意给定瞬时时刻t_i＝iΔt，极限状态将退化为一个区间量G(t_i,x,X(t_i))，有限个离散区间量的可行范围被约束在一个超立方体域Ωⁿ内。引入非概率区间过程理论，实现任意离散时刻极限状态中心值G^c、半径值G^r以及任意微小时间增量[iΔt,(i+1)Δt]内协方差Cov_G(iΔt,(i+1)Δt)和相关系数ρ_G(iΔt,(i+1)Δt)的显式表达。其中，Cov_G(iΔt,(i+1)Δt)和ρ_G(iΔt,(i+1)Δt)的定义需借助标准化处理手段，并转换工作坐标系至(U₁,U₂)，即：

G(t_i,x,X(t_i))∈[G^L(t_i,x,X(t_i)),G^U(t_i,x,X(t_i))]

＝G^c(t_i,x,X(t_i))+G^r(t_i,x,X(t_i))×U₁

G(t_i+1,x,X(t_i+1))∈[G^L(t_i+1i,x,X(t_i+1)),G^U(t_i+1,x,X(t_i+1))]

＝G^c(t_i+1,x,X(t_i+1))+G^r(t_i+1,x,X(t_i+1))×U₂

其中，G^c(t_i,x,X(t_i))和G^r(t_i,x,X(t_i))表示t_i时刻极限状态函数的均值和半径，G^c(t_i+1,x,X(t_i+1))和G^r(t_i+1,x,X(t_i+1))表示t_i+1时刻极限状态函数的均值和半径。

第六步：将首次穿越理论与步骤五中建立的复合材料层合结构极限状态的应力-强度区间过程干涉模型相结合，可以得到任意时间区间内的穿越可能度：其中，Pos{·}表示事件发生的可能性度量，E_iΔt表示穿越事件，即事件A:(iΔt)时刻结构安全G(iΔt)>0与事件B:((i+1)Δt)时刻结构失效G((i+1)Δt)<0的交事件，Δt表示微小时间增量。

第七步：遍历所有时间段内的穿越可能度Pos{E_iΔt}，定义复合材料层板结构的非概率动力可靠度计算指标：

其中，R_s(T)表示整个生命周期T内的动力可靠度，G(iΔt)为时刻iΔt的极限状态函数，和P^I(iΔt)分别对应时刻iΔt的临界许用载荷和真实加载，求解上式即实现复材层合结构动力安全态势的有效评估。这里，R_s(T)的求解计算需借助时间离散化方法进行简化近似，通过遍历每一个微小时间增量内结构发生穿越破坏的可能性指标，经叠加运算后得到可靠度的等效表达式如下：

$R_s\left(T\right)=1-P_f\left(T\right)=1-\left(Pos\right(0)+\underset{i=1}{\overset{t_i=T}{\Sigma }}Pos\{E_{i\Delta t}\left\}\right)$

其中，P_f(T)表示整个生命周期内的失效度，表示结构在初始时刻即发生失效的可能度，Δt＝t_i+1-t_i表示微小时间增量，其取值设定为整个服役寿命的1/1000。

本发明与现有技术相比的优点在于：本发明提供了含静动力不确定性复合材料层合结构可靠性分析的新思路，弥补和完善了传统基于概率理论的静态可靠性分析方法的局限性。所构建的动力可靠性度量模型，一方面可大幅减小对样本信息的依赖性，另一方面可有效计及并量化材料退化与载荷交变作用下结构安全性能的折减效应。在对结构极限载荷进行定界计算时，高效稳健的迭代算法与区间差分格式相结合，可确保载荷信息输入条件的可信性；基于首次穿越的可靠度计算方法更为真实地反映了时间相关性对结构安全的综合影响，为复合材料结构精细化强度分析理论的完善作出了积极贡献。

附图说明

图1是本发明针对复合材料层合结构的非概率动力可靠性评估流程图；

图2是本发明针对复合材料层合结构有限元建模过程中坐标变换示意图；

图3是本发明定义的极限状态相关性函数所对应的几何可行域示意图；

图4是本发明采用的首次穿越方法示意图；

图5是本发明提出的微小时间段内穿越失效可能度计算方法示意图；

图6是本发明实施例中复合材料层合结构的几何模型示意图；

图7是本发明实施例中临界许用载荷与真实载荷的区间过程干涉情况示意图。

具体实施方式

如图1所示，本发明提出了一种复合材料层合结构的非概率动力可靠性评估方法，包括以下步骤：

(1)根据复材层合结构的本构特征及载荷边界条件，基于最小势能原理构建有限元列式，推导出如下动力方程：其中，δ分别表示节点加速度、速度和位移，M,C,K分别为总体质量矩阵、总体阻尼矩阵和总体刚度矩阵，F为载荷列向量，N为节点个数。这里，复合材料层合板结构总体刚度矩阵的计算如下积分式表达：

$\left[K\right]={\int }_V{\left[B\right]}^T\left[D\right]\left[B\right]dV=\underset{k=1}{\overset{n}{\Sigma }}{\int }_{V_k}{\left[B\right]}^T\left[D\right]\left[B\right]{\mathrm{dV}}_k=\underset{k=1}{\overset{n}{\Sigma }}{\int }_{z_{k-1}}^{z_k}{\int }_A{\left[B\right]}^T\left[D\right]\left[B\right]dAdz$

其中，V代表全域体积，V_k表示第k层板的体积，n为层板总数，A为层板面积，B为应变矩阵。上式中，应变矩阵B＝{B_i}如下表达式：

ε＝Bδ^e＝[B₁>2>3>4]δ^e

和

$B_i=\left(\begin{array}{ccccc}\frac{\partial N_i}{\partial x}& 0& 0& z\frac{\partial N_i}{\partial x}& 0\\ 0& \frac{\partial N_i}{\partial y}& 0& 0& z\frac{\partial N_i}{\partial y}\\ \frac{\partial N_i}{\partial y}& \frac{\partial N_i}{\partial x}& 0& z\frac{\partial N_i}{\partial y}& z\frac{\partial N_i}{\partial x}\\ 0& 0& \frac{\partial N_i}{\partial x}& N_i& 0\\ 0& 0& \frac{\partial N_i}{\partial y}& 0& N_i\end{array}\right)$

其中，ε为应变向量，δ^e为单元位移场，N_i表示第i个节点自由度对应的形函数。弹性矩阵D可通过下式获得：

σ_xyz＝[T]^-1σ＝[T]^-1[C]ε＝[T]^-1[C][T]^-Tε_xyz＝[D]ε_xyz

其中，σ_xyz和ε_xyz表示单元在全局坐标系下的应力场和应变场，T为坐标转置矩阵(如图2所示)。

(2)引入强度比R，结合Tsai-Wu强度失效判定准则，代入到步骤一中的动力方程中，可得到临界失效载荷的显式表达式：

F_iσ_imax+F_ijσ_imaxσ_jmax＝1(i,j＝1,2,…6)

$[F_1{\sigma }_1+F_2{\sigma }_2]R+[F_{11}{\sigma }_1^2+F_{22}{\sigma }_2^2+F_{66}{\tau }_{12}^2+2F_{12}{\sigma }_1{\sigma }_2]R^2=1$

其中，F_i和F_ij表示强度特征量，σ＝[σ₁,σ₂,τ₁₂]^T为应力向量，σ_imax对应最危险点处的应力值，角标i和j分别表示单元对应的自由度计数指标。这里所述的强度比R其物理意义为极限载荷与真实载荷的比值，具体公式为：

$R=\frac{-Coe\_B+\sqrt{Coe\_B^2+4Coe\_A}}{2Coe\_A}$

其中，分别定义中间参量Coe_A＝F_ijσ_iσ_j和Coe_B＝F_iσ_i，σ_i和σ_j分别对应计数指标i和j处的单元应力。强度特征量F_i和F_ij取决于材料工程强度参数{X_t,X_c,Y_t,Y_c,S₂}，分别表示材料纵向拉伸强度、纵向压缩强度、横向拉伸强度、横向压缩强度以及剪切强度，即：

$F_1=\frac{1}{X_t}-\frac{1}{X_c},F_2=\frac{1}{Y_t}-\frac{1}{Y_c},F_{11}=\frac{1}{X_t{X}_c}$

$F_{22}=\frac{1}{Y_t{Y}_c},F_{66}=\frac{1}{S_2},F_{12}=-\frac{1}{2\sqrt{X_t{X}_c{Y}_t{Y}_c}}$

(3)综合考虑存在于第一步和第二步中结构参数中的不确定性效应，不确定性效应包括材料分散性引起的M,C,K的变动、载荷F的未确知性引起响应δ的变化、强度判定准则的模糊性导致强度特征量F_i和F_ij的不一致。引入区间向量x∈x^I＝(x₁,x₂,...,x_m)以及区间过程模型X(t)∈X(t)^I＝(X₁(t),X₂(t),...,X_n(t))，得到有限样本信息条件下结构静动力不确定性参数的数学表达：

$x^{U/L}=({x_1}^{U/L},{x_2}^{U/L},...,x_m^{U/L})=({x_1}^c\pm {x_1}^r,{x_2}^c\pm {x_2}^r,...,x_m^c\pm x_m^r)$

X(t)^U/L＝(X₁(t)^U/L,X₂(t)^U/L,...,X_n(t)^U/L)＝(X₁(t)^c±X₁(t)^r,X₂(t)^c±X₂(t)^r,...,X_m(t)^c±X_m(t)^r)其中，上标U代表参量的取值上界，上标L代表参量的取值下界，上标c代表中心值，上标r代表半径。静动力不确定性参数包括：载荷边界参数、材料特征参数、测量精度参数以及设计准则(设计许用值)参数等，为便于计算，通常表示为标准化形式，具体如下：

x＝[x^L,x^U]＝[x^c-x^r,x^c+x^r]

＝x^c+x^r[-1,1]

＝x^c+x^r×ξ₁

和

X(t)＝[X^L(t),X^U(t)]＝[X^c(t)-X^r(t),X^c(t)+X^r(t)]

＝X^c(t)+X^r(t)[-1,1]

＝X^c(t)+X^r(t)×ξ₂

其中，x^c＝(x₁^c,x₂^c,...,x_m^c)和x^r＝(x₁^r,x₂^r,...,x_m^r)表示区间向量x的均值和半径，X^c(t)＝(X₁^c(t),X₂^c(t),...,x_n^c(t))和X^c(t)＝(X₁^r(t),X₂^r(t),...,x_n^r(t))表示区间过程X(t)的均值和半径；ξ₁,ξ₂∈Ξ⁴，Ξ⁴定义为所有元素包含在[-1,1]内的4维向量集合，符号“×”定义为两个向量各对应元素相乘的算子，乘积仍为维数为4的向量。

(4)将第三步中表征的静动力不确定参数代入到第一步和第二步的方程中，形成考虑时间效应的区间集格式，即给定任意时刻有：

$\left(\begin{array}{c}{\left[M^I(x,X(t\left)\right)\right]}_{5N\times 5N}{\left\{{\stackrel{··}{\delta }}^I(x,X(t\left)\right)\right\}}_{5N}+{\left[C^I(x,X(t\left)\right)\right]}_{5N\times 5N}{\left\{{\stackrel{·}{\delta }}^I(x,X(t\left)\right)\right\}}_{5N}\\ +{\left[K^I(x,X(t\left)\right)\right]}_{5N\times 5N}{\left\{{\delta }^I(x,X(t\left)\right)\right\}}_{5N}={\left\{F^I(x,X(t\left)\right)\right\}}_{5N}\end{array}\right)$

其中，δ^I(x,X(t))分别为考虑不确定性影响下的加速度、速度和位移历程，M^I(x,X(t)),C^I(x,X(t)),K^I(x,X(t))分别表示不确定性影响下的总体质量、总体阻尼和总体刚度矩阵历程。

对上述连续时刻下的平衡方程进行时间离散化处理，结合Newmark迭代求解算法和区间差分格式，实现任意离散点处临界许用载荷的上下界计算，即：

进而，获得完整寿命期内复材层合结构临界许用载荷的区间历程，即上界和下界上式中，δ^I(x,X(t))的计算采用节点位移的差分迭代格式，针对确定性问题，δ^I(x,X(t))采用退化的δ(t_k)表示，计算如下：

$\left(\begin{array}{c}\delta \left(t_k\right)=\underset{j=1}{\overset{k}{\Sigma }}\frac{1}{2}\Delta t·\left(\stackrel{·}{\delta }\right(t_{j-1})+\stackrel{·}{\delta }(t_j\left)\right)+\delta \left(0\right)=\underset{j=1}{\overset{k}{\Sigma }}\frac{1}{2}\Delta t·\left(\underset{i=1}{\overset{j-1}{\Sigma }}\right(\stackrel{··}{\delta }\left(t_i\right)·\Delta t)+\underset{i=1}{\overset{j}{\Sigma }}(\stackrel{··}{\delta }\left(t_i\right)·\Delta t\left)\right)+\stackrel{·}{\delta }\left(0\right)·\Delta t+\delta \left(0\right)\\ =\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j-1}{\Sigma }}\stackrel{··}{\delta }\left(t_i\right)·{\mathrm{\Delta t}}^2+\frac{1}{2}\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j-1}{\Sigma }}\stackrel{··}{\delta }\left(t_i\right)·{\mathrm{\Delta t}}^2+\stackrel{·}{\delta }\left(0\right)·\Delta t+\delta \left(0\right)\end{array}\right)$

其中，t_k表示第k个离散时刻，即t_k＝kΔt，δ(t_k)表示t_k时刻下的节点位移向量，δ(0)，表示初始位移向量、速度向量和加速度向量，j和k分别表示迭代离散过程的计数指标。并且：

${\delta }^L\left(t_k\right)=\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j-1}{\Sigma }}{\stackrel{··}{\delta }}^L\left(t_i\right)·{\mathrm{\Delta t}}^2+\frac{1}{2}\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j}{\Sigma }}{\stackrel{··}{\delta }}^L\left(t_i\right)·{\mathrm{\Delta t}}^2+\stackrel{·}{\delta }\left(0\right)·\Delta t+\delta \left(0\right)$

和

${\delta }^U\left(t_k\right)=\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j-1}{\Sigma }}{\stackrel{··}{\delta }}^U\left(t_i\right)·{\mathrm{\Delta t}}^2+\frac{1}{2}\underset{j=1}{\overset{k}{\Sigma }}\underset{i=1}{\overset{j}{\Sigma }}{\stackrel{··}{\delta }}^U\left(t_i\right)·{\mathrm{\Delta t}}^2+\stackrel{·}{\delta }\left(0\right)·\Delta t+\delta \left(0\right)$

根据节点位移的动态边界信息，可知应力场界限为：

σ^L(t_k)＝min{D(x,X(t))B(x,X(t))δ^L(t_k)}和σ^U(t_k)＝max{D(x,X(t))B(x,X(t))δ^U(t_k)}

因此，任意离散点处临界许用载荷的上下界可进一步表示为：

和

(5)根据实际加载历程P^I(t)，构建与步骤四获得的临界许用载荷历程间的应力-强度区间过程干涉模型，并定义时变极限状态函数如下：

引入非概率区间过程理论，实现任意离散时刻极限状态中心值G^c、半径值G^r以及任意微小时间增量[iΔt,(i+1)Δt]内协方差Cov_G(iΔt,(i+1)Δt)和相关系数ρ_G(iΔt,(i+1)Δt)的显式表达。这里，时变极限状态函数G(t,x,X(t))的定义及相关性特征的求解过程如下：

极限状态函数G(t,x,X(t))同样被定义为一个区间过程，即G(t,x,X(t))∈G^I(t,x,X(t))。对于任意给定瞬时时刻t_i＝iΔt，极限状态将退化为一个区间量G(t_i,x,X(t_i))，有限个离散区间量的可行范围被约束在一个超立方体域Ωⁿ内。

基于区间数学理论，分别定义出任意给定时刻t_i下中心值G^c(t_i,x,X(t_i))和半径G^r(t_i,x,X(t_i))如下：

$G^c(t_i,x,X(t_i\left)\right)=\frac{G^U(t_i,x,X(t_i\left)\right)+G^L(t_i,x,X(t_i\left)\right)}{2}$

$G^r(t_i,x,X(t_i\left)\right)=\frac{G^U(t_i,x,X(t_i\left)\right)-G^L(t_i,x,X(t_i\left)\right)}{2}$

为方便起见，定义方差函数为：

$D_G\left(t_i\right)={\left(G^r\right(t_i,x,X\left(t_i\right)\left)\right)}^2={\left(\frac{G^U(t_i,x,X(t_i\left)\right)-G^L(t_i,x,X(t_i\left)\right)}{2}\right)}^2$

利用正则化手段，首先得到：

G(t_i,x,X(t_i))∈[G^L(t_i,x,X(t_i)),G^U(t_i,x,X(t_i))]

＝G^c(t_i,x,X(t_i))+G^r(t_i,x,X(t_i))×U₁

G(t_i+1,x,X(t_i+1))∈[G^L(t_i+1i,x,X(t_i+1)),G^U(t_i+1,x,X(t_i+1))]

＝G^c(t_i+1,x,X(t_i+1))+G^r(t_i+1,x,X(t_i+1))×U₂

其中，G^c(t_i,x,X(t_i))和G^r(t_i,x,X(t_i))表示t_i时刻极限状态函数的均值和半径，G^c(t_i+1,x,X(t_i+1))和G^r(t_i+1,x,X(t_i+1))表示t_i+1时刻极限状态函数的均值和半径。

从几何角度不难发现，存在无数多个不同形状的偏转矩形域包含于标准方形域内，而这些矩形域形状的改变与其对应的极限状态相关性具有映射关系(如图3所示)。于是，定义协方差函数Cov_G(iΔt,(i+1)Δt)和相关系数函数ρ_G(iΔt,(i+1)Δt)如下：

$\begin{array}{c}{\mathrm{Cov}}_G(i\Delta t,(i+1\left)\Delta t\right)=Cov(U_1,U_2)\times G^r(t_i,x,X(t_i\left)\right)\times G^r(t_{i+1},x,X(t_{i+1}\left)\right)\\ =(1-\sqrt{2}d)\times G^r(t_i,x,X(t_i\left)\right)\times G^r(t_{i+1},x,X(t_{i+1}\left)\right)\end{array},0\le d\le \sqrt{2}$

$\begin{array}{c}{\rho }_G(i\Delta t,(i+1\left)\Delta t\right)=\frac{{\mathrm{Cov}}_G(i\Delta t,(i+1\left)\Delta t\right)}{\sqrt{D_G\left(t_i\right)}\sqrt{D_G\left(t_{i+1}\right)}}\\ =\frac{Cov(U_1,U_2)}{\sqrt{D_{U_1}}\sqrt{D_{U_2}}}={\rho }_{U_1{U}_2}=1-\sqrt{2}d\end{array},0\le d\le \sqrt{2}$

其中，d表示如图3所示矩形域边长的一半，和分别是标准区间变量U₁和U₂的方差ρ_G(iΔt,(i+1)Δt)是一个无量纲量，其大小代表了G(t_i)和G(t_i+1)的线性相关度。

综上，本发明实现了对复合材料层合结构强度失效的定量表征，为后续动力可靠性建模及求解提供了必要的理论依据。

(6)将首次穿越理论(如图4所示)与与步骤五中建立的复合材料层合结构极限状态的应力-强度区间过程干涉模型相结合，可以得到任意时间区间内的穿越可能度：其中，Pos{·}表示事件发生的可能性度量，E_iΔt表示穿越事件，即事件A:(iΔt)时刻结构安全G(iΔt)>0与事件B:((i+1)Δt)时刻结构失效G((i+1)Δt)<0的交事件，Δt表示微小时间增量_。这里_，引入面积比思想(如图5所示)，Pos{E_iΔt}可定义为穿越几何条件与极限状态可行域的干涉面积与过程中总可行域(偏转矩形)面积之比，即：

上式中，几何边界条件为：和的计算通常是一个分段函数，需要结合几何边界与可行域的相交条件分类讨论。

(7)遍历所有时间段内的穿越可能度Pos{E_iΔt}，定义复合材料层板结构的非概率时变可靠度计算指标：

其中，R_s(T)表示整个生命周期T内的动力可靠度，G(iΔt)为时刻iΔt的极限状态函数，和P^I(iΔt)分别对应时刻iΔt的临界许用载荷和真实加载，求解上式即实现复材层合结构动力安全态势的有效评估。这里，R_s(T)的求解计算需借助时间离散化方法进行简化近似，通过遍历每一个微小时间增量内结构发生穿越破坏的可能性指标，经叠加运算后得到可靠度的等效表达式如下：

$R_s\left(T\right)=1-P_f\left(T\right)=1-\left(Pos\right(0)+\underset{i=1}{\overset{t_i=T}{\Sigma }}Pos\{E_{i\Delta t}\left\}\right)$

其中，P_f(T)表示整个生命周期内的失效度，表示结构在初始时刻即发生失效的可能度，Δt＝t_i+1-t_i表示微小时间增量，其取值设定为整个服役寿命的1/1000。

实施例：

为了更充分地了解该发明的特点及其对工程实际的适用性，本发明针对如图6所示24层复合材料层合板结构进行非概率动力可靠性分析。其采用对称铺层方案[θ/θ/θ/θ/θ/θ/-θ/-θ/-θ/-θ/-θ/-θ]_对称。该结构受到作用于几何中心的集中载荷，并采用四边固支加以约束。板的长度l和宽度b均为100mm，每层板厚度为t＝0.147mm。层板材料具有横向各向同性，其密度为ρ＝1.38×10³kg/m³。该层合板工程强度参数信息如表1所示，表2列出了层板模量的非概率静动力不确定性特征。

表1

表2

根据本发明提出的方法可以快速获得临界许用载荷的区间过程表达式，即：分别讨论θ＝15°和θ＝45°两种铺层形式，对应的载荷工况分别为：

于是，可获得上述两种动态载荷工况复材层合结构强度约束下的极限状态函数G(t)|_θ＝15°和G(t)|_θ＝45°。图7显示了临界许用载荷与真实载荷的区间过程干涉情况。基于前述的动力可靠性评价方法，两种工况对应的可靠度结果分别为：(1)θ＝15°,R_s(T)＝0.53；(2)θ＝45°,R_s(T)＝0.85。其中，全寿命周期为T＝20年。从结果可以看出，不同铺层形式下复合材料层合结构的强度性能差异明显，也为后续最优设计提供了较大的控制裕度。此外，动力可靠度可以有效量化结构的安全性，可靠度越高，结构越安全，反之越危险。

综上所述，本发明提出了一种复合材料层合结构的非概率动力可靠性评估方法。首先，基于有限元思想和最小势能原理，构建并解析结构的平衡方程；其次，引入非概率静动力不确定性参量，结合数值迭代算法与有限差分思想，实现复材层合结构临界许用载荷的区间过程表述，进而构建极限状态函数的显式表达式；基于首次穿越理论和离散化策略，完成任意微小时间段内穿越失效可能度的几何定义与求解；最后，通过遍历所有时间区间的穿越可能度函数，叠加运算后计算动力可靠度指标，构造出复合材料层合结构动力安全的量化判据。

以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制；其可扩展应用于含缺陷结构的优化设计领域，凡采用等同变换或者等效替换而形成的技术方案，均落在本发明权利保护范围之内。

本发明未详细阐述部分属于本领域技术人员的公知技术。