CN105426343A - 一种基于傅里叶级数的复杂结构振动解析分析方法 - Google Patents

Abstract

本发明公开了一种基于傅里叶级数的复杂结构振动解析分析方法，步骤包括：将被分析复杂结构划分为规则的子结构；为各子结构设定控制微分方程的位移傅里叶级数；将子结构控制微分方程的位移傅里叶级数代入各子结构的控制微分方程，采用傅里叶级数展开方法获得各子结构包含质量、刚度矩阵的振动方程；根据振动方程求取各子结构的联接点至激励点或响应点的位移阻抗函数结果；依据各子结构激励点与耦合点之间的阻抗关系、耦合点处的力平衡和位移相容条件，推导得到整个复杂结构系统的位移方程，将各子结构位移阻抗函数结果代入位移方程，求得被分析复杂结构的位移响应解。

Description

一种基于傅里叶级数的复杂结构振动解析分析方法

技术领域

本发明涉及复杂结构的振动求解方法，具体涉及一种基于傅里叶级数的复杂结构振动解析分析方法，用于将傅里叶级数方法扩展到复杂结构振动的求解方法。

背景技术

傅里叶级数解法作为一种解析分析方法，是近年来受到重视的方法之一，它可以适用于各种边界条件、运算方便，并由于级数间具有正交性，可使计算量大为减小并保证很高精度。

美国韦恩州立大学李文龙提出一种广义傅里叶级数方法，该方法成功解决了任意边界条件下梁结构的弹性振动问题(W.L.Li,FREEVIBRATIONSOFBEAMSWITHGENERALBOUNDARYCONDITIONS,JournalofSoundandVibration,237(2000)709-725.)，求解出的固有频率和振型的精度及级数收敛速度都达到了十分理想的效果，该方法应用于多跨度梁结构振动问题也获得了理想的结果(W.L.Li,H.G.Xu,AnExactFourierSeriesMethodfortheVibrationAnalysisofMulti-spanBeamSystems,JournalofComputationalandNonlinearDynamics,4(2009)1-9.)。并且，该方法也被扩展到任意边界板结构振动问题的求解(W.L.Li,VIBRATIONANALYSISOFRECTANGULARPLATESWITHGENERALELASTICBOUNDARYSUPPORTS,JournalofSoundandVibration273(2004)619-635)。

文献“VIBRATIONSOFRECTANGULARPLATESWITHARBITRARYNONUNIFORMELASTICEDGERESTRAINTS(X.Zhang,WenL.Li,JournalofSoundandVibration326(2009)221-234)”中，国内的杜敬涛等人将傅里叶级数方法解决了任意边界弹性边界矩形板振动问题，不均匀边界问题，并且扩展到了板与板耦合振动分析问题。对于梁、板、圆柱壳等规则结构，这些结构的微分方程是四阶的，展开级数具有四阶(或更高阶)的逐项可导的性质，这些结构均可用广义傅里叶级数方法进行结构振动求解。但是对于不规则的复杂结构而言，广义傅里叶级数方法还无法求解，原因是一般结构的高阶逐项可导的条件很难满足，所以，傅里叶级数解法解决问题的范围受到很大限制。

发明内容

本发明要解决的技术问题：针对现有技术的上述问题，提供一种能够发挥傅里叶级数展开方法作为解析方法的收敛快速、计算所需资源少的优点，通用性好，能够推广到各种可规则化子结构的复杂结构，不需要通过复杂推导求解，能够利用现有傅里叶级数展开方法的结果进行耦合即可获得与傅里叶级数方法同等精度的复杂结构振动解析解的基于傅里叶级数的复杂结构振动解析分析方法。

为了解决上述技术问题，本发明采用的技术方案为：

一种基于傅里叶级数的复杂结构振动解析分析方法，其特征在于步骤包括：

1)将被分析复杂结构进行子结构划分为规则的子结构；

2)根据子结构的形状为各子结构设定控制微分方程的位移傅里叶级数；

3)将所述子结构控制微分方程的位移傅里叶级数代入各子结构的控制微分方程，采用傅里叶级数展开方法获得各子结构包含质量、刚度矩阵的振动方程；

4)根据各子结构包含质量、刚度矩阵的振动方程求取各子结构的联接点至激励点或响应点的位移阻抗函数结果；

5)依据各子结构激励点与耦合点之间的阻抗关系、耦合点处的力平衡和位移相容条件，推导得到整个复杂结构系统关于子结构位移阻抗函数的激励点与响应点之间的位移方程，将所述各子结构位移阻抗函数结果代入所述位移方程，求得被分析复杂结构的位移响应解。

优选地，所述步骤1)进行子结构划分时，所述规则的子结构的边界为固支边界、简支边界、自由边界、弹性边界中的一种。

优选地，所述步骤3)的详细步骤包括：

3.1)将所述子结构控制微分方程的位移傅里叶级数代入各子结构的控制微分方程；

3.2)将各子结构代入位移傅里叶级数后的控制微分方程的两边同时乘以预设的因子，并在子结构对应的范围进行积分，得到各子结构包含质量、刚度矩阵的振动方程。

优选地，所述步骤4)的具体步骤是指：根据所述包含质量、刚度矩阵的振动方程求取各子结构的位移傅里叶级数的未知幅值系数，将所述未知幅值系数代入各子结构的位移傅里叶级数、并令施加的外力幅值为单位1进行求解，求解得到的结果即为各子结构的联接点至激励点或响应点的位移阻抗函数结果。

本发明基于傅里叶级数的复杂结构振动解析分析方法具有下述优点：本发明将被分析复杂结构进行子结构划分为规则的子结构，根据子结构的形状为各子结构设定控制微分方程的位移傅里叶级数，将子结构控制微分方程的位移傅里叶级数代入各子结构的控制微分方程，采用傅里叶级数展开方法获得各子结构包含质量、刚度矩阵的振动方程，根据各子结构包含质量、刚度矩阵的振动方程求取各子结构的联接点至激励点或响应点的位移阻抗函数结果，依据各子结构激励点与耦合点之间的阻抗关系、耦合点处的力平衡和位移相容条件，推导得到整个复杂结构系统关于子结构位移阻抗函数的激励点与响应点之间的位移方程，将各子结构位移阻抗函数结果代入位移方程，求得被分析复杂结构的位移响应解，本发明能够发挥傅里叶级数展开方法作为解析方法的收敛快速、计算所需资源少的优点，通用性好，能够推广到各种可规则化子结构的复杂结构，不需要通过复杂推导求解，能够利用现有傅里叶级数展开方法的结果进行耦合即可获得与傅里叶级数方法同等精度的复杂结构振动解析解。

附图说明

图1为本发明实施例中被分析复杂结构的示意图。

图2为本发明实施例方法的基本流程示意图。

图3为本发明实施例中梁-圆柱壳耦合结构的抽象结构示意图。

图4为本发明实施例中两个子结构的耦合条件示意图。

图5为本发明实施例和现有技术关于子结构A的跨点阻抗结果对比曲线图。

图6为本发明实施例和现有技术关于子结构A的原点阻抗结果对比曲线图。

图7为本发明实施例和现有技术关于子结构B的跨点阻抗结果对比曲线图。

图8为本发明实施例和现有技术关于子结构B的原点阻抗结果对比曲线图。

图9为本发明实施例和现有技术的激励点位移响应对比曲线图。

图10为本发明实施例和现有技术的耦合点位移响应对比曲线图。

具体实施方式

下文将以如图1所示被分析复杂结构为例，对本发明基于傅里叶级数的复杂结构振动解析分析方法进行进一步的详细说明。

参见图1，本实施例中的被分析复杂结构为梁-圆柱壳耦合结构，梁结构通过左右两端点处的弹性支撑k_l、k_r与圆柱壳结构耦合，梁结构上作用有幅值为F、圆频率ω为的简谐激励Fe^jwt。本实施例中，梁的结构参数为：横截面积A＝3.1416×10^-4m²，长度L₁＝1m，两端弹性支撑刚度取为k_l＝k_r＝k＝1×10¹⁰N/m；圆柱壳的结构参数为：半径R＝0.25m，长度L₂＝2m，壁厚h＝0.005m，杨氏模量E＝2.1×10¹¹N/m²，密度ρ＝7850kg/m³，泊松比μ＝0.3，边界条件为两端简支。耦合点在梁两端端点，在圆柱壳长度角度方向的(0.5m,0)和(1.5m,0)处，为对称结构。在梁结构中点施加的简谐激励幅值F＝1N。

如图2所示，本实施例基于傅里叶级数的复杂结构振动解析分析方法的步骤包括：

1)将被分析复杂结构进行子结构划分为规则的子结构。

本实施例步骤1)进行子结构划分时，规则的子结构的边界为固支边界、简支边界、自由边界、弹性边界中的一种，规则的子结构可为集中质量点、梁、板、壳等特征。如图3和图4所示，本实施例中梁-圆柱壳耦合结构在进行子结构划分后，被划分为梁(子结构A)和两端简支的圆柱壳(子结构B)两个规则的子结构，两个子结构通过耦合点2、耦合点3联接，且在子结构A的1点上有一简谐激励F₁e^jwt。参见图4，梁(子结构A)和两端简支的圆柱壳(子结构B)两个规则的子结构通过力F和位移X在耦合点联接，其中F₂ ^A表示作用于子结构A、耦合点2上的力，F₂ ^B表示作用于子结构B、耦合点2上的力，F₃ ^A表示作用于子结构A、耦合点3上的力，F₃ ^B表示作用于子结构B、耦合点3上的力；表示作用于子结构A、耦合点2上的位移，表示作用于子结构B、耦合点2上的位移，表示作用于子结构A、耦合点3上的位移，表示作用于子结构B、耦合点3上的位移。

2)根据子结构的形状为各子结构设定控制微分方程的位移傅里叶级数。

本实施例中，子结构A(梁)的控制微分方程的位移傅里叶级数的展开形式可表示为式(1)；

$w\left(x\right)=\underset{m=0}{\overset{\mathrm{\∞}}{\Σ}}{A}_m({\mathrm{cos\λ}}_{m}x+\mathrm{\ζ}{\left(x\right)}^T{\mathrm{HQ}}^m)+\underset{m=0}{\overset{\mathrm{\∞}}{\Σ}}{A}_m\mathrm{\ζ}{\left(x\right)}^T{\mathrm{HQ}}^m---\left(1\right)$

式(1)中，w(x)为子结构A(梁)的振动横向位移，A_m为收敛因子取m时的未知幅值系数，x为局部坐标，x的起点为梁的最左端点，m为收敛因子；其中，λ_m的计算式如式(1-1)所示，ζ(x)^T的计算式如式(1-2)所示，H的计算式如式(1-3)所示，Q^m的计算式如式(1-4)所示；

λ_m＝mπ/L₁,(m＝0,1,2…)(1-1)

$\mathrm{\ζ}{\left(x\right)}^T=\left\{\begin{array}{c}-(15x^4-60L_1{x}^3+60{L_1}^2{x}^2-8{L_1}^4)/360L_1\\ (15x^4-30{L_1}^2{x}^2+7{L_1}^4)/360L_1\\ (60L_{1}x-2{L_1}^2-3x^2)/6L_1\\ (3x^2-{L_1}^2)/6L_1\end{array}\right\}---(1-2)$

$H=\left[\begin{array}{cccc}\frac{8k_l{L_1}^3}{360}+1& \frac{7k_l{L_1}^3}{360}& \frac{-k_l{L}_1}{3}& \frac{-k_l{L}_1}{6}\\ \frac{7k_r{L_1}^3}{360}& \frac{8k_r{L_1}^3}{360}+1& \frac{-k_r{L}_1}{6}& \frac{-k_r{L}_1}{3}\\ \frac{L_1}{3}& \frac{L_1}{6}& \frac{1}{L_1}& \frac{-1}{L_1}\\ \frac{L_1}{6}& \frac{L_1}{3}& \frac{-1}{L_1}& \frac{1}{L_1}\end{array}\right]---(1-3)$

$Q^m={\left\{\begin{array}{cccc}-k_l& {(-1)}^m{k}_r& -{\mathrm{\λ}}_m^2& {(-1)}^m{\mathrm{\λ}}_m^2\end{array}\right\}}^T---(1-4)$

式(1-1)～(1-4)中，m为收敛因子，L₁为子结构A(梁)的长度，x为局部坐标，k_l和k_r分别为子结构A(梁)两端的弹性支撑刚度(1×10¹⁰N/m)。

本实施例中，子结构B(圆柱壳)三个方向的控制微分方程的位移傅里叶级数的展开形式可表示为式(2)～(4)所示函数式；

$u(s,\mathrm{\θ},t)=e^{j\mathrm{\ω}t}\underset{m=1}{\overset{M}{\Σ}}\underset{n=0}{\overset{N}{\Σ}}{U}_{mn}{\mathrm{cos\λ}}_{m}scos\left(n\mathrm{\θ}\right)---\left(2\right)$

$v(s,\mathrm{\θ},t)=e^{j\mathrm{\ω}t}\underset{m=1}{\overset{M}{\Σ}}\underset{n=0}{\overset{N}{\Σ}}{V}_{mn}{\mathrm{sin\λ}}_{m}ssin\left(n\mathrm{\θ}\right)---\left(3\right)$

$w(s,\mathrm{\θ},t)=e^{j\mathrm{\ω}t}\underset{m=1}{\overset{M}{\Σ}}\underset{n=0}{\overset{N}{\Σ}}{W}_{mn}{\mathrm{sin\λ}}_{m}scos\left(n\mathrm{\θ}\right)---\left(4\right)$

式(2)～(4)中，u(s,θ,t)、v(s,θ,t)、w(s,θ,t)分别为子结构B(圆柱壳)在轴向、切向、径向的位移，(s,θ)为子结构B(圆柱壳)在柱坐标系统中的点坐标，s起点为子结构B(圆柱壳)的最左端位置，ω为圆频率，U_mn、V_mn、W_mn分别为子结构B(圆柱壳)在轴向、切向、径向三个方向位移的未知幅值系数，λ_m＝mπR/L₂，m为收敛因子，n为子结构B(圆柱壳)的周向模态数，R为子结构B(圆柱壳)的半径，L₂为子结构B(圆柱壳)的长度。

3)将各子结构控制微分方程的位移傅里叶级数代入子结构的控制微分方程，采用傅里叶级数展开方法获得各子结构包含质量、刚度矩阵的振动方程。

本实施例中，步骤3)的详细步骤包括：

3.1)将子结构控制微分方程的位移傅里叶级数代入各子结构的控制微分方程；

3.2)将各子结构代入位移傅里叶级数后的控制微分方程的两边同时乘以预设的因子，并在子结构对应的范围进行积分，得到各子结构包含质量、刚度矩阵的振动方程。

本实施例中，子结构A(梁)的控制微分方程如式(5)所示；

$EI\frac{d^{4}w\left(x\right)}{{\mathrm{dx}}^4}-{\mathrm{\ρA\ω}}^{2}w\left(x\right)-{\mathrm{\ω}}^2\underset{k=1}{\overset{N_m}{\mathrm{\Σ}}}{m}_k\mathrm{\δ}(x-x_k)w\left(x\right)=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{F}_j\mathrm{\δ}(x-{x_j}^f)---\left(5\right)$

式(5)中，EI为子结构A(梁)的弯曲刚度(抗弯刚度)，ρA为子结构A(梁)的单位长度质量，w(x)为子结构A(梁)的振动横向位移，ω为圆频率，m_k为子结构A(梁)上第k个集中质量点的质量，x为局部坐标，x_k为子结构A(梁)上第k个集中质量点的位置，N_m为子结构A(梁)上集中质量点个数，F_j为子结构A(梁)上第j个简谐外力幅值，x_j ^f为子结构A(梁)上第j个简谐外力位置，J为简谐外力个数，δ为狄拉克delta函数。

采用傅里叶级数展开方法获得子结构A(梁)包含质量、刚度矩阵的振动方程时，具体是指将子结构A(梁)的控制微分方程的位移傅里叶级数代入子结构A(梁)的控制微分方程，即将式(1)代入式(5)，再将控制微分方程的两边同时乘以预设的因子2/L₁[cosλ_mx+ζ(x)^THQ^m]，再在0到L₁上积分，得到子结构A(梁)包含质量、刚度矩阵的振动方程如式(6)所示；

([K]-ω²[M]){A}＝{F}(6)

式(6)中，ω为圆频率，矩阵K的形式如式(6-1)所示，M的形式如式(6-2)所示，矩阵A为傅里叶级数的未知幅值系数矩阵，矩阵A的形式如式(6-3)所示，矩阵F的形式如式(6-4)所示。

$K=\left[\begin{array}{cccccc}{c}^T{H}^{-1}{Q}^0& c^T{H}^{-1}{Q}^1& ...& c^T{H}^{-1}{Q}^{m^{\′}}& ...& c^T{H}^{-1}{Q}^M\\ 0& {\mathrm{\λ}}_1^4& ...& 0& ...& 0\\ \·& \·& & \·& & \·\\ \·& \·& ...& \·& ...& \·\\ \·& \·& & \·& & \·\\ 0& 0& ...& {\mathrm{\λ}}_m^4{\mathrm{\δ}}_{{\mathrm{mm}}^{\′}}& ...& 0\\ \·& \·& & \·& & \·\\ \·& \·& ...& \·& ...& \·\\ \·& \·& & \·& & \·\\ 0& 0& ...& 0& ...& {\mathrm{\λ}}_M^4\end{array}\right]---(6-1)$

$M=\left[\begin{array}{cccccc}1& 0& ...& 0& ...& 0\\ S_{10}& 1+S_1& ...& S_{1m^{\′}}& ...& S_{1M}\\ \·& \·& & \·& & \·\\ \·& \·& ...& \·& ...& \·\\ \·& \·& & \·& & \·\\ S_{m0}& S_{m1}& ...& {\mathrm{\δ}}_{{\mathrm{mm}}^{\′}}+S_{{\mathrm{mm}}^{\′}}& ...& S_{mM}\\ \·& \·& & \·& & \·\\ \·& \·& ...& \·& ...& \·\\ \·& \·& & \·& & \·\\ S_{M0}& S_{M1}& ...& S_{{\mathrm{Mm}}^{\′}}& ...& 1+S_{MM}\end{array}\right]---(6-2)$

A＝{A₁…A_m…A_M}^T(6-3)

F＝{f₁…f_m…f_M}^T(6-4)

式(6-1)～(6-4)中，c＝1/L₁{-1100}^T，其中L₁为子结构A(梁)的长度，H的计算式如式(1-3)所示，Q^m′的计算式如式(1-4)所示，δ_mm'为克罗内克delta函数；S_mm'的计算式为 $S_{{\mathrm{mm}}^{\′}}=P_m^T{H}^{-1}{Q}^{m^{\′}}$ ，其中 $P_m=2/L_1{\left\{\begin{array}{cccc}1/{\mathrm{\λ}}_m^4& {(-1)}^{m+1}/{\mathrm{\λ}}_m^4& -1/{\mathrm{\λ}}_m^2& {(-1)}^m/{\mathrm{\λ}}_m^2\end{array}\right\}}^T,$ λ_m＝mπ/L₁,(m＝0,1,2…)，m为收敛因子，L₁为子结构A(梁)的长度，H的计算式如式(1-3)所示，Q^m'为的计算式如式(1-4)所示；A_m为收敛因子取m时的未知幅值系数；f_m的表达式如式(6-5)所示；

$f_m=\frac{2}{L_{1}EI}{F}_0({\mathrm{cos\λ}}_m{x_j}^f+\mathrm{\ζ}{\left({x_j}^f\right)}^T{H}^T{Q}^m)---(6-5)$

式(6-5)中，L₁为子结构A(梁)的长度，EI为子结构A(梁)的弯曲刚度(抗弯刚度)，F₀为激励幅值系数，x_j ^f为子结构A(梁)上第j个简谐外力位置，ζ(x_j ^f)的计算式如式(1-2)所示，H的计算式如式(1-3)所示，Q^m的计算式如式(1-4)所示；λ_m＝mπ/L₁,(m＝0,1,2…)，m为收敛因子，L₁为子结构A(梁)的长度；

本实施例中，子结构B(圆柱壳)的控制微分方程如式(7)～(9)所示；

$\frac{{\∂}^{2}u}{\∂s^2}+\frac{(1-\mathrm{\μ})}{2}(1+{\mathrm{\β}}^2)\frac{{\∂}^{2}u}{\∂{\mathrm{\θ}}^2}+\frac{(1+\mathrm{\μ})}{2}\frac{{\∂}^{2}v}{\∂s\∂\mathrm{\θ}}+\mathrm{\μ}\frac{\∂w}{\∂s}-{\mathrm{\β}}^2\frac{{\∂}^{3}w}{\∂s^3}+{\mathrm{\β}}^2\frac{(1-\mathrm{\μ})}{2}\frac{{\∂}^{3}w}{\∂s\∂{\mathrm{\θ}}^2}-\frac{1}{{c_L}^2}\frac{{\∂}^{2}u}{\∂t^2}=0---\left(7\right)$

$\frac{(1+\mathrm{\μ})}{2}\frac{{\∂}^{2}u}{\∂s\∂\mathrm{\θ}}+\frac{(1+\mathrm{\μ})}{2}\frac{{\∂}^{2}v}{\∂s^2}+\frac{{\∂}^{2}v}{\∂{\mathrm{\θ}}^2}+\frac{\∂w}{\∂\mathrm{\θ}}+{\mathrm{\β}}^2(\frac{3(1-\mathrm{\μ})}{2}\frac{{\∂}^{2}v}{\∂s^2}-\frac{(3-\mathrm{\μ})}{2}\frac{{\∂}^{3}w}{\∂s^2\∂\mathrm{\θ}})-\frac{1}{{c_L}^2}\frac{{\∂}^2}{\∂t^2}=0---\left(8\right)$

$\begin{array}{c}{\mathrm{\β}}^2(\frac{{\∂}^{4}w}{\∂s^4}+2\frac{{\∂}^{4}w}{\∂s^2\∂{\mathrm{\θ}}^2}+\frac{{\∂}^{4}w}{\∂{\mathrm{\θ}}^2}-\frac{{\∂}^{4}u}{\∂s^3}+\frac{(1-\mathrm{\μ})}{2}\frac{{\∂}^{3}u}{\∂s\∂{\mathrm{\θ}}^2}-\frac{(3-\mathrm{\μ})}{2}\frac{{\∂}^{2}v}{\∂s^2\∂\mathrm{\θ}}+2\frac{{\∂}^{2}w}{\∂{\mathrm{\θ}}^2})+\mathrm{\μ}\frac{\∂u}{\∂s}+\frac{\∂v}{\∂\mathrm{\θ}}\\ +w(1+{\mathrm{\β}}^2)+\frac{1}{{c_L}^2}\frac{{\∂}^{2}w}{\∂t^2}=\frac{(1-{\mathrm{\μ}}^2)}{Eh}{p}_r\end{array}---\left(9\right)$

式(7)～(9)中，u、v、w分别为子结构B(圆柱壳)在轴向、切向、径向三个方向的位移，(s,θ)为子结构B(圆柱壳)在柱坐标系统中的点坐标，s起点为子结构B(圆柱壳)的最左端位置，μ为泊松比，β为子结构B(圆柱壳)的壳厚因子，c_L为子结构B(圆柱壳)的壳体纵向波的传播速度。壳厚因子β的计算方式为其中h为子结构B(圆柱壳)的壁厚，R为子结构B(圆柱壳)的半径；传播速度c_L的计算表达式为c_L＝[E/ρ(1-μ²)]^1/2，其中ρ为密度，E为杨氏模量，p_r为径向分布力。

采用傅里叶级数展开方法获得子结构B(圆柱壳)包含质量、刚度矩阵的振动方程时，具体是指将子结构B(圆柱壳)的控制微分方程的位移傅里叶级数代入子结构B(圆柱壳)的控制微分方程，即将式(2)～(4)代入式(7)～(9)，再将式(7)的控制微分方程的两边同时乘以预设的因子cosλ_mscos(nθ)、将式(8)的控制微分方程的两边同时乘以sinλ_mssin(nθ)、将式(9)的控制微分方程的两边同时乘以sinλ_mscos(nθ)，再在在0～2π及0～L₂/R上积分，得到子结构B(圆柱壳)包含质量、刚度矩阵的振动方程如式(10)所示；

$\{\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ A_{21}& A_{22}& A_{23}\\ A_{31}& A_{32}& A_{33}\end{array}\right]-{\mathrm{\ω}}^2{m}_s\left[\begin{array}{ccc}I& 0& 0\\ 0& I& 0\\ 0& 0& I\end{array}\right]\}\left\{\begin{array}{c}{U}_{mn}\\ V_{mn}\\ W_{mn}\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ p_{r,mn}\end{array}\right\}---\left(10\right)$

式(10)中，A₁₁～A₁₁的表达式如式(10-1)～(10-9)所示，ω为圆频率，m_s＝2πRhL₂ρ/4，其中R为子结构B(圆柱壳)的半径，h为子结构B(圆柱壳)的壁厚，L₂为子结构B(圆柱壳)的长度，ρ为密度，矩阵I为对应维度的单位矩阵，U_mn、V_mn、W_mn分别为三个方向位移的未知幅值系数，p_r,mn为径向分布力元素，p_r,mn的表达式如式(10-10)所示；

$A_{11}={\mathrm{\δ}}_{{\mathrm{mm}}^{\′}}{\mathrm{\δ}}_{{\mathrm{nn}}^{\′}}(1+{\mathrm{\δ}}_{n0}){c_L}^2{m}_s\[{\left({\mathrm{\λ}}_m\right)}^2+\frac{(1-\mathrm{\μ})(1+{\mathrm{\β}}^2)}{2}{n}^2\]---(10-1)$

$A_{12}=-{\mathrm{\δ}}_{{\mathrm{mm}}^{\′}}{\mathrm{\δ}}_{{\mathrm{nn}}^{\′}}(1+{\mathrm{\δ}}_{n0}){c_L}^2{m}_s\[\frac{(1+\mathrm{\μ})}{2}{\mathrm{n\λ}}_m\]---(10-2)$

$A_{13}=-{\mathrm{\δ}}_{{\mathrm{mm}}^{\′}}{\mathrm{\δ}}_{{\mathrm{nn}}^{\′}}(1+{\mathrm{\δ}}_{n0}){c_L}^2{m}_s\[{\mathrm{\μ\λ}}_m+{\mathrm{\β}}^2{{\mathrm{\λ}}_m}^3-\frac{{\mathrm{\β}}^2(1-\mathrm{\μ})}{2}{\mathrm{\λ}}_m{n}^2\]---(10-3)$

A₂₁＝A₁₂(10-4)

$A_{22}={\mathrm{\δ}}_{{\mathrm{mm}}^{\′}}{\mathrm{\δ}}_{{\mathrm{nn}}^{\′}}(1+{\mathrm{\δ}}_{n0}){c_L}^2{m}_s\[\frac{(1-\mathrm{\μ})}{2}{{\mathrm{\λ}}_m}^2+n^2+\frac{3{\mathrm{\β}}^2(1-\mathrm{\μ})}{2}{{\mathrm{\λ}}_m}^2\]---(10-5)$

$A_{23}={\mathrm{\δ}}_{{\mathrm{mm}}^{\′}}{\mathrm{\δ}}_{{\mathrm{nn}}^{\′}}(1+{\mathrm{\δ}}_{n0}){c_L}^2{m}_s\[n+\frac{{\mathrm{\β}}^2(3-\mathrm{\μ})}{2}{{\mathrm{n\λ}}_m}^2\]---(10-6)$

A₃₁＝A₁₃(10-7)

A₃₂＝A₂₃(10-8)

A₃₃＝δ_mm'δ_nn'(1+δ_n0)c_L ²m_s[β²(λ_m ⁴+2λ_m ²n²+n⁴-2n²)+1+β²](10-9)

式(10-1)～(10-9)中，δ_mm'、δ_nn'、δ_n0均为克罗内克delta函数，μ为泊松比，β为子结构B(圆柱壳)的壳厚因子，n为子结构B(圆柱壳)的周向模态数；β为子结构B(圆柱壳)的壳厚因子，c_L为子结构B(圆柱壳)的壳体纵向波的传播速度；m_s＝2πRhL₂ρ/4，其中R为子结构B(圆柱壳)的半径，h为子结构B(圆柱壳)的壁厚，L₂为子结构B(圆柱壳)的长度，ρ为密度。λ_m＝mπR/L₂，m为收敛因子，R为子结构B(圆柱壳)的半径，L₂为子结构B(圆柱壳)的长度；壳厚因子β的计算方式为其中h为子结构B(圆柱壳)的壁厚，R为子结构B(圆柱壳)的半径；传播速度c_L的计算表达式为c_L＝[E/ρ(1-μ²)]^1/2，其中ρ为密度，E为杨氏模量，p_r为径向分布力。

$p_{r,mn}={\∫}_0^{L_2/R}{\∫}_0^{2\mathrm{\π}}{p}_r(s,\mathrm{\θ})\sin \left({\mathrm{\λ}}_{m}s\right)cos\left(n\mathrm{\θ}\right)dsd\mathrm{\θ}---(10-10)$

式(10-10)中，p_r(s,θ)为作用于点(s,θ)的径向激励力，(s,θ)为子结构B(圆柱壳)在柱坐标系统中的点坐标，s起点为子结构B(圆柱壳)的最左端位置，λ_m＝mπR/L₂，m为收敛因子，R为子结构B(圆柱壳)的半径，L₂为子结构B(圆柱壳)的长度。当外激励为径向点激励时，则有p_r(s,θ)＝f_rδ(s-s₀)δ(θ-θ₀)，其中f_r为激励力，(s₀,θ₀)为激励力位置，δ为DiracDelta函数，此时式(10-10)可转换为式(10-11)所示函数式；

p_r,mn＝f_rsin(λ_ms₀)cos(nθ₀)(10-11)

式(10-11)中，p_r,mn为径向分布力元素，f_r为激励力，(s₀,θ₀)为激励力位置，n为子结构B(圆柱壳)的周向模态数，λ_m＝mπR/L₂，m为收敛因子，R为子结构B(圆柱壳)的半径，L₂为子结构B(圆柱壳)的长度。

4)根据各子结构包含质量、刚度矩阵的振动方程求取各子结构的联接点至激励点或响应点的位移阻抗函数结果。

本实施例中，步骤4)的具体步骤是指：根据包含质量、刚度矩阵的振动方程求取各子结构的位移傅里叶级数的未知幅值系数，将未知幅值系数代入各子结构的位移傅里叶级数、并令施加的外力幅值为单位1进行求解，求解得到的结果即为各子结构的联接点至激励点或响应点的位移阻抗函数结果。

令施加的外力幅值为单位1，将式(6)所示子结构A(梁)包含质量、刚度矩阵的振动方程的傅里叶级数的未知幅值系数矩阵A代入式(1)所示子结构A(梁)的控制微分方程的位移傅里叶级数就可以求得子结构A(梁)的位移结果，当施加的外力幅值为单位1时，所求得的位移结果即可视为子结构A(梁)的阻抗函数结果。因此，对上述参数下的子结构A(梁)进行位移阻抗函数结果求取，收敛因子取m＝20，可以求得子结构A(梁)的激励点(梁结构中点)、耦合点(梁两个端点)等的原点阻抗函数跨点阻抗函数结果，其左右两个端点(耦合点)之间的原点阻抗函数(由于对称结构：跨点阻抗函数(同理：将上述结果与有限元模型(FEA：ANSYS)谐响应分析进行对比，所得跨点阻抗结果对比曲线图如图5所示，所得原点阻抗结果对比曲线图如图6所示。由图5、6可知，本实施例对于子结构A(梁)的阻抗函数处理是正确的。

通过式(10)可以求解在径向点激励下子结构B(圆柱壳)三个方向的未知幅值系数{…U_mn…V_mn…W_mn…}，代回式(2)～(4)即可求得子结构B(圆柱壳)三个方向的位移结果，当施加的外力幅值为单位1时，所求得的位移结果即可视为子结构B(圆柱壳)的阻抗函数结果。对上述参数下的子结构B(圆柱壳)进行位移阻抗函数结果求取，收敛因子m及子结构B(圆柱壳)的周向模态数n均取20(m＝n＝20)，求得子结构B(圆柱壳)在两个耦合点(s_l，θ_l)＝(0.5，0)、(s_r，θ_r)＝(1.5，0)处的原点阻抗函数(由于对称结构：)、跨点阻抗函数(同理：)结果与有限元模型(FEA：ANSYS)谐响应分析进行对比，所得跨点阻抗结果对比曲线图如图7所示，所得原点阻抗结果对比曲线图如图8所示。由图7、8可知，本实施例对于子结构B(圆柱壳)的阻抗函数处理是正确的。

5)依据各子结构激励点与耦合点之间的阻抗关系、耦合点处的力平衡和位移相容条件，推导得到整个复杂结构系统关于子结构位移阻抗函数的激励点与响应点之间的位移方程，将各子结构位移阻抗函数结果代入位移方程，求得被分析复杂结构的位移响应解。

本实施例中，子结构A(梁)的激励点与耦合点之间的阻抗关系如式(11)或式(12)所示；

$\left\{\begin{array}{c}{X}_1^A=H_{11}^A{F}_1^A+H_{12}^A{F}_2^A+H_{13}^A{F}_3^A\\ X_2^A=H_{21}^A{F}_1^A+H_{22}^A{F}_2^A+H_{23}^A{F}_3^A\\ X_3^A=H_{31}^A{F}_1^A+H_{32}^A{F}_2^A+H_{33}^A{F}_3^A\end{array}\right\}---\left(11\right)$

$\left\{\begin{array}{c}{X}_1^A\\ X_2^A\\ X_3^A\end{array}\right\}=\left\{\begin{array}{ccc}{H}_{11}^A& H_{12}^A& H_{13}^A\\ H_{21}^A& H_{22}^A& H_{23}^A\\ H_{31}^A& H_{32}^A& H_{33}^A\end{array}\right\}\left\{\begin{array}{c}{F}_1^A\\ F_2^A\\ F_3^A\end{array}\right\}---\left(12\right)$

式(11)和式(12)中，分别为子结构A(梁)在激励点1、耦合点2和3处的位移，F₁ ^A、F₂ ^A、F₃ ^A分别为子结构A(梁)在激励点1、耦合点2和3处的力， 分别为子结构A(梁)在激励点1、耦合点2和3处的原点阻抗函数， 分别为子结构A(梁)在激励点1、耦合点2和3三者之间的跨点阻抗函数，例如表示子结构A(梁)在激励点1、耦合点2之间的跨点阻抗函数，依次类推，在此不再赘述。

本实施例中，子结构B(圆柱壳)的激励点与耦合点之间的阻抗关系如式(13)或式(14)所示；

$\left\{\begin{array}{c}{X}_2^B=H_{22}^B{F}_2^B+H_{23}^B{F}_3^B\\ X_3^B=H_{32}^B{F}_2^B+H_{33}^B{F}_3^B\end{array}\right\}---\left(13\right)$

$\left\{\begin{array}{c}{X}_2^B\\ X_3^B\end{array}\right\}=\left\{\begin{array}{cc}{H}_{22}^B& H_{23}^B\\ H_{32}^B& H_{33}^B\end{array}\right\}\left\{\begin{array}{c}{F}_2^B\\ F_3^B\end{array}\right\}---\left(14\right)$

式(13)和式(14)中，分别为子结构B(圆柱壳)B在耦合点2和3处的位移，F₂ ^B、F₃ ^B分别为子结构B(圆柱壳)B在耦合点2和3处的力，分别为子结构B(圆柱壳)B在耦合点2和3处的原点阻抗函数，分别为子结构B(圆柱壳)B在耦合点2和3之间的跨点阻抗函数。

参见图4可知，本实施例中耦合点处的力平衡和位移相容条件(耦合条件)为：F₂ ^A＝-F₂ ^B，F₃ ^A＝-F₃ ^B，，其中各参数符号的符号和式(11)～式(14)中的含义相同。根据耦合点处的力平衡和位移相容条件，子结构B(圆柱壳)的激励点与耦合点之间的阻抗关系式(13)和式(14)可以推导为式(15)和式(16)；

$\left\{\begin{array}{c}{X}_2^A\\ X_3^A\end{array}\right\}=\left\{\begin{array}{c}{X}_2^B\\ X_3^B\end{array}\right\}=\left\{\begin{array}{cc}{H}_{22}^B& H_{23}^B\\ H_{32}^B& H_{33}^B\end{array}\right\}\left\{\begin{array}{c}{F}_2^B\\ F_3^B\end{array}\right\}=-\left\{\begin{array}{cc}{H}_{22}^B& H_{23}^B\\ H_{32}^B& H_{33}^B\end{array}\right\}\left\{\begin{array}{c}{F}_2^A\\ F_3^A\end{array}\right\}---\left(15\right)$

$-\left\{\begin{array}{cc}{\mathrm{\χ}}_{22}& {\mathrm{\χ}}_{23}\\ {\mathrm{\χ}}_{32}& {\mathrm{\χ}}_{33}\end{array}\right\}\left\{\begin{array}{c}{X}_2^A\\ X_3^A\end{array}\right\}=\left\{\begin{array}{c}{F}_2^A\\ F_3^A\end{array}\right\}---\left(16\right)$

式(15)和式(16)中，F₂ ^A表示作用于子结构A、耦合点2上的力，F₂ ^B表示作用于子结构B、耦合点2上的力，F₃ ^A表示作用于子结构A、耦合点3上的力，F₃ ^B表示作用于子结构B、耦合点3上的力；表示作用于子结构A、耦合点2上的位移，表示作用于子结构B、耦合点2上的位移，表示作用于子结构A、耦合点3上的位移，表示作用于子结构B、耦合点3上的位移；分别为子结构B(圆柱壳)B在耦合点2和3处的原点阻抗函数；分别为子结构B(圆柱壳)B在耦合点2和3之间的跨点阻抗函数；其中，式(16)中的 $\left\{\begin{array}{cc}{\mathrm{\χ}}_{22}& {\mathrm{\χ}}_{23}\\ {\mathrm{\χ}}_{32}& {\mathrm{\χ}}_{33}\end{array}\right\}$ 的函数表达式如式(17)所示；

$\left\{\begin{array}{cc}{\mathrm{\χ}}_{22}& {\mathrm{\χ}}_{23}\\ {\mathrm{\χ}}_{32}& {\mathrm{\χ}}_{33}\end{array}\right\}={\left\{\begin{array}{cc}{H}_{22}^B& H_{23}^B\\ H_{32}^B& H_{33}^B\end{array}\right\}}^{-1}---\left(17\right)$

式(17)中，分别为子结构B(圆柱壳)B在耦合点2和3处的原点阻抗函数；分别为子结构B(圆柱壳)B在耦合点2和3之间的跨点阻抗函数。

将式(15)或式(16)代入子结构A(梁)的阻抗关系式式(11)或式(12)中，可以重新推导得到式(18)；可以最终求解得到式(19)～式(21)。

$\left\{\begin{array}{c}{X}_1^A\\ X_2^A\\ X_3^A\end{array}\right\}=\left\{\begin{array}{ccc}{H}_{11}^A& H_{12}^A& H_{13}^A\\ H_{21}^A& H_{22}^A& H_{23}^A\\ H_{31}^A& H_{32}^A& H_{33}^A\end{array}\right\}\left\{\begin{array}{c}{F}_1^A\\ F_2^A\\ F_3^A\end{array}\right\}=\left\{\begin{array}{ccc}{H}_{11}^A& H_{12}^A& H_{13}^A\\ H_{21}^A& H_{22}^A& H_{23}^A\\ H_{31}^A& H_{32}^A& H_{33}^A\end{array}\right\}\left\{\begin{array}{c}{F}_1^A\\ -{\mathrm{\χ}}_{22}{X}_2^A-{\mathrm{\χ}}_{23}{X}_3^A\\ -{\mathrm{\χ}}_{32}{X}_2^A-{\mathrm{\χ}}_{33}{X}_3^A\end{array}\right\}---\left(18\right)$

式(18)中，各个参数符号与式(11)～式(17)相同。

$X_1^A=(H_{11}^A-(H_{12}^A{\mathrm{\χ}}_{22}+H_{13}^A{\mathrm{\χ}}_{32})\frac{H_{21}^A(1-A)-(H_{22}^A{\mathrm{\χ}}_{23}+H_{23}^A{\mathrm{\χ}}_{33})B}{(1+H_{22}^A{\mathrm{\χ}}_{22}+H_{23}^A{\mathrm{\χ}}_{32})(1-A)}-(H_{12}^A{\mathrm{\χ}}_{23}+H_{13}^A{\mathrm{\χ}}_{33}\left)\frac{B}{(1-A)}\right)F_1^A---\left(19\right)$

$X_2^A=\frac{H_{21}^A(1-A)-(H_{22}^A{\mathrm{\χ}}_{23}+H_{23}^A{\mathrm{\χ}}_{33})B}{(1+H_{22}^A{\mathrm{\χ}}_{22}+H_{23}^A{\mathrm{\χ}}_{32})(1-A)}{F}_1^A---\left(20\right)$

$X_3^A=\frac{B}{(1-A)}{F}_1^A---\left(21\right)$

式(19)～式(21)中，A的表达式如式(22)所示，B的表达式如式(23)所示，其余各个参数符号与式(11)～式(17)相同。

$A=\frac{(H_{32}^A{\mathrm{\χ}}_{22}+H_{33}^A{\mathrm{\χ}}_{32})(H_{22}^A{\mathrm{\χ}}_{23}+H_{23}^A{\mathrm{\χ}}_{33})}{(1+H_{22}^A{\mathrm{\χ}}_{22}+H_{23}^A{\mathrm{\χ}}_{32})}-(H_{32}^A{\mathrm{\χ}}_{23}+H_{33}^A{\mathrm{\χ}}_{33})---\left(22\right)$

$B=\frac{H_{31}^A(1+H_{22}^A{\mathrm{\χ}}_{22}+H_{23}^A{\mathrm{\χ}}_{32})-H_{21}^A(H_{32}^A{\mathrm{\χ}}_{22}+H_{33}^A{\mathrm{\χ}}_{32})}{(1+H_{22}^A{\mathrm{\χ}}_{22}+H_{23}^A{\mathrm{\χ}}_{32})}---\left(23\right)$

式(22)～式(23)中，各个参数符号与式(11)～式(17)相同。

最终求解得到的式(19)～式(21)即为耦合结构响应点与激振点之间关于子结构阻抗函数的位移关系式。在此基础上，将上文已经求得的子结构A(梁)原点阻抗函数( )、跨点阻抗函数()结果，子结构B(圆柱壳)原点阻抗函数()、跨点阻抗函数()结果代入式(19)～(21)，即可求得各子结构的对应阻抗函数结果即可求得耦合结构的位移响应解。将本实施例求得的耦合结构的位移响应解和现有技术的有限元模型(FEA：ANSYS)所得结果进行对比，得到的对比曲线图如图9和图10所示，其中图9为激励点位移响应的对比曲线图，图10为耦合点位移响应的对比曲线图。通过图9、10对比，验证了本实施例方法的正确性。

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

Claims (4)

1.一种基于傅里叶级数的复杂结构振动解析分析方法，其特征在于步骤包括：

1)将被分析复杂结构进行子结构划分为规则的子结构；

2)根据子结构的形状为各子结构设定控制微分方程的位移傅里叶级数；

3)将所述子结构控制微分方程的位移傅里叶级数代入各子结构的控制微分方程，采用傅里叶级数展开方法获得各子结构包含质量、刚度矩阵的振动方程；

4)根据各子结构包含质量、刚度矩阵的振动方程求取各子结构的联接点至激励点或响应点的位移阻抗函数结果；

5)依据各子结构激励点与耦合点之间的阻抗关系、耦合点处的力平衡和位移相容条件，推导得到整个复杂结构系统关于子结构位移阻抗函数的激励点与响应点之间的位移方程，将所述各子结构位移阻抗函数结果代入所述位移方程，求得被分析复杂结构的位移响应解。

2.根据权利要求1所述的基于傅里叶级数的复杂结构振动解析分析方法，其特征在于，所述步骤1)进行子结构划分时，所述规则的子结构的边界为固支边界、简支边界、自由边界、弹性边界中的一种。

3.根据权利要求2所述的基于傅里叶级数的复杂结构振动解析分析方法，其特征在于，所述步骤3)的详细步骤包括：

3.1)将所述子结构控制微分方程的位移傅里叶级数代入各子结构的控制微分方程；

3.2)将各子结构代入位移傅里叶级数后的控制微分方程的两边同时乘以预设的因子，并在子结构对应的范围进行积分，得到各子结构包含质量、刚度矩阵的振动方程。

4.根据权利要求3所述的基于傅里叶级数的复杂结构振动解析分析方法，其特征在于，所述步骤4)的具体步骤是指：根据所述包含质量、刚度矩阵的振动方程求取各子结构的位移傅里叶级数的未知幅值系数，将所述未知幅值系数代入各子结构的位移傅里叶级数、并令施加的外力幅值为单位1进行求解，求解得到的结果即为各子结构的联接点至激励点或响应点的位移阻抗函数结果。