CN104809301A - 一种反映非线性刚柔混合连接特性的时域子结构方法 - Google Patents

Abstract

本发明公开的一种反映非线性刚柔混合连接特性的时域子结构方法，涉及动力学模拟方法及算法，属于结构动力学技术领域。本发明基本思路为：获得各个子结构的脉冲响应函数矩阵；根据子结构间的连接关系建立子结构界面的相容条件以及连接件的运动方程；利用界面力相容条件和脉冲响应函数矩阵建立子结构的运动方程；利用位移相容条件和连接件的运动方程将所有子结构的运动方程综合起来，解得各个子结构的响应和界面力，完成结构动力学模拟。本发明实现使传统IBS方法真实反应子结构连接件的非线性刚柔混合连接特性，避免整个系统的动力学响应分析失真，进而提高航天技术领域结构动力学模拟精度。此外，本发明拓展了传统IBS方法的应用范围。

Description

一种反映非线性刚柔混合连接特性的时域子结构方法

技术领域

本发明涉及一种动力学模拟方法，特别涉及一种反映非线性刚柔混合连接特性的动力学模拟方法，属于结构动力学技术领域。

背景技术

随着航天技术的发展和对工程结构动态设计要求的提高，航天器结构系统变得日益复杂和庞大。在动态分析和优化设计过程中，由于模型自由度较多，不得不耗费大量的时间进行计算。另一方面，许多航天器在研制过程中往往需要不同地域甚至不同国家之间的合作，考虑到技术保护问题，合作双方无法直接共享有限元模型。动态子结构方法就是为解决上述问题而发展起来的一种较为理想的方法。

自从1960年Hurty首次实现动态子结构(Dynamic Substituting，DS)技术以来，经过半个世纪的发展动态子结构技术已经广泛应用于工程领域，先后形成了三类方法：模态综合(Component Mode Synthesis，CMS)法、频域子结构(Frequency Based Substructuring，FBS)方法和基于脉冲的子结构(Impulse BasedSubstituting，IBS)即经典时域子结构方法。前两种方法中，子结构的动力学特性分别由模态和频率响应函数描述。IBS方法是Rixen于2010年提出的新型时域子结构方法，子结构的动力学特性由脉冲响应函数描述。IBS相对于其它两种方法更适于瞬态冲击动力学问题，目前已在海上风力发电机的振动分析和月球探测器的软着陆动响应计算中得到了应用。

需要注意的是IBS方法中假设子结构之间的连接是刚性的，即两个子结构的连接自由度的位移始终是相等的。连接件在系统的动响应中起重要作用，一个高精度的结构动力学分析过程不仅需要建立准确的子结构模型，也要如实地反映子结构间的连接。在实践中，所有类型的连接都是半刚性或有弹性的，如果刚性连接不能充分描述真实的连接，那么仿真计算得到的动响应结果必然与试验结果相差甚远。董威利于2012年提出了一种基于脉冲响应函数的弹性连接子结构综合方法，但是对弹性连接件的非线性特性没有进行考虑。因此连接件的阻尼和非线性特性必须被进一步了解以更好的设计结构以满足任务要求。

发明内容

针对航天技术领域结构动力学模拟使用的传统IBS方法不能真实反应子结构连接件的非线性刚柔混合连接特性，从而导致整个系统的动力学响应分析失真的问题，本发明公开的一种反映非线性刚柔混合连接特性的时域子结构方法，要解决的技术问题是在航天技术领域结构动力学模拟过程中，实现使传统IBS方法真实反应子结构连接件的非线性刚柔混合连接特性，避免整个系统的动力学响应分析失真，进而提高航天技术领域结构动力学模拟精度。

本发明的目的是通过下述技术方案实现的：

本发明基本思路为：首先获得各个子结构的脉冲响应函数矩阵；然后根据子结构间的连接关系建立子结构界面的相容条件(包括界面位移相容条件和界面力相容条件)以及连接件的运动方程；接下来利用界面力相容条件和脉冲响应函数矩阵建立子结构的运动方程；最后利用位移相容条件和连接件的运动方程将所有子结构的运动方程综合起来，解得各个子结构的响应和界面力。

本发明公开的一种反映非线性刚柔混合连接特性的时域子结构方法，包括以下步骤：

步骤1：定义基本未知量，分别是系统的第s个子结构的位移函数列向量u^(s)(t)、速度函数列向量加速度函数列向量ü^(s)(t)以及连接件对系统的作用力函数列向量λ(t)。通过数值积分方法Newmark法或者试验方法获得每个子结构的脉冲响应函数矩阵H^(s)(t)，其中s表示子结构的序号，H^(s)(t)矩阵的元素表示系统第s个子结构的第j个自由度在脉冲激励下第i个自由度的位移响应。

步骤2：建立连接件的动力学方程，包括以下步骤：

步骤2.1：将系统位移向量u分块为描述各自子结构的位移向量u^(s)，具体分块方法为：

系统的位移向量用全体子结构的位移向量表示为公式(1)：

$u={\left[\begin{array}{cccc}{u}^{{\left(1\right)}^T}& u^{{\left(2\right)}^T}& ...& u^{{\left(N_s\right)}^T}\end{array}\right]}^T---\left(1\right)$

式中，N_s为系统中子结构的个数，T表示矩阵转置。根据子结构间的连接关系，将子结构自由度u^(s)分为界面自由度和内部自由度为公式(2)：

$u^{\left(s\right)}={\left[\begin{array}{cc}{u}_c^{{\left(s\right)}^T}& u_i^{{\left(s\right)}^T}\end{array}\right]}^T---\left(2\right)$

将界面自由度分为刚性连接自由度弹性连接自由度为公式(3)：

$u^{\left(s\right)}=\left[\begin{array}{ccc}{u}_{\mathrm{cr}}^{\left(s\right)}& u_{\mathrm{ce}}^{\left(s\right)}& u_i^{\left(s\right)}\end{array}\right]---\left(3\right)$

步骤2.2：确定界面位移相容条件：

定义符号型布尔矩阵B，将系统自由度映射到界面自由度，具体方法为：子结构间所有弹性连接件的自由度u_e全部分布于各个连接件的边界上，位移相容条件要求u_e＝u_c；通过公式(4)求得布尔矩阵B：

$\mathrm{Bu}=\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}^{\left(s\right)}{u}^{\left(s\right)}=u_c---\left(4\right)$

则界面位移相容条件最终表示为公式(5)：

$\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}^{\left(s\right)}{u}^{\left(s\right)}=u_e---\left(5\right)$

步骤2.3：确定界面力相容条件：

连接件自由度u_e对全体子结构边界的作用力向量用λ表示，则系统第s个子结构受到的来自连接件的作用力g^(S)为公式(6)：

$g^{\left(s\right)}=B^{{\left(s\right)}^T}\mathrm{\λ}---\left(6\right)$

界面力相容条件(即牛顿第三定律)要求子结构边界对连接件自由度u_e的作用力向量为λ_e＝-λ；

步骤2.4：建立连接件的运动方程描述，具体方法为：

令M_j、g_j是总体质量矩阵和所有弹性部件的作用力矢量，则所有连接件的运动方程为公式(7)：

$M_j{\stackrel{\·\·}{u}}_e+g_j({\stackrel{\·}{u}}_e,u_e)=-{\mathrm{\λ}}_e---\left(7\right)$

式中：ü_e、和u_e分别是连接件自由度的加速度、速度和位移向量。

步骤3：建立子结构的运动方程，具体实现方法，包括以下步骤：

步骤3.1：建立时间连续形式的子结构运动方程：

由Duhamel积分可知，系统第s个子结构Ω^(s)在外载荷f^(s)、刚性连接界面力和弹性连接界面力作用下的位移u^(s)可表示为公式(8)：

$u^{\left(s\right)}\left(t\right)={\∫}_0^t{H}^{\left(s\right)}(t-\mathrm{\τ})(f^{\left(s\right)}\left(\mathrm{\τ}\right)+g_{\mathrm{cr}}^{\left(s\right)}\left(\mathrm{\τ}\right)+g_{\mathrm{ce}}^{\left(s\right)}\left(\mathrm{\τ}\right))\mathrm{d\τ}---\left(8\right)$

步骤3.2：将步骤3.1中的子结构运动方程式(8)进行时间离散，得到公式(9)：

$\begin{array}{c}{u}_n^{\left(s\right)}=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{h}_{n-i}^{\left(s\right)}[f_i^{\left(s\right)}+B_{\mathrm{cr}}^{{\left(s\right)}^T}({\mathrm{\λ}}_{\mathrm{vr},i}+{\mathrm{\λ}}_{\mathrm{cr},i+1})]\frac{\mathrm{dt}}{2}\\ +\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{H}_{n-i}^{\left(s\right)}[f_{i+1}^{\left(s\right)}+B_{\mathrm{ce}}^{{\left(s\right)}^T}({\mathrm{\λ}}_{\mathrm{ce},i}+{\mathrm{\λ}}_{\mathrm{ce},i+1})]\frac{\mathrm{dt}}{2}\end{array}---\left(9\right)$

式中，dt为积分步长，角标代表时刻(如u_n＝u(ndt))；

步骤3.3：利用Newmark方法描述子结构的速度和加速度ü^(s)，如公式(10)：

$\left\{\begin{array}{c}{\stackrel{\·}{u}}_n=\frac{\mathrm{\γ}}{\mathrm{\βdt}}(u_n-u_{n-1})+(1-\frac{\mathrm{\γ}}{\mathrm{\β}}){\stackrel{\·}{u}}_{n-1}+(1-\frac{\mathrm{\γ}}{2\mathrm{\β}})\mathrm{dt}{\stackrel{\·\·}{u}}_{n-1}\\ {\stackrel{\·\·}{u}}_n=\frac{1}{\mathrm{\β}{\mathrm{dt}}^2}(u_n-u_{n-1})-\frac{1}{\mathrm{\βdt}}{\stackrel{\·}{u}}_{n-1}+(1-\frac{1}{2\mathrm{\β}}){\stackrel{\·\·}{u}}_{n-1}\end{array}\right.---\left(10\right)$

式中，γ和β是Newmark法的无量纲参数；

步骤4：根据公式(5)、(7)、(9)、(10)求解各个子结构子结构间界面力λ_n和位移响应的时间递推公式(11)和递推迭代收敛条件，根据时间递推公式(11)和递推迭代收敛条件完成结构动力学模拟。具体方法为：

定义

$B^{\left(s\right)}=\left[\begin{array}{c}{B}_{\mathrm{cr}}^{\left(s\right)}\\ B_{\mathrm{ce}}^{\left(s\right)}\end{array}\right],\mathrm{\λ}=\left[\begin{array}{c}{\mathrm{\λ}}_{\mathrm{cr}}\\ {\mathrm{\λ}}_{\mathrm{ce}}\end{array}\right]$

其中为界面刚性连接件和非线性弹性连接件自由度的映射矩阵；λ_cr、λ_ce为界面刚性和非线性弹性连接件的界面力。

各个子结构间界面力λ_n和子结构位移响应的时间递推公式为公式(11)：

$\left\{\begin{array}{c}{\mathrm{\λ}}_n=-G^{-1}\left[\begin{array}{c}{p}_{n-1}\\ q_n\end{array}\right]\\ u_n^{\left(s\right)}={\tilde{u}}_{n-1}^{\left(s\right)}+H_1^{\left(s\right)}{B}^{{\left(s\right)}^T}{\mathrm{\λ}}_n\frac{\mathrm{dt}}{2}\end{array}\right.---\left(11\right)$

式中：

$G=\frac{\mathrm{dt}}{2}\left[\begin{array}{cc}\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_{\mathrm{cr}}^{\left(s\right)}{H}_1^{\left(s\right)}{B}_{\mathrm{cr}}^{{\left(s\right)}^T}& \underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_{\mathrm{cr}}^{\left(s\right)}{H}_1^{\left(s\right)}{B}_{\mathrm{ce}}^{{\left(s\right)}^T}\\ 0_{\mathrm{er}}& \frac{2I_{\mathrm{ee}}}{\mathrm{dt}}\end{array}\right]---\left(12\right)$

$\begin{array}{c}{\tilde{u}}_{n-1}^{\left(s\right)}=\underset{i=0}{\overset{n-2}{\mathrm{\Σ}}}{H}_{n-i}^{\left(s\right)}[f_i^{\left(s\right)}+f_{i+1}^{\left(s\right)}+B^{{\left(s\right)}^T}({\mathrm{\λ}}_i+{\mathrm{\λ}}_{i+1})]\frac{\mathrm{dt}}{2}\\ +H_1^{\left(s\right)}(f_{n-1}^{\left(s\right)}+f_n^{\left(s\right)}+B^{{\left(s\right)}^T}{\mathrm{\λ}}_{n-1})\frac{\mathrm{dt}}{2}\end{array}---\left(13\right)$

$p_{n-1}=\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_{\mathrm{cr}}^{\left(s\right)}{\tilde{u}}_{n-1}^{\left(s\right)}---\left(14\right)$

$q_u=M_e{\stackrel{\·\·}{u}}_{\mathrm{ce},n}+g_e\left({\stackrel{\·}{u}}_{\mathrm{ce},n}{u}_{\mathrm{ce},n}\right)---\left(15\right)$

其中I_ee是单位矩阵；

求出各子结构的位移后，各子结构的速度及加速度ü^(s)响应可以根据式(10)求出，即完成各个子结构间界面力λ_n、子结构位移u^(s)、速度及加速度ü^(s)的求解，完成结构动力学模拟。

本发明的一种反映非线性刚柔混合连接特性的时域子结构方法，实现使传统IBS方法真实反应子结构连接件的非线性刚柔混合连接特性，避免整个系统的动力学响应分析失真，进而提高航天技术领域结构动力学模拟精度。

根据式(11)求解界面力λ_n，为确保求解精度，结果需满足迭代收敛条件。步骤4所述的递推迭代收敛条件为递推迭代收敛条件中各参数的求解过程为：

式(11)是关于λ_n的非线性方程，改写为残差形式和其线性化模型，如公式(16)、(17)：

$r_n={\mathrm{G\λ}}_n+\left[\begin{array}{c}{p}_{n-1}\\ q_n\end{array}\right]---\left(16\right)$

r_n+A_nΔλ_n＝0      (17)

式中：A_n称为等效切向刚度矩阵或Jacobian矩阵，Δλ_n是λ_n的增量

$A_n=\frac{\∂r_n}{{\∂\mathrm{\λ}}_n}=G+\left[\begin{array}{c}0\\ \∂q_n/\∂{\mathrm{\λ}}_n\end{array}\right]---\left(18\right)$

${\mathrm{\Δ\λ}}_n={\mathrm{\λ}}_n^{k+1}-{\mathrm{\λ}}_n^k---\left(19\right)$

$\frac{\∂q_n}{\∂{\mathrm{\λ}}_n}=(\frac{M_e}{\mathrm{\β}{\mathrm{dt}}^2}+\frac{\∂g_e}{\∂{\stackrel{\·}{u}}_{\mathrm{ce},n}}\frac{\mathrm{\γ}}{\mathrm{\βdt}}+\frac{\∂g_e}{{\∂u}_{\mathrm{ce},n}})\left(\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_{\mathrm{ce}}^{\left(s\right)}{H}_1^{\left(s\right)}{B}^{{\left(s\right)}^T}\frac{\mathrm{dt}}{2}\right)---\left(20\right)$

式中：上标表示迭代次数，比如表示第n个时间步中第k次迭代得到的界面力。

通过式(16)和式(18)分别得到第k次迭代的和后，第k+1次迭代的界面力可由式求得

${\mathrm{\λ}}_n^{k+1}={\mathrm{\λ}}_n^k-{\left(A_n^k\right)}^{-1}{r}_n^k---\left(21\right)$

迭代计算直到满足收敛条件后停止迭代，进入下一个时间步的计算。

有益效果：

1、以脉冲响应函数为基础，本发明的一种反映非线性刚柔混合连接特性的时域子结构方法，克服了经典时域子结构法仅适于分析子结构间为纯刚性或弹性连接的缺陷，实现使传统IBS方法真实反应子结构连接件的非线性刚柔混合连接特性，避免整个系统的动力学响应分析失真，进而提高航天技术领域结构动力学模拟精度。

2、本发明的一种反映非线性刚柔混合连接特性的时域子结构方法，拓展了传统IBS方法(时域子结构方法)的应用范围。

附图说明

图1为本发明的一种反映非线性刚柔混合连接特性的时域子结构方法的流程示意图。

图2为具体实施方式中的12自由度弹簧-阻尼器-质量系统。

图3为具体实施方式中质量块m₁₃的外载荷-时间曲线。

图4为具体实施方式中质量块m₁₁的加速度响应-时间曲线。

具体实施方式

实施例1：

为了更好的阐述本发明的目的和作用，下面通过对一个12自由度的弹簧-阻尼器-质量系统的动力学分析对本发明进行详细解释。分析步长取1ms，分析时间为1s。

考虑一个12自由度的弹簧-阻尼-质量块系统如图2所示。它由两个子结构组成，子结构1含7个自由度，子结构2含5个自由度，质量块m₅与m₁₁为刚性连接，m₆和m₇通过弹性界面与质量块m₁₂相连。系统的所有物理参数见表1。

表1弹簧-阻尼-质量块系统的物理参数

本实施例的一种反映非线性刚柔混合连接特性的时域子结构方法，流程图如图1所示，包括如下步骤：

步骤1：定义基本未知量，分别是系统的第s个子结构的位移函数列向量u^(s)(t)、速度函数列向量加速度函数列向量ü^(s)(t)以及连接件对系统的作用力函数列向量λ(t)；初始化基本未知量，即t＝0的初始时刻系统不受外力作用，基本未知量全部为零向量；

通过Newmark法计算每个子结构的脉冲响应函数矩阵H^(s)(t)。以H^(s)(t)中某一个元素的计算为例，就是在m₇上施加迪拉克函数δ(t)，然后通过Newmark法计算得出的m₂的位移响应；

步骤2：根据子结构之间的连接关系，求得布尔矩阵B＝[B⁽¹⁾ B⁽²⁾]，B的作用是从系统全体自由度中筛选出系统中两个子结构的界面自由度，该“界面自由度”指子结构与弹性连接件相连接的自由度；布尔矩阵B通过步骤21和步骤2.2求得：

步骤2.1：系统位移向量的分块。

系统的位移向量u用全体子结构的位移向量表示为公式(22)：

$u={\left[\begin{array}{cc}{u}^{{\left(1\right)}^T}& u^{{\left(2\right)}^T}\end{array}\right]}^T={\left[\begin{array}{cccc}{u}^1& u^2& ...& u^{15}\end{array}\right]}^T---\left(22\right)$

式中：u⁽¹⁾代表子结构1的所有自由度的位移向量，u¹代表质量块m₁的位移。

根据子结构间的连接关系，将子结构的位移向量u^(s)分块为刚性连接自由度弹性连接自由度和内部自由度如公式(23)：

$u^{\left(s\right)}={\left[\begin{array}{ccc}{u}_{\mathrm{cr}}^{\left(s\right)}& u_{\mathrm{ce}}^{\left(s\right)}& u_i^{\left(s\right)}\end{array}\right]}^T---\left(23\right)$

定义柔性连接自由度u_ce，如公式(24)：

$u_{\mathrm{ce}}={\left[\begin{array}{ccc}{u}_{\mathrm{ce}}^{\left(6\right)}& u_{\mathrm{ce}}^{\left(7\right)}& u_{\mathrm{ce}}^{\left(12\right)}\end{array}\right]}^T---\left(24\right)$

步骤2.2：确定界面位移相容条件。

弹性连接件的自由度u_e全部分布于各个连接件的边界上，由于弹性连接件与子结构之间是刚性连接的，因此位移相容条件要求弹性连接件的自由度等于子结构的边界自由度，即u_e＝u_c。通过公式(25)求得布尔矩阵B：

$\mathrm{Bu}=\underset{s=1}{\overset{2}{\mathrm{\Σ}}}{B}^{\left(s\right)}{u}^{\left(s\right)}=u_c---\left(25\right)$

求得布尔矩阵如公式(26)：

$B=\left[\begin{array}{c}{B}_{\mathrm{cr}}\\ B_{\mathrm{ce}}\end{array}\right]---\left(26\right)$

其中布尔矩阵B_cr和B_ce为

B_cr＝[0 0 0 0 1 0 0 -1 0 0 0 0]

$B_{\mathrm{ce}}=\left[\begin{array}{cccccccccccc}0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\end{array}\right]---\left(27\right)$

则界面位移相容条件最终表示为公式(28)：

$\underset{s=1}{\overset{2}{\mathrm{\Σ}}}{B}^{\left(s\right)}{u}^{\left(s\right)}=u_e---\left(28\right)$

步骤3：根据系统结构及参数计算得出系统弹性连接件的质量矩阵M_e、连接元件柔性部分的内力g_e。本实例中使用的具体计算方法为结构动力学中的“影响系数法”，M_e如式(29)；

$M_e=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 5\end{array}\right]\mathrm{Kg}---\left(29\right)$

弹簧k₈和k₉具有式(30)所示立方刚度

g_s＝kx³       (30)

阻尼器c₈和c₉具有式(31)所示平方阻尼

$g_d=c\stackrel{\·}{x}\left|\stackrel{\·}{x}\right|---\left(31\right)$

连接元件柔性部分的内力可表示为公式(32)：

$g_e=y(\left[\begin{array}{c}{k}_8{\left({\mathrm{au}}_{\mathrm{ce}}\right)}^3\\ k_9{\left({\mathrm{bu}}_{\mathrm{ce}}\right)}^3\end{array}\right]+\left[\begin{array}{c}{c}_8{a\stackrel{\·}{u}}_{\mathrm{ce}}|{a\stackrel{\·}{u}}_{\mathrm{ce}}\\ c_{9}b{\stackrel{\·}{u}}_{\mathrm{ce}}|b{\stackrel{\·}{u}}_{\mathrm{ce}}\end{array}\right])---\left(32\right)$

式中：

$a=\left[\begin{array}{ccc}1& 0& -1\end{array}\right],b=\left[\begin{array}{ccc}0& 1& -1\end{array}\right],y=\left[\begin{array}{cc}1& 0\\ 0& 1\\ -1& -1\end{array}\right]$

步骤4：确定各个子结构间界面力λ_n和子结构位移响应的时间递推公式。取Newmark法的无量纲参数β＝0.25，γ＝0.5，步骤1中四组基本未知量的时间递推公式如式(33)：

$\left\{\begin{array}{c}{u}_{e,n}=\left(\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_e^{\left(s\right)}{H}_1^{\left(s\right)}{B}^{{\left(s\right)}^T}\frac{\mathrm{dt}}{2}\right){\mathrm{\λ}}_n+\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_e^{\left(s\right)}{\tilde{u}}_{n-1}^{\left(s\right)}\\ {\stackrel{\·}{u}}_{e,n}=\left(\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_e^{\left(s\right)}{H}_1^{\left(s\right)}{B}^{{\left(s\right)}^T}\right){\mathrm{\λ}}_n+\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_e^{\left(s\right)}\left[\frac{2}{\mathrm{dt}}({\tilde{u}}_{n-1}^{\left(s\right)}-u_{n-1}^{\left(s\right)})\right]\\ -\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_e^{\left(s\right)}{\stackrel{\·}{u}}_{n-1}^{\left(s\right)}\\ {\stackrel{\·\·}{u}}_{e,n}=\left(\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_e^{\left(s\right)}{H}_1^{\left(s\right)}{B}^{{\left(s\right)}^T}\frac{2}{\mathrm{dt}}\right){\mathrm{\λ}}_n+\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_e^{\left(s\right)}\left[\frac{4}{{\mathrm{dt}}^2}({\tilde{u}}_{n-1}^{\left(s\right)}-u_{n-1}^{\left(s\right)})\right]\\ -\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_e^{\left(s\right)}[\frac{4}{\mathrm{dt}}{\stackrel{\·}{u}}_{n-1}^{\left(s\right)}+{\stackrel{\·\·}{u}}_{n-1}^{\left(s\right)}]\\ u_n^{\left(s\right)}={\tilde{u}}_{n-1}^{\left(s\right)}+H_1^{\left(s\right)}{B}^{{\left(s\right)}^T}{\mathrm{\λ}}_n\frac{\mathrm{dt}}{2}\end{array}\right.---\left(33\right)$

其中，dt为0.001s，下角标代表时刻，如β和γ是Newmark法的无量纲参数，G、q_n以及p_n-1分别由下式确定；

$\begin{array}{c}{\tilde{u}}_{n-1}^{\left(s\right)}=\underset{i=0}{\overset{n-2}{\mathrm{\Σ}}}{H}_{n-i}^{\left(s\right)}[f_i^{\left(s\right)}+f_{i+1}^{\left(s\right)}+B^{{\left(s\right)}^T}({\mathrm{\λ}}_i+{\mathrm{\λ}}_{i+1})]\frac{\mathrm{dt}}{2}\\ +H_1^{\left(s\right)}(f_{n-1}^{\left(s\right)}+f_n^{\left(s\right)}+B^{{\left(s\right)}^T}{\mathrm{\λ}}_{n-1})\frac{\mathrm{dt}}{2}\end{array}$

$G=\left[\begin{array}{cc}{G}_{\mathrm{rr}}& G_{\mathrm{re}}\\ 0_{\mathrm{er}}& I_{\mathrm{ee}}\end{array}\right]$

$q_n=M_j{\stackrel{\·\·}{u}}_{e,n}+g_j({\stackrel{\·}{u}}_{e,n},u_{e,n})$

$p_{n-1}=\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_r^{\left(s\right)}{\tilde{u}}_{n-1}^{\left(s\right)}---\left(34\right)$

其中，

$G_{\mathrm{re}}=\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_r^{\left(s\right)}{H}_1^{\left(s\right)}{B}_e^{{\left(s\right)}^T}\frac{\mathrm{dt}}{2}$

$G_{\mathrm{rr}}=\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_r^{\left(s\right)}{H}_1^{\left(s\right)}{B}_e^{{\left(s\right)}^T}\frac{\mathrm{dt}}{2}$

步骤5：求解系统在外力f(t)作用下的响应，本实例中f(t)仅第13个元素非零，第13个元素如图2所示。外力f(t)如公式(35)。

在质量块m₁₃上施加如下载荷(单位N，作用时间1s)，分别以N刚DM和Newmark法对系统进行动力学分析，m₁₁的加速度响应对比如图4所示。

$f\left(t\right)=\left\{\begin{array}{cc}5000t& t\∈\left[\mathrm{0,0.2}\right)\\ 1000& t\∈\left[\mathrm{0.2,0.8}\right)\\ 5000-5000t& t\∈\left[\mathrm{0.8.1}\right]\end{array}\right.---\left(35\right)$

对响应点加速度进行误差分析，以Newmark方法为参照(Newmark方法为本领域公认计算精度很高的算法，但该算法计算效率低)，本节所提方法对质量块m₁₁加速度的相对峰值误差(RAAE)为0.0015。由此可见，本发明不仅能够处理非线性刚柔混合连接的子结构问题，而且拥有极高的精度。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种反映非线性刚柔混合连接特性的时域子结构方法，包括如下步骤，

步骤1：定义基本未知量，分别是系统的第s个子结构的位移函数列向量u^(s)(t)、速度函数列向量加速度函数列向量以及连接件对系统的作用力函数列向量λ(t)；通过数值积分方法Newmark法或者试验方法获得每个子结构的脉冲响应函数矩阵H^(s)(t)，其中s表示子结构的序号，H^(s)(t)矩阵的元素表示系统第s个子结构的第j个自由度在脉冲激励下第i个自由度的位移响应；

步骤2：建立连接件的动力学方程，包括以下步骤，

步骤2.1：将系统位移向量u分块为描述各自子结构的位移向量u^(s)，具体分块方法为：

系统的位移向量用全体子结构的位移向量表示为公式(1)：

$u={\left[\begin{array}{cccc}{u}^{{\left(1\right)}^T}& u^{{\left(2\right)}^T}& ...& u^{{\left(N_s\right)}^T}\end{array}\right]}^T---\left(1\right)$

式中，N_s为系统中子结构的个数，T表示矩阵转置；根据子结构间的连接关系，将子结构自由度u(s)分为界面自由度和内部自由度为公式(2)：

$u^{\left(s\right)}={\left[\begin{array}{cc}{u}_c^{{\left(s\right)}^T}& u_i^{{\left(s\right)}^T}\end{array}\right]}^T---\left(2\right)$

将界面自由度分为刚性连接自由度弹性连接自由度为公式(3)：

$u^{\left(s\right)}=\left[\begin{array}{ccc}{u}_{\mathrm{cr}}^{\left(s\right)}& u_{\mathrm{ce}}^{\left(s\right)}& u_i^{\left(s\right)}\end{array}\right]---\left(3\right)$

步骤2.2：确定界面位移相容条件：

定义符号型布尔矩阵B，将系统自由度映射到界面自由度，具体方法为：子结构间所有弹性连接件的自由度u_e全部分布于各个连接件的边界上，位移相容条件要求u_e＝u_c；通过公式(4)求得布尔矩阵B：

$\mathrm{Bu}=\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}^{\left(s\right)}{u}^{\left(s\right)}=u_c---\left(4\right)$

则界面位移相容条件最终表示为公式(5)：

$\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}^{\left(s\right)}{u}^{\left(s\right)}=u_e---\left(5\right)$

步骤2.3：确定界面力相容条件：

连接件自由度u_e对全体子结构边界的作用力向量用λ表示，则系统第s个子结构受到的来自连接件的作用力g^(s)为公式(6)：

$g^{\left(s\right)}=B^{{\left(s\right)}^T}\mathrm{\λ}---\left(6\right)$

界面力相容条件要求子结构边界对连接件自由度u_e的作用力向量为λ_e＝-λ；

步骤2.4：建立连接件的运动方程描述，具体方法为：

令M_j、g_j是总体质量矩阵和所有弹性部件的作用力矢量，则所有连接件的运动方程为公式(7)：

$M_j{\stackrel{\·\·}{u}}_e+g_j({\stackrel{\·}{u}}_e,u_e)=-{\mathrm{\λ}}_e---\left(7\right)$

式中：和u_e分别是连接件自由度的加速度、速度和位移向量；

其特征在于：还包括如下步骤，

步骤3：建立子结构的运动方程，具体实现方法，包括以下步骤，

步骤3.1：建立时间连续形式的子结构运动方程：

由Duhamel积分可知，系统第s个子结构Ω^(s)在外载荷f^(s)、刚性连接界面力和弹性连接界面力作用下的位移u^(s)可表示为公式(8)：

$u^{\left(s\right)}\left(t\right)={\∫}_0^t{H}^{\left(s\right)}(t-\mathrm{\τ})(f^{\left(s\right)}\left(\mathrm{\τ}\right)+g_{\mathrm{cr}}^{\left(s\right)}\left(\mathrm{\τ}\right)+g_{\mathrm{ce}}^{\left(s\right)}\left(\mathrm{\τ}\right))\mathrm{d\τ}---\left(8\right)$

步骤3.2：将步骤3.1中的子结构运动方程式(8)进行时间离散，得到公式(9)：

$\begin{array}{c}{u}_n^{\left(s\right)}=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{H}_{n-i}^{\left(s\right)}[f_i^{\left(s\right)}+B_{\mathrm{cr}}^{{\left(s\right)}^T}({\mathrm{\λ}}_{\mathrm{cr},i}+{\mathrm{\λ}}_{\mathrm{cr},i+1})]\frac{\mathrm{dt}}{2}\\ +\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{H}_{n-i}^{\left(s\right)}[f_{i+1}^{\left(s\right)}+B_{\mathrm{ce}}^{{\left(s\right)}^T}({\mathrm{\λ}}_{\mathrm{ce},i}+{\mathrm{\λ}}_{\mathrm{ce},i+1})]\frac{\mathrm{dt}}{2}\end{array}---\left(9\right)$

式中，dt为积分步长，角标代表时刻(如u_n＝u(ndt))；

步骤3.3：利用Newmark方法描述子结构的速度和加速度如公式(10)：

$\left\{\begin{array}{c}{\stackrel{\·}{u}}_n=\frac{\mathrm{\γ}}{\mathrm{\βdt}}(u_n-u_{n-1})+(1-\frac{\mathrm{\γ}}{\mathrm{\β}}){\stackrel{\·}{u}}_{n-1}+(1-\frac{\mathrm{\γ}}{2\mathrm{\β}})\mathrm{dt}{\stackrel{\·\·}{u}}_{n-1}\\ {\stackrel{\·\·}{u}}_n=\frac{1}{\mathrm{\β}{\mathrm{dt}}^2}(u_n-u_{n-1})-\frac{1}{\mathrm{\βdt}}{\stackrel{\·}{u}}_{n-1}+(1-\frac{1}{2\mathrm{\β}}){\stackrel{\·\·}{u}}_{n-1}\end{array}\right.---\left(10\right)$

式中，γ和β是Newmark法的无量纲参数；

步骤4：根据公式(5)、(7)、(9)、(10)求解各个子结构子结构间界面力λ_n和位移响应的时间递推公式(11)和递推迭代收敛条件，根据时间递推公式(11)和递推迭代收敛条件完成结构动力学模拟；具体方法为：

定义

$B^{\left(s\right)}=\left[\begin{array}{c}{B}_{\mathrm{cr}}^{\left(s\right)}\\ B_{\mathrm{ce}}^{\left(s\right)}\end{array}\right],\mathrm{\λ}=\left[\begin{array}{c}{\mathrm{\λ}}_{\mathrm{cr}}\\ {\mathrm{\λ}}_{\mathrm{ce}}\end{array}\right]$

其中为界面刚性连接件和非线性弹性连接件自由度的映射矩阵；λ_cr、λ_ce为界面刚性和非线性弹性连接件的界面力；

各个子结构间界面力λ_n和子结构位移响应的时间递推公式为公式(11)：

$\left\{\begin{array}{c}{\mathrm{\λ}}_n=-G^{-1}\left[\begin{array}{c}{p}_{n-1}\\ q_n\end{array}\right]\\ u_n^{\left(s\right)}={\tilde{u}}_{n-1}^{\left(s\right)}+H_1^{\left(s\right)}{B}^{{\left(s\right)}^T}{\mathrm{\λ}}_n\frac{\mathrm{dt}}{2}\end{array}\right.---\left(11\right)$

式中：

$G=\frac{\mathrm{dt}}{2}\left[\begin{array}{cc}\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_{\mathrm{cr}}^{\left(s\right)}{H}_1^{\left(s\right)}{B}_{\mathrm{cr}}^{{\left(s\right)}^T}& \underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_{\mathrm{cr}}^{\left(s\right)}{H}_1^{\left(s\right)}{B}_{\mathrm{ce}}^{{\left(s\right)}^T}\\ 0_{\mathrm{er}}& \frac{2I_{\mathrm{ee}}}{\mathrm{dt}}\end{array}\right]---\left(12\right)$

$\begin{array}{c}{\tilde{u}}_{n-1}^{\left(s\right)}=\underset{i=0}{\overset{n-2}{\mathrm{\Σ}}}{H}_{n-i}^{\left(s\right)}[f_i^{\left(s\right)}+f_{i+1}^{\left(s\right)}+B^{{\left(s\right)}^T}({\mathrm{\λ}}_i+{\mathrm{\λ}}_{i+1})]\frac{\mathrm{dt}}{2}\\ +H_1^{\left(s\right)}(f_{n-1}^{\left(s\right)}+f_n^{\left(s\right)}+B^{{\left(s\right)}^T}{\mathrm{\λ}}_{n-1})\frac{\mathrm{dt}}{2}\end{array}---\left(13\right)$

$p_{n-1}=\underset{s=1}{\overset{N_s}{\mathrm{\Σ}}}{B}_{\mathrm{cr}}^{\left(s\right)}{\tilde{u}}_{n-1}^{\left(s\right)}---\left(14\right)$

$q_n=M_e{\stackrel{\·\·}{u}}_{\mathrm{ce},n}+g_e({\stackrel{\·}{u}}_{\mathrm{ce},n},u_{\mathrm{ce},n})---\left(15\right)$

其中I_ee是单位矩阵；

求出各子结构的位移后，各子结构的速度及加速度响应可以根据式(10)求出，即完成各个子结构间界面力λ_n、子结构位移u^(s)、速度及加速度的求解，完成结构动力学模拟。

2.根据权利要求1所述的一种反映非线性刚柔混合连接特性的时域子结构方法，其特征在于：步骤4所述的递推迭代收敛条件为递推迭代收敛条件中各参数的求解过程为：

式(11)是关于λ_n的非线性方程，改写为残差形式和其线性化模型，如公式(16)、(17)：

$r_n=G{\mathrm{\λ}}_n+\left[\begin{array}{c}{p}_{n-1}\\ q_n\end{array}\right]---\left(16\right)$

r_n+A_nΔλ_n＝0   (17)

式中：A_n称为等效切向刚度矩阵或Jacobian矩阵，Δλ_n是λ_n的增量

$A_n=\frac{\∂r_n}{\∂{\mathrm{\λ}}_n}=G+\left[\begin{array}{c}0\\ \∂q_n/\∂{\mathrm{\λ}}_n\end{array}\right]---\left(18\right)$

$\mathrm{\Δ}{\mathrm{\λ}}_n={\mathrm{\λ}}_n^{k+1}-{\mathrm{\λ}}_n^k---\left(19\right)$

$\frac{\∂q_n}{\∂{\mathrm{\λ}}_n}=(\frac{M_e}{\mathrm{\β}{\mathrm{dt}}^2}+\frac{\∂g_e}{\∂{\stackrel{\·}{u}}_{\mathrm{ce},n}}\frac{\mathrm{\γ}}{\mathrm{\βdt}}+\frac{\∂g_e}{\∂u_{\mathrm{ce},n}})\left(\underset{s=1}{\overset{\stackrel{}{}{N}_s}{\mathrm{\Σ}}}{B}_{\mathrm{ce}}^{\left(s\right)}{H}_1^{\left(s\right)}{B}^{{\left(s\right)}^T}\frac{\mathrm{dt}}{2}\right)---\left(20\right)$

式中：上标表示迭代次数，比如表示第n个时间步中第k次迭代得到的界面力；

通过式(16)和式(18)分别得到第k次迭代的和后，第k+1次迭代的界面力可由式求得

${\mathrm{\λ}}_n^{k+1}={\mathrm{\λ}}_n^k-{\left(A_n^k\right)}^{-1}{r}_n^k---\left(21\right)$

迭代计算直到满足收敛条件后停止迭代，进入下一个时间步的计算。