一种基于谱有限元的矩形板振动模态计算方法

摘要

本发明公开了一种基于谱有限元的矩形板振动模态计算方法，包括以下步骤：基于能量函数变分原理得到板振动的边界值问题模型，再基于变量分离方法得到两个方向的自由振动方程；确定矩形板两个方向的边界条件，在预设定的模态阶数下采用不同方向相互迭代来计算振动频率和振型，直至两个方向计算得到的振动频率误差在某个指定范围时停止；利用最终得到的两个方向的振动模态进行叠加得到矩形板振动模态。本发明针对多种组合边界条件进行板振动模态的计算，并考虑板内任意位置点两个方向的转角对板振动模态计算效果的影响，得到以振幅和两个转角为变量的振动方程，计算出来的结果更加准确，具有实用性广、方便应用的优点。

说明书

技术领域

本发明涉及一种基于谱有限元的矩形板振动模态计算方法。

背景技术

板结构在工程应用中占有非常重要的地位，其在航空、土木、电子工程中均有广泛的应用。对于板结构进行振动分析是工程应用设计中的关键，也是振动响应分析的基础和重要组成部分。

针对板结构的自由振动模态计算的问题，众多学者已经开展了卓有成效的研究。Y.Xiang针对一边带有简支边界条件的阶梯矩形Mindlin板采用域分解计算方法对其振动频率进行了计算，得到其分析解。Yufeng Xing利用直接分离变量的方法获得了矩形Mindlin板自由振动闭环形式的解。J.M.Lee基于铁木辛柯梁函数采用Kantorovich方法获得了Mindlin板特征函数，利用迭代法对矩形均质板进行了自由振动分析。Gang Wang对板的自由振动问题和振动模态计算问题进行了一系列的研究，主要有采用Kantorovich–Krylov变分方法针对矩形板的板内振动模态，弯曲振动问题进行了研究，也采用谱有限元法对阶梯板的自由振动问题进行了研究，但是其在计算过程中并没有考虑板中任意一点两个方向剪切角的影响，针对相对较厚的板进行计算时精度有待提高。M.Boscolo采用一阶剪切理论得到板动态刚度元对矩形板的自由振动进行了准确分析，然而此方法只考虑了矩形板某对边是简支边界条件的情况，尚未考虑应用到矩形板具有其他边界条件的情形。

发明内容

为了解决上述技术问题，本发明提供一种准确性高、实用性强的基于谱有限元的矩形板振动模态计算方法。

本发明解决上述问题的技术方案是：一种基于谱有限元的矩形板振动模态计算方法，包括以下步骤：

步骤一：基于能量函数变分原理得到板振动的边界值问题模型，再基于变量分离方法得到两个方向的自由振动方程；

步骤二：确定矩形板两个方向的边界条件，在预设定的模态阶数下采用不同方向相互迭代来计算振动频率和振型，直至两个方向计算得到的振动频率误差在某个指定范围时停止；

步骤三：利用最终得到的两个方向的振动模态进行叠加得到矩形板振动模态。

上述基于谱有限元的矩形板振动模态计算方法，所述步骤一具体步骤为

建立矩形板直角坐标系(x,y,z)，分别定义x,y,z三个方向的位移多元函数u(x,y,t)，v(x,y,t)，w(x,y,t)如下:

u(x,y,t)＝zθ_y(x,y,t)；

v(x,y,t)＝-zθ_x(x,y,t)；

w(x,y,t)＝w₀(x,y,t)；

其中w₀(x,y,t)表示板(x,y)位置在z方向t时刻的振动幅值，θ_x表示板(x,y)位置在t时刻绕x轴的旋转角，θ_y表示板(x,y)位置在t时刻绕y轴的旋转角；

定义能量函数：

Π＝U+T，

其中表示w(x,y,t)的对时间的导数表示θ_x(x,y,t)对时间的导数表示θ_y(x,y,t)对时间的导数常数常数

参数定义为：E：板抗弯刚度，a：板长度，b：板宽度，h：板厚度，v：泊松比，ρ:板材料密度；

对w(x,y,t)，θ_x(x,y,t)和θ_y(x,y,t)进行时空变量分离：

w(x,y,t)＝W(x,y)e^jωt,

其中二元变量W(x,y)表示板(x,y)位置的振动幅值，二元变量表示板(x,y)位置绕x轴的旋转角，二元变量表示板(x,y)位置绕y轴的旋转角；ω:圆频率(rad/s)；

求Π的变分且令变分等于0，即

对变量W(x,y),进行二元变量分离有：

W(x,y)＝W_x(x)W_y(y)

其中W_x(x)，W_y(y)分别表示W(x,y)在x,y方向的分离模态函数，分别表示在x,y方向的分离模态函数，分别表示在x,y方向的分离模态函数

假定已知y向振动模态函数W_y,利用如下变分公式：

δW(x,y)＝W_yδW_x

得到关于x向三个变量的常微分方程组Ⅰ：

其中

假定已知x向振动模态函数W_x,则利用如下变分公式：

δW(x,y)＝W_xδW_y

得到关于y向三个变量的常微分方程组Ⅱ：

其中

上述基于谱有限元的矩形板振动模态计算方法，所述步骤二中，

针对x方向的边界条件定义为三种：

简支边界：W_x＝0,M_xx＝0；

夹紧边界：W_x＝0,

自由边界：V_xz＝0,M_xx＝0,M_xy＝0；

其中力矩经过计算后得到如下的计算用边界条件：

简支边界：W_x＝0,

夹紧边界：W_x＝0,

自由边界：

针对y方向的边界条件定义为三种：

简支边界：W_y＝0,M_yy＝0；

夹紧边界：W_y＝0,

自由边界：V_yz＝0,M_yy＝0,M_xy＝0；

其中力矩

经过计算后得到如下的计算用边界条件：

简支边界：W_y＝0,

夹紧边界：W_y＝0,

自由边界：

上述基于谱有限元的矩形板振动模态计算方法，所述步骤二中，采用不同方向相互迭代计算振动频率和振型函数的具体步骤如下

2-1)设定y向为第m阶模态，根据y向边界条件选择相应自由振动梁的模态函数作为y向的模态函数；

2-2)利用给定的y向模态函数计算得到x向的常微分方程组Ⅰ，根据常微分方程特征值的情况确定x向模态函数的形式，根据边界条件建立以模态函数未知系数为变量的线性方程组，对线性方程组6×6系数矩阵行列式值为0的非线性方程进行优化求解，得到频率值w_x,再利用w_x代入线性方程组进行求解得到模态函数系数，从而得到x向的模态函数；

2-3)利用步骤2-2中计算得到的x向的模态函数，计算得到y向的常微分方程组Ⅱ，根据常微分方程特征值的情况确定y向模态函数的形式，根据边界条件建立以模态函数未知系数为变量的线性方程组，利用系数矩阵的行列式为0进行优化求解得到频率值w_y，再利用w_y代入线性方程组进行求解得到模态函数系数，从而得到y向的模态函数；

2-4)比较步骤2-2和2-3得到的频率值的大小，如果满足|w_x-w_y|≤ε，则退出迭代，ε为误差值，取值为0.0001；如果不满足条件，则将步骤2-3中得到的y向的模态函数作为给定的步骤2-2中模态函数进行迭代计算。

上述基于谱有限元的矩形板振动模态计算方法，所述步骤2-2中，基于已知y向模态函数计算得到频率值w_x和x向的模态函数方法如下：

根据常微分方程组Ⅰ，令表示微分算子d/dx，则有如下的式子：

展开式(3)中的行列式得到如下的方程:

其中ψ＝w_x,or

d₁＝a₃+b₃+c₃-b₁c₁-a₁c₂；

d₂＝a₃b₃+a₃c₃+b₃c₃+a₂b₁c₂+a₁b₂c₁-a₃b₁c₁-a₂b₂-a₁b₃c₂；

d₃＝a₃b₃c₃-a₂b₂c₃；

将试验解ψ＝e^λ代入方程(4)，由于表示微分算子，可产生如下的方程：

λ⁶+d₁λ⁴+d₂λ²+d₃＝0(5)

令μ＝λ²，则

μ³+d₁μ²+d₂μ+d₃＝0(6)

令判别式Δ＝18d₁d₂d₃-4d₁³d₃+d₁²d₂²-4d₂³-27d₃²，

定义w_x(x)与函数表达式中的系数关系数值δ₁,δ₂,δ₃，与函数表达式中的系数关系数值γ₁,γ₂,γ₃如下

其中i＝1,2,3，k为调节系数，为展开式用到的正弦函数频率值，根据展开式不同的项来进行取整数值，m＝1,2,…,∞；r_i为根据式(6)的解计算得到的值；根据判别式的符号以及式(6)中解μ₁,μ₂,μ₃的情形得到如下10种x向振动模态函数形式：

(1)判别式小于0，μ₁>0：实根，μ₂,μ₃：共轭复根

r₂,r₃为共轭复数，Re表示实部，Im表示虚部；

w_x(x)＝-δ₁B₂cosh(r₁x)-δ₁B₁sinh(r₁x)-δ₂B₄cosh(Re(r₂)x)-δ₂B₃sinh(Re(r₂)x)

-δ₃B₆cos(Im(r₂)x)-δ₃B₅sin(Im(r₂)x)

(2)判别式小于0，μ₁<0：实根，μ₂,μ₃：共轭复根

令r₂,r₃为共轭复数，Re表示实部，Im表示虚部；

w_x(x)＝-δ₁B₂cos(r₁x)-δ₁B₁sin(r₁x)-δ₂B₄cosh(Re(r₂)x)-δ₂B₃sinh(Re(r₂)x)

-δ₃B₆cos(Im(r₂)x)-δ₃B₅sin(Im(r₂)x)

(3)判别式大于0，μ₁>0：实根，μ₂>0：实根，μ₃>0：实根

令

(4)判别式大于0，μ₁>0：实根，μ₂>0：实根，μ₃<0：实根

(5)判别式大于0，μ₁>0：实根，μ₂<0：实根，μ₃>0：实根

(6)判别式大于0，μ₁<0：实根，μ₂>0：实根，μ₃>0：实根

(7)判别式大于0，μ₁>0：实根，μ₂<0：实根，μ₃<0：实根

(8)判别式大于0，μ₁<0：实根，μ₂>0：实根，μ₃<0：实根

(9)判别式大于0，μ₁<0：实根，μ₂<0：实根，μ₃>0：实根

(10)判别式大于0，μ₁<0：实根，μ₂<0：实根，μ₃<0：实根

将得到的模态函数表达式代入边界条件，整理以后得到如下形式的方程组：

其中表示线性方程组系数矩阵，A₁₁～A₆₆表示对应位置的系数值；令采用matlab中的非线性优化函数fsolve进行求解得到ω_x，再将ω_x代入式(7)中得到线性方程组，对线性方程组求解得到B₁,B₂,…,B₆的值，从而得到x向的振动模态函数的表达式。

上述基于谱有限元的矩形板振动模态计算方法，所述步骤2-3中基于已知x向模态函数计算得到频率值w_y和y向模态函数方法如下：

根据常微分方程组Ⅱ，令表示微分算子d/dy，则有如下的式子：

展开式(8)中的行列式得到如下的方程:

其中ψ＝w_y,or

j₁＝e₃+f₃+g₃-e₁f₂-g₁f₁；

j₂＝g₃f₃+e₃g₃+e₃f₃+e₁f₁g₂+e₂f₂g₁-e₃g₁f₁-e₂g₂

j₃＝e₃g₃f₃-e₂g₂f₃-e₁g₃f₂；

将试验解ψ＝e^λ代入方程(9)，由于表示微分算子，则产生如下的方程：

λ⁶+j₁λ⁴+j₂λ²+j₃＝0(10)

令μ＝λ²，则

μ³+j₁μ²+j₂μ+j₃＝0(11)

令判别式Δ＝18j₁j₂j₃-4j₁³j₃+j₁²j₂²-4j₂³-27j₃²，

定义w_y(y)与函数表达式中的系数关系数值仍表示为δ₁,δ₂,δ₃，与函数表达式中的系数关系数值仍表示为γ₁,γ₂,γ₃如下

其中i＝1,2,3，k为调节系数，为展开式用到的正弦函数频率值，根据展开式不同的项来进行取整数值，m＝1,2,…,∞；r_i为根据式(11)的解计算得到的值；根据判别式的符号以及式(11)中解μ₁,μ₂,μ₃的情形得到如下10种y向振动模态函数形式：

(11)判别式小于0，μ₁>0：实根，μ₂,μ₃：共轭复根

r₂,r₃为共轭复数，Re表示实部，Im表示虚部；

w_y(y)＝-δ₁B₂cosh(r₁y)-δ₁B₁sinh(r₁y)-δ₂B₄cosh(Re(r₂)y)-δ₂B₃sinh(Re(r₂)y)

-δ₃B₆cos(Im(r₂)y)-δ₃B₅sin(Im(r₂)y)

(12)判别式小于0，μ₁<0：实根，μ₂,μ₃：共轭复根

令r₂,r₃为共轭复数，Re表示实部，Im表示虚部；

w_y(y)＝-δ₁B₂cos(r₁y)-δ₁B₁sin(r₁y)-δ₂B₄cosh(Re(r₂)y)-δ₂B₃sinh(Re(r₂)y)

-δ₃B₆cos(Im(r₂)y)-δ₃B₅sin(Im(r₂)y)

(13)判别式大于0，μ₁>0：实根，μ₂>0：实根，μ₃>0：实根

令

(14)判别式大于0，μ₁>0：实根，μ₂>0：实根，μ₃<0：实根

(15)判别式大于0，μ₁>0：实根，μ₂<0：实根，μ₃>0：实根

(16)判别式大于0，μ₁<0：实根，μ₂>0：实根，μ₃>0：实根

(17)判别式大于0，μ₁>0：实根，μ₂<0：实根，μ₃<0：实根

(18)判别式大于0，μ₁<0：实根，μ₂>0：实根，μ₃<0：实根

(19)判别式大于0，μ₁<0：实根，μ₂<0：实根，μ₃>0：实根

(20)判别式大于0，μ₁<0：实根，μ₂<0：实根，μ₃<0：实根

将得到的表达式代入边界条件，整理以后得到如下形式的方程组：

令系数矩阵的行列式采用matlab中的非线性优化函数fsolve进行求解得到ω_y，再将ω_y代入方程(12)中得到线性方程组，对线性方程组求解得到B₁,B₂,…,B₆的值，从而得到y向的振动模态函数的表达式。

本发明的有益效果在于：本发明可针对多种组合边界条件进行板振动模态的计算，并考虑板内任意位置点两个方向的旋转角对板振动模态计算的影响，得到以振幅和两个旋转角为变量的振动方程，更加适应计算对象板厚度变化时等实际情况，计算出来的结果更加准确，具有实用性广、方便应用的优点。

附图说明

图1为本发明的流程图。

图2为矩形板坐标系图。

图3为预选定的y向初始模态函数及其一、二阶导数图。

图4为x向振幅一阶模态图。

图5为x向绕y轴剪切角一阶模态图。

图6为x向绕x轴剪切角一阶模态图。

图7为y向振幅一阶模态图。

图8为y向绕x轴剪切角一阶模态图。

图9为y向绕y轴剪切角一阶模态图。

图10为矩形板振动幅值一阶空间模态图。

图11为矩形板绕x轴剪切角一阶空间模态图。

图12为矩形板绕y轴剪切角一阶空间模态图。

具体实施方式

下面结合附图和实施例对本发明作进一步的说明。

一种基于谱有限元的矩形板振动模态计算方法，包括以下步骤：

步骤一：基于能量函数变分原理得到板振动的边界值问题模型，再基于变量分离方法得到两个方向的自由振动方程。

具体步骤为

建立矩形板直角坐标系(x,y,z)，分别定义x,y,z三个方向的位移多元函数u(x,y,t)，v(x,y,t)，w(x,y,t)如下:

u(x,y,t)＝zθ_y(x,y,t)；

v(x,y,t)＝-zθ_x(x,y,t)；

w(x,y,t)＝w₀(x,y,t)；

其中w₀(x,y,t)表示板(x,y)位置在z方向t时刻的振动幅值，θ_x表示板(x,y)位置在t时刻绕x轴的旋转角，θ_y表示板(x,y)位置在t时刻绕y轴的旋转角；

定义能量函数：

Π＝U+T，

其中表示w(x,y,t)的对时间的导数表示θ_x(x,y,t)对时间的导数表示θ_y(x,y,t)对时间的导数常数常数

参数定义为：E：板抗弯刚度，a：板长度，b：板宽度，h：板厚度，v：泊松比，ρ:板材料密度；

对w(x,y,t)，θ_x(x,y,t)和θ_y(x,y,t)进行时空变量分离：

w(x,y,t)＝W(x,y)e^jωt,

其中二元变量W(x,y)表示板(x,y)位置的振动幅值，二元变量表示板(x,y)位置绕x轴的旋转角，二元变量表示板(x,y)位置绕y轴的旋转角。ω:圆频率(rad/s)；

求Π的变分且令变分等于0，即

对变量W(x,y),进行变量分离有：

W(x,y)＝W_x(x)W_y(y)

其中W_x(x)，W_y(y)分别表示W(x,y)在x,y方向的分离模态函数，分别表示在x,y方向的分离模态函数，分别表示在x,y方向的分离模态函数。

假定已知y向振动模态函数W_y,利用如下变分公式：

δW(x,y)＝W_yδW_x

得到关于x向三个变量的常微分方程组Ⅰ：

其中

假定已知x向振动模态函数W_x,则利用如下变分公式：

δW(x,y)＝W_xδW_y

得到关于y向三个变量的常微分方程组Ⅱ：

其中

步骤二：确定矩形板两个方向的边界条件，在预设定的模态阶数下采用不同方向相互迭代来计算振动频率和振型，直至两个方向计算得到的振动频率误差在某个指定范围时停止。

针对x方向的边界条件定义为三种：

简支边界：W_x＝0,M_xx＝0；

夹紧边界：W_x＝0,

自由边界：V_xz＝0,M_xx＝0,M_xy＝0；

其中力矩

经过计算后得到如下的计算用边界条件：

简支边界：W_x＝0,

夹紧边界：W_x＝0,

自由边界：

针对y方向的边界条件定义为三种：

简支边界：W_y＝0,M_yy＝0；

夹紧边界：W_y＝0,

自由边界：V_yz＝0,M_yy＝0,M_xy＝0；

其中力矩

经过计算后得到如下的计算用边界条件：

简支边界：W_y＝0,

夹紧边界：W_y＝0,

自由边界：

采用不同方向相互迭代计算振动频率和振型函数的具体步骤如下

2-1)设定y(或者x)向为第m阶模态，根据y(或者x)向边界条件选择相应自由振动梁的模态函数作为y(或者x)向的模态函数；

2-2)利用给定的y(或者x)向模态函数计算得到x(或者y)向的常微分方程组Ⅰ(或者Ⅱ)，根据常微分方程特征值的情况确定x(或者y)向模态函数的形式，根据边界条件建立以模态函数未知系数为变量的线性方程组，对线性方程组6×6系数矩阵行列式值为0的非线性方程进行优化求解，得到频率值w_x,再利用w_x代入线性方程组进行求解得到模态函数系数，从而得到x(或者y)向的模态函数；

2-3)利用步骤2-2中计算得到的x(或者y)向的模态函数，计算得到y(或者x)向的常微分方程组Ⅱ(或者Ⅰ)，根据常微分方程特征值的情况确定y(或者x)向模态函数的形式，根据边界条件建立以模态函数未知系数为变量的线性方程组，利用系数矩阵的行列式为0进行优化求解得到频率值w_y，再利用w_y代入线性方程组进行求解得到模态函数系数，从而得到y(或者x)向的模态函数；

2-4)比较步骤2-2和2-3得到的频率值的大小，如果满足|w_x-w_y|≤ε，则退出迭代，ε为误差值，取值为0.0001；如果不满足条件，则将步骤2-3中得到的y向或者x向的模态函数作为给定的步骤2-2中模态函数进行迭代计算。

所述步骤2-2中，基于已知y向模态函数计算得到频率值w_x和x向的模态函数方法如下：

根据常微分方程组Ⅰ，令表示微分算子d/dx，则有如下的式子：

展开式(3)中的行列式得到如下的方程:

其中ψ＝w_x,or

d₁＝a₃+b₃+c₃-b₁c₁-a₁c₂；

d₂＝a₃b₃+a₃c₃+b₃c₃+a₂b₁c₂+a₁b₂c₁-a₃b₁c₁-a₂b₂-a₁b₃c₂；

d₃＝a₃b₃c₃-a₂b₂c₃；

将试验解ψ＝e^λ代入方程(4)，由于D表示微分算子，可产生如下的方程：

λ⁶+d₁λ⁴+d₂λ²+d₃＝0(5)

令μ＝λ²，则

μ³+d₁μ²+d₂μ+d₃＝0(6)

令判别式Δ＝18d₁d₂d₃-4d₁³d₃+d₁²d₂²-4d₂³-27d₃²，

定义w_x(x)与函数表达式中的系数关系数值δ₁,δ₂,δ₃，与函数表达式中的系数关系数值γ₁,γ₂,γ₃如下

其中i＝1,2,3，k为调节系数，为展开式用到的正弦函数频率值，根据展开式不同的项来进行取整数值，m＝1,2,…,∞；r_i为根据式(6)的解计算得到的值；根据判别式的符号以及式(6)中解μ₁,μ₂,μ₃的情形得到如下10种x向振动模态函数形式：

(1)判别式小于0，μ₁>0：实根，μ₂,μ₃：共轭复根

r₂,r₃为共轭复数，Re表示实部，Im表示虚部；

w_x(x)＝-δ₁B₂cosh(r₁x)-δ₁B₁sinh(r₁x)-δ₂B₄cosh(Re(r₂)x)-δ₂B₃sinh(Re(r₂)x)

-δ₃B₆cos(Im(r₂)x)-δ₃B₅sin(Im(r₂)x)

(2)判别式小于0，μ₁<0：实根，μ₂,μ₃：共轭复根

令r₂,r₃为共轭复数，Re表示实部，Im表示虚部；

w_x(x)＝-δ₁B₂cos(r₁x)-δ₁B₁sin(r₁x)-δ₂B₄cosh(Re(r₂)x)-δ₂B₃sinh(Re(r₂)x)

-δ₃B₆cos(Im(r₂)x)-δ₃B₅sin(Im(r₂)x)

(3)判别式大于0，μ₁>0：实根，μ₂>0：实根，μ₃>0：实根

令

(4)判别式大于0，μ₁>0：实根，μ₂>0：实根，μ₃<0：实根

(5)判别式大于0，μ₁>0：实根，μ₂<0：实根，μ₃>0：实根

(6)判别式大于0，μ₁<0：实根，μ₂>0：实根，μ₃>0：实根

(7)判别式大于0，μ₁>0：实根，μ₂<0：实根，μ₃<0：实根

(8)判别式大于0，μ₁<0：实根，μ₂>0：实根，μ₃<0：实根

(9)判别式大于0，μ₁<0：实根，μ₂<0：实根，μ₃>0：实根

(10)判别式大于0，μ₁<0：实根，μ₂<0：实根，μ₃<0：实根

将得到的模态函数表达式代入边界条件，整理以后得到如下形式的方程组：

其中表示线性方程组系数矩阵，A₁₁～A₆₆表示对应位置的系数值；令采用matlab中的非线性优化函数fsolve进行求解得到ω_x，再将ω_x代入式(7)中得到线性方程组，对线性方程组求解得到B₁,B₂,…,B₆的值，从而得到x向的振动模态函数的表达式。

所述步骤2-3中基于已知x向模态函数计算得到频率值w_y和y向模态函数方法如下：

根据常微分方程组Ⅱ，令表示微分算子d/dy，则有如下的式子：

展开式(8)中的行列式得到如下的方程:

其中ψ＝w_y,or

j₁＝e₃+f₃+g₃-e₁f₂-g₁f₁；

j₂＝g₃f₃+e₃g₃+e₃f₃+e₁f₁g₂+e₂f₂g₁-e₃g₁f₁-e₂g₂

j₃＝e₃g₃f₃-e₂g₂f₃-e₁g₃f₂；

将试验解ψ＝e^λ代入方程(9)，由于表示微分算子，则产生如下的方程：

λ⁶+j₁λ⁴+j₂λ²+j₃＝0(10)

令μ＝λ²，则

μ³+j₁μ²+j₂μ+j₃＝0(11)

令判别式Δ＝18j₁j₂j₃-4j₁³j₃+j₁²j₂²-4j₂³-27j₃²，

定义w_y(y)与函数表达式中的系数关系数值仍表示为δ₁,δ₂,δ₃，与函数表达式中的系数关系数值仍表示为γ₁,γ₂,γ₃如下

其中i＝1,2,3，k为调节系数，为展开式用到的正弦函数频率值，根据展开式不同的项来进行取整数值，m＝1,2,…,∞；r_i为根据式(11)的解计算得到的值；根据判别式的符号以及式(11)中解μ₁,μ₂,μ₃的情形得到如下10种y向振动模态函数形式：

(11)判别式小于0，μ₁>0：实根，μ₂,μ₃：共轭复根

r₂,r₃为共轭复数，Re表示实部，Im表示虚部；

w_y(y)＝-δ₁B₂cosh(r₁y)-δ₁B₁sinh(r₁y)-δ₂B₄cosh(Re(r₂)y)-δ₂B₃sinh(Re(r₂)y)

-δ₃B₆cos(Im(r₂)y)-δ₃B₅sin(Im(r₂)y)

(12)判别式小于0，μ₁<0：实根，μ₂,μ₃：共轭复根

令r₂,r₃为共轭复数，Re表示实部，Im表示虚部；

w_y(y)＝-δ₁B₂cos(r₁y)-δ₁B₁sin(r₁y)-δ₂B₄cosh(Re(r₂)y)-δ₂B₃sinh(Re(r₂)y)

-δ₃B₆cos(Im(r₂)y)-δ₃B₅sin(Im(r₂)y)

(13)判别式大于0，μ₁>0：实根，μ₂>0：实根，μ₃>0：实根

令

(14)判别式大于0，μ₁>0：实根，μ₂>0：实根，μ₃<0：实根

(15)判别式大于0，μ₁>0：实根，μ₂<0：实根，μ₃>0：实根

(16)判别式大于0，μ₁<0：实根，μ₂>0：实根，μ₃>0：实根

(17)判别式大于0，μ₁>0：实根，μ₂<0：实根，μ₃<0：实根

(18)判别式大于0，μ₁<0：实根，μ₂>0：实根，μ₃<0：实根

(19)判别式大于0，μ₁<0：实根，μ₂<0：实根，μ₃>0：实根

(20)判别式大于0，μ₁<0：实根，μ₂<0：实根，μ₃<0：实根

将得到的表达式代入边界条件，整理以后得到如下形式的方程组：

令采用matlab中的非线性优化函数fsolve进行求解得到ω_y再将ω_y代入方程(12)中得到线性方程组，对线性方程组求解得到B₁,B₂,…,B₆的值，从而得到y向的振动模态函数的表达式。

步骤三：利用最终得到的两个方向的振动模态进行叠加得到矩形板振动模态。

具体实例选择铝合金矩形板作为计算对象，具体的参数如下：

板长a＝1m，板宽b＝1m，板厚h＝0.1m,板材料密度ρ＝2700kg/m³，杨氏模量E＝70×10⁹N/m²,泊松比v＝0.3，调节系数k＝5/6，D＝6.4103×10⁶,G＝2.6923×10¹⁰。

边界条件设定为y向均为简支边界条件，x向也为简支边界条件。已知y向第一阶模态函数及其导数，需要求解x向一阶模态函数。

本发明的矩形板振动模态计算方法，具体包括如下步骤：

步骤一：建立如图2坐标系，采用变分法和变量分离方法计算得到如下的两个方程组：

步骤二：确定矩形板两个方向的边界条件，在预设定的模态阶数下采用不同方向相互迭代来计算振动频率和振型，直至两个方向计算得到的振动频率误差在某个指定范围时停止。具体步骤如下：

(Ⅰ)设定y(或者x)向为第m阶模态，根据y(或者x)向边界条件选择相应自由振动梁的模态函数作为y(或者x)向的模态函数；

如图所示，图3中模态函数为先确定矩形板y向的边界条件为简支时得到的板带y向振动第一阶模态函数及其前两阶导数。

(Ⅱ)利用给定的y(或者x)向模态函数计算得到x(或者y)向的常微分方程组(Ⅰ)(或者(Ⅱ))，根据常微分方程特征值的情况确定x(或者y)向模态函数的形式，根据边界条件建立以模态函数未知系数为变量的线性方程组，对线性方程组6×6系数矩阵行列式值为0的非线性方程进行优化求解，得到频率值w_x,再利用w_x代入线性方程组进行求解得到模态函数系数，从而得到x(或者y)向的模态函数。

则在优化时输入初始值20hz时，得到如下的公式计算。

(D⁶+d₁D⁴+d₂D²+d₃)ψ＝0

其中ψ＝w_x,or

d₁＝a₃+b₃+c₃-b₁c₁-a₁c₂；

d₂＝a₃b₃+a₃c₃+b₃c₃+a₂b₁c₂+a₁b₂c₁-a₃b₁c₁-a₂b₂-a₁b₃c₂；

d₃＝a₃b₃c₃-a₂b₂c₃。

其中d₁＝-1027.4，d₂＝18284.5，d₃＝281500；

用e^λ代入(D⁶+d₁D⁴+d₂D²+d₃)ψ＝0产生如下的方程：

λ⁶+d₁λ⁴+d₂λ²+d₃＝0

令μ＝λ²，则

μ³+d₁μ²+d₂μ+d₃＝0

则判别式Δ＝18d₁d₂d₃-4d₁³d₃+d₁²d₂²-4d₂³-27d₃²>0，

此时特征根分别为：

μ₁＝1009,μ₂＝28.2672，μ₃＝-9.8696

由于判别式大于0，μ₁>0：实根，μ₂>0：实根，μ₃<0：实根

其中模态函数表达式为：

其中δ₁＝0,δ₂＝0.1780,δ₃＝0.3174，γ₁＝10.1111,γ₂＝0.0989,γ₃＝0.8951。

将得到的表达式代入边界条件，整理以后得到如下形式的方程组：

令

采用matlab中的非线性优化函数进行求解得到ω_x。

再将ω_x代入式(3)到线性方程组，对线性方程组求解得到B₁,B₂,…,B₆的值，从而得到x向的振动模态函数的表达式。具体的归一化后的x向幅值、绕x轴转角、绕y轴转角的第一阶模态如图4-6所示。

(Ⅲ)利用步骤(Ⅱ)中计算得到的x(或者y)的模态函数，计算得到y(或者x)向的常微分方程组(Ⅱ)(或者(Ⅰ))，根据常微分方程特征值的情况确定y(或者x)向模态函数的形式，根据边界条件建立以模态函数未知系数为变量的线性方程组，利用系数矩阵的行列式为0进行优化求解得到频率值w_y，再利用w_y代入线性方程组进行求解得到模态函数系数，从而得到y(或者x)向的模态函数。

按照如(Ⅱ)中相同步骤进行计算得到y向一阶模态函数如图7-9所示。

(Ⅳ)比较步骤(Ⅱ)和(Ⅲ)得到的频率值的大小，如果满足|w_x-w_y|≤ε(ε为误差值，取0.0001)，则退出迭代。如果不满足条件，则将步骤(Ⅲ)中得到的y(或者x)向的模态函数作为给定的步骤(Ⅱ)模态函数进行迭代计算。

步骤三：利用最终得到的两个方向的一阶振动模态进行叠加得到矩形板振动模态，如图10-12所示。