CN106777825A - 一种基于谱有限元的矩形板振动模态计算方法 - Google Patents
一种基于谱有限元的矩形板振动模态计算方法 Download PDFInfo
- Publication number
- CN106777825A CN106777825A CN201710059235.1A CN201710059235A CN106777825A CN 106777825 A CN106777825 A CN 106777825A CN 201710059235 A CN201710059235 A CN 201710059235A CN 106777825 A CN106777825 A CN 106777825A
- Authority
- CN
- China
- Prior art keywords
- theta
- overbar
- gamma
- delta
- sinh
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于谱有限元的矩形板振动模态计算方法,包括以下步骤:基于能量函数变分原理得到板振动的边界值问题模型,再基于变量分离方法得到两个方向的自由振动方程;确定矩形板两个方向的边界条件,在预设定的模态阶数下采用不同方向相互迭代来计算振动频率和振型,直至两个方向计算得到的振动频率误差在某个指定范围时停止;利用最终得到的两个方向的振动模态进行叠加得到矩形板振动模态。本发明针对多种组合边界条件进行板振动模态的计算,并考虑板内任意位置点两个方向的转角对板振动模态计算效果的影响,得到以振幅和两个转角为变量的振动方程,计算出来的结果更加准确,具有实用性广、方便应用的优点。
Description
技术领域
本发明涉及一种基于谱有限元的矩形板振动模态计算方法。
背景技术
板结构在工程应用中占有非常重要的地位,其在航空、土木、电子工程中均有广泛的应用。对于板结构进行振动分析是工程应用设计中的关键,也是振动响应分析的基础和重要组成部分。
针对板结构的自由振动模态计算的问题,众多学者已经开展了卓有成效的研究。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)=w0(x,y,t);
其中w0(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)ejωt,
其中二元变量W(x,y)表示板(x,y)位置的振动幅值,二元变量表示板(x,y)位置绕x轴的旋转角,二元变量表示板(x,y)位置绕y轴的旋转角;ω:圆频率(rad/s);
求Π的变分且令变分等于0,即
对变量W(x,y),进行二元变量分离有:
W(x,y)=Wx(x)Wy(y)
其中Wx(x),Wy(y)分别表示W(x,y)在x,y方向的分离模态函数,分别表示在x,y方向的分离模态函数,分别表示在x,y方向的分离模态函数
假定已知y向振动模态函数Wy,利用如下变分公式:
δW(x,y)=WyδWx
得到关于x向三个变量的常微分方程组Ⅰ:
其中
假定已知x向振动模态函数Wx,则利用如下变分公式:
δW(x,y)=WxδWy
得到关于y向三个变量的常微分方程组Ⅱ:
其中
上述基于谱有限元的矩形板振动模态计算方法,所述步骤二中,
针对x方向的边界条件定义为三种:
简支边界:Wx=0,Mxx=0;
夹紧边界:Wx=0,
自由边界:Vxz=0,Mxx=0,Mxy=0;
其中力矩经过计算后得到如下的计算用边界条件:
简支边界:Wx=0,
夹紧边界:Wx=0,
自由边界:
针对y方向的边界条件定义为三种:
简支边界:Wy=0,Myy=0;
夹紧边界:Wy=0,
自由边界:Vyz=0,Myy=0,Mxy=0;
其中力矩
经过计算后得到如下的计算用边界条件:
简支边界:Wy=0,
夹紧边界:Wy=0,
自由边界:
上述基于谱有限元的矩形板振动模态计算方法,所述步骤二中,采用不同方向相互迭代计算振动频率和振型函数的具体步骤如下
2-1)设定y向为第m阶模态,根据y向边界条件选择相应自由振动梁的模态函数作为y向的模态函数;
2-2)利用给定的y向模态函数计算得到x向的常微分方程组Ⅰ,根据常微分方程特征值的情况确定x向模态函数的形式,根据边界条件建立以模态函数未知系数为变量的线性方程组,对线性方程组6×6系数矩阵行列式值为0的非线性方程进行优化求解,得到频率值wx,再利用wx代入线性方程组进行求解得到模态函数系数,从而得到x向的模态函数;
2-3)利用步骤2-2中计算得到的x向的模态函数,计算得到y向的常微分方程组Ⅱ,根据常微分方程特征值的情况确定y向模态函数的形式,根据边界条件建立以模态函数未知系数为变量的线性方程组,利用系数矩阵的行列式为0进行优化求解得到频率值wy,再利用wy代入线性方程组进行求解得到模态函数系数,从而得到y向的模态函数;
2-4)比较步骤2-2和2-3得到的频率值的大小,如果满足|wx-wy|≤ε,则退出迭代,ε为误差值,取值为0.0001;如果不满足条件,则将步骤2-3中得到的y向的模态函数作为给定的步骤2-2中模态函数进行迭代计算。
上述基于谱有限元的矩形板振动模态计算方法,所述步骤2-2中,基于已知y向模态函数计算得到频率值wx和x向的模态函数方法如下:
根据常微分方程组Ⅰ,令表示微分算子d/dx,则有如下的式子:
展开式(3)中的行列式得到如下的方程:
其中ψ=wx,or
d1=a3+b3+c3-b1c1-a1c2;
d2=a3b3+a3c3+b3c3+a2b1c2+a1b2c1-a3b1c1-a2b2-a1b3c2;
d3=a3b3c3-a2b2c3;
将试验解ψ=eλ代入方程(4),由于表示微分算子,可产生如下的方程:
λ6+d1λ4+d2λ2+d3=0 (5)
令μ=λ2,则
μ3+d1μ2+d2μ+d3=0 (6)
令判别式Δ=18d1d2d3-4d1 3d3+d1 2d2 2-4d2 3-27d3 2,
定义wx(x)与函数表达式中的系数关系数值δ1,δ2,δ3,与函数表达式中的系数关系数值γ1,γ2,γ3如下
其中i=1,2,3,k为调节系数,为展开式用到的正弦函数频率值,根据展开式不同的项来进行取整数值,m=1,2,…,∞;ri为根据式(6)的解计算得到的值;根据判别式的符号以及式(6)中解μ1,μ2,μ3的情形得到如下10种x向振动模态函数形式:
(1)判别式小于0,μ1>0:实根,μ2,μ3:共轭复根
r2,r3为共轭复数,Re表示实部,Im表示虚部;
wx(x)=-δ1B2cosh(r1x)-δ1B1sinh(r1x)-δ2B4cosh(Re(r2)x)-δ2B3sinh(Re(r2)x)
-δ3B6cos(Im(r2)x)-δ3B5sin(Im(r2)x)
(2)判别式小于0,μ1<0:实根,μ2,μ3:共轭复根
令r2,r3为共轭复数,Re表示实部,Im表示虚部;
wx(x)=-δ1B2cos(r1x)-δ1B1sin(r1x)-δ2B4cosh(Re(r2)x)-δ2B3sinh(Re(r2)x)
-δ3B6cos(Im(r2)x)-δ3B5sin(Im(r2)x)
(3)判别式大于0,μ1>0:实根,μ2>0:实根,μ3>0:实根
令
(4)判别式大于0,μ1>0:实根,μ2>0:实根,μ3<0:实根
(5)判别式大于0,μ1>0:实根,μ2<0:实根,μ3>0:实根
(6)判别式大于0,μ1<0:实根,μ2>0:实根,μ3>0:实根
(7)判别式大于0,μ1>0:实根,μ2<0:实根,μ3<0:实根
(8)判别式大于0,μ1<0:实根,μ2>0:实根,μ3<0:实根
(9)判别式大于0,μ1<0:实根,μ2<0:实根,μ3>0:实根
(10)判别式大于0,μ1<0:实根,μ2<0:实根,μ3<0:实根
将得到的模态函数表达式代入边界条件,整理以后得到如下形式的方程组:
其中表示线性方程组系数矩阵,A11~A66表示对应位置的系数值;令采用matlab中的非线性优化函数fsolve进行求解得到ωx,再将ωx代入式(7)中得到线性方程组,对线性方程组求解得到B1,B2,…,B6的值,从而得到x向的振动模态函数的表达式。
上述基于谱有限元的矩形板振动模态计算方法,所述步骤2-3中基于已知x向模态函数计算得到频率值wy和y向模态函数方法如下:
根据常微分方程组Ⅱ,令表示微分算子d/dy,则有如下的式子:
展开式(8)中的行列式得到如下的方程:
其中ψ=wy,or
j1=e3+f3+g3-e1f2-g1f1;
j2=g3f3+e3g3+e3f3+e1f1g2+e2f2g1-e3g1f1-e2g2
j3=e3g3f3-e2g2f3-e1g3f2;
将试验解ψ=eλ代入方程(9),由于表示微分算子,则产生如下的方程:
λ6+j1λ4+j2λ2+j3=0 (10)
令μ=λ2,则
μ3+j1μ2+j2μ+j3=0 (11)
令判别式Δ=18j1j2j3-4j1 3j3+j1 2j2 2-4j2 3-27j3 2,
定义wy(y)与函数表达式中的系数关系数值仍表示为δ1,δ2,δ3,与函数表达式中的系数关系数值仍表示为γ1,γ2,γ3如下
其中i=1,2,3,k为调节系数,为展开式用到的正弦函数频率值,根据展开式不同的项来进行取整数值,m=1,2,…,∞;ri为根据式(11)的解计算得到的值;根据判别式的符号以及式(11)中解μ1,μ2,μ3的情形得到如下10种y向振动模态函数形式:
(11)判别式小于0,μ1>0:实根,μ2,μ3:共轭复根
r2,r3为共轭复数,Re表示实部,Im表示虚部;
wy(y)=-δ1B2cosh(r1y)-δ1B1sinh(r1y)-δ2B4cosh(Re(r2)y)-δ2B3sinh(Re(r2)y)
-δ3B6cos(Im(r2)y)-δ3B5sin(Im(r2)y)
(12)判别式小于0,μ1<0:实根,μ2,μ3:共轭复根
令r2,r3为共轭复数,Re表示实部,Im表示虚部;
wy(y)=-δ1B2cos(r1y)-δ1B1sin(r1y)-δ2B4cosh(Re(r2)y)-δ2B3sinh(Re(r2)y)
-δ3B6cos(Im(r2)y)-δ3B5sin(Im(r2)y)
(13)判别式大于0,μ1>0:实根,μ2>0:实根,μ3>0:实根
令
(14)判别式大于0,μ1>0:实根,μ2>0:实根,μ3<0:实根
(15)判别式大于0,μ1>0:实根,μ2<0:实根,μ3>0:实根
(16)判别式大于0,μ1<0:实根,μ2>0:实根,μ3>0:实根
(17)判别式大于0,μ1>0:实根,μ2<0:实根,μ3<0:实根
(18)判别式大于0,μ1<0:实根,μ2>0:实根,μ3<0:实根
(19)判别式大于0,μ1<0:实根,μ2<0:实根,μ3>0:实根
(20)判别式大于0,μ1<0:实根,μ2<0:实根,μ3<0:实根
将得到的表达式代入边界条件,整理以后得到如下形式的方程组:
令系数矩阵的行列式采用matlab中的非线性优化函数fsolve进行求解得到ωy,再将ωy代入方程(12)中得到线性方程组,对线性方程组求解得到B1,B2,…,B6的值,从而得到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)=w0(x,y,t);
其中w0(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)ejωt,
其中二元变量W(x,y)表示板(x,y)位置的振动幅值,二元变量表示板(x,y)位置绕x轴的旋转角,二元变量表示板(x,y)位置绕y轴的旋转角。ω:圆频率(rad/s);
求Π的变分且令变分等于0,即
对变量W(x,y),进行变量分离有:
W(x,y)=Wx(x)Wy(y)
其中Wx(x),Wy(y)分别表示W(x,y)在x,y方向的分离模态函数,分别表示在x,y方向的分离模态函数,分别表示在x,y方向的分离模态函数。
假定已知y向振动模态函数Wy,利用如下变分公式:
δW(x,y)=WyδWx
得到关于x向三个变量的常微分方程组Ⅰ:
其中
假定已知x向振动模态函数Wx,则利用如下变分公式:
δW(x,y)=WxδWy
得到关于y向三个变量的常微分方程组Ⅱ:
其中
步骤二:确定矩形板两个方向的边界条件,在预设定的模态阶数下采用不同方向相互迭代来计算振动频率和振型,直至两个方向计算得到的振动频率误差在某个指定范围时停止。
针对x方向的边界条件定义为三种:
简支边界:Wx=0,Mxx=0;
夹紧边界:Wx=0,
自由边界:Vxz=0,Mxx=0,Mxy=0;
其中力矩
经过计算后得到如下的计算用边界条件:
简支边界:Wx=0,
夹紧边界:Wx=0,
自由边界:
针对y方向的边界条件定义为三种:
简支边界:Wy=0,Myy=0;
夹紧边界:Wy=0,
自由边界:Vyz=0,Myy=0,Mxy=0;
其中力矩
经过计算后得到如下的计算用边界条件:
简支边界:Wy=0,
夹紧边界:Wy=0,
自由边界:
采用不同方向相互迭代计算振动频率和振型函数的具体步骤如下
2-1)设定y(或者x)向为第m阶模态,根据y(或者x)向边界条件选择相应自由振动梁的模态函数作为y(或者x)向的模态函数;
2-2)利用给定的y(或者x)向模态函数计算得到x(或者y)向的常微分方程组Ⅰ(或者Ⅱ),根据常微分方程特征值的情况确定x(或者y)向模态函数的形式,根据边界条件建立以模态函数未知系数为变量的线性方程组,对线性方程组6×6系数矩阵行列式值为0的非线性方程进行优化求解,得到频率值wx,再利用wx代入线性方程组进行求解得到模态函数系数,从而得到x(或者y)向的模态函数;
2-3)利用步骤2-2中计算得到的x(或者y)向的模态函数,计算得到y(或者x)向的常微分方程组Ⅱ(或者Ⅰ),根据常微分方程特征值的情况确定y(或者x)向模态函数的形式,根据边界条件建立以模态函数未知系数为变量的线性方程组,利用系数矩阵的行列式为0进行优化求解得到频率值wy,再利用wy代入线性方程组进行求解得到模态函数系数,从而得到y(或者x)向的模态函数;
2-4)比较步骤2-2和2-3得到的频率值的大小,如果满足|wx-wy|≤ε,则退出迭代,ε为误差值,取值为0.0001;如果不满足条件,则将步骤2-3中得到的y向或者x向的模态函数作为给定的步骤2-2中模态函数进行迭代计算。
所述步骤2-2中,基于已知y向模态函数计算得到频率值wx和x向的模态函数方法如下:
根据常微分方程组Ⅰ,令表示微分算子d/dx,则有如下的式子:
展开式(3)中的行列式得到如下的方程:
其中ψ=wx,or
d1=a3+b3+c3-b1c1-a1c2;
d2=a3b3+a3c3+b3c3+a2b1c2+a1b2c1-a3b1c1-a2b2-a1b3c2;
d3=a3b3c3-a2b2c3;
将试验解ψ=eλ代入方程(4),由于D表示微分算子,可产生如下的方程:
λ6+d1λ4+d2λ2+d3=0 (5)
令μ=λ2,则
μ3+d1μ2+d2μ+d3=0 (6)
令判别式Δ=18d1d2d3-4d1 3d3+d1 2d2 2-4d2 3-27d3 2,
定义wx(x)与函数表达式中的系数关系数值δ1,δ2,δ3,与函数表达式中的系数关系数值γ1,γ2,γ3如下
其中i=1,2,3,k为调节系数,为展开式用到的正弦函数频率值,根据展开式不同的项来进行取整数值,m=1,2,…,∞;ri为根据式(6)的解计算得到的值;根据判别式的符号以及式(6)中解μ1,μ2,μ3的情形得到如下10种x向振动模态函数形式:
(1)判别式小于0,μ1>0:实根,μ2,μ3:共轭复根
r2,r3为共轭复数,Re表示实部,Im表示虚部;
wx(x)=-δ1B2cosh(r1x)-δ1B1sinh(r1x)-δ2B4cosh(Re(r2)x)-δ2B3sinh(Re(r2)x)
-δ3B6cos(Im(r2)x)-δ3B5sin(Im(r2)x)
(2)判别式小于0,μ1<0:实根,μ2,μ3:共轭复根
令r2,r3为共轭复数,Re表示实部,Im表示虚部;
wx(x)=-δ1B2cos(r1x)-δ1B1sin(r1x)-δ2B4cosh(Re(r2)x)-δ2B3sinh(Re(r2)x)
-δ3B6cos(Im(r2)x)-δ3B5sin(Im(r2)x)
(3)判别式大于0,μ1>0:实根,μ2>0:实根,μ3>0:实根
令
(4)判别式大于0,μ1>0:实根,μ2>0:实根,μ3<0:实根
(5)判别式大于0,μ1>0:实根,μ2<0:实根,μ3>0:实根
(6)判别式大于0,μ1<0:实根,μ2>0:实根,μ3>0:实根
(7)判别式大于0,μ1>0:实根,μ2<0:实根,μ3<0:实根
(8)判别式大于0,μ1<0:实根,μ2>0:实根,μ3<0:实根
(9)判别式大于0,μ1<0:实根,μ2<0:实根,μ3>0:实根
(10)判别式大于0,μ1<0:实根,μ2<0:实根,μ3<0:实根
将得到的模态函数表达式代入边界条件,整理以后得到如下形式的方程组:
其中表示线性方程组系数矩阵,A11~A66表示对应位置的系数值;令采用matlab中的非线性优化函数fsolve进行求解得到ωx,再将ωx代入式(7)中得到线性方程组,对线性方程组求解得到B1,B2,…,B6的值,从而得到x向的振动模态函数的表达式。
所述步骤2-3中基于已知x向模态函数计算得到频率值wy和y向模态函数方法如下:
根据常微分方程组Ⅱ,令表示微分算子d/dy,则有如下的式子:
展开式(8)中的行列式得到如下的方程:
其中ψ=wy,or
j1=e3+f3+g3-e1f2-g1f1;
j2=g3f3+e3g3+e3f3+e1f1g2+e2f2g1-e3g1f1-e2g2
j3=e3g3f3-e2g2f3-e1g3f2;
将试验解ψ=eλ代入方程(9),由于表示微分算子,则产生如下的方程:
λ6+j1λ4+j2λ2+j3=0 (10)
令μ=λ2,则
μ3+j1μ2+j2μ+j3=0 (11)
令判别式Δ=18j1j2j3-4j1 3j3+j1 2j2 2-4j2 3-27j3 2,
定义wy(y)与函数表达式中的系数关系数值仍表示为δ1,δ2,δ3,与函数表达式中的系数关系数值仍表示为γ1,γ2,γ3如下
其中i=1,2,3,k为调节系数,为展开式用到的正弦函数频率值,根据展开式不同的项来进行取整数值,m=1,2,…,∞;ri为根据式(11)的解计算得到的值;根据判别式的符号以及式(11)中解μ1,μ2,μ3的情形得到如下10种y向振动模态函数形式:
(11)判别式小于0,μ1>0:实根,μ2,μ3:共轭复根
r2,r3为共轭复数,Re表示实部,Im表示虚部;
wy(y)=-δ1B2cosh(r1y)-δ1B1sinh(r1y)-δ2B4cosh(Re(r2)y)-δ2B3sinh(Re(r2)y)
-δ3B6cos(Im(r2)y)-δ3B5sin(Im(r2)y)
(12)判别式小于0,μ1<0:实根,μ2,μ3:共轭复根
令r2,r3为共轭复数,Re表示实部,Im表示虚部;
wy(y)=-δ1B2cos(r1y)-δ1B1sin(r1y)-δ2B4cosh(Re(r2)y)-δ2B3sinh(Re(r2)y)
-δ3B6cos(Im(r2)y)-δ3B5sin(Im(r2)y)
(13)判别式大于0,μ1>0:实根,μ2>0:实根,μ3>0:实根
令
(14)判别式大于0,μ1>0:实根,μ2>0:实根,μ3<0:实根
(15)判别式大于0,μ1>0:实根,μ2<0:实根,μ3>0:实根
(16)判别式大于0,μ1<0:实根,μ2>0:实根,μ3>0:实根
(17)判别式大于0,μ1>0:实根,μ2<0:实根,μ3<0:实根
(18)判别式大于0,μ1<0:实根,μ2>0:实根,μ3<0:实根
(19)判别式大于0,μ1<0:实根,μ2<0:实根,μ3>0:实根
(20)判别式大于0,μ1<0:实根,μ2<0:实根,μ3<0:实根
将得到的表达式代入边界条件,整理以后得到如下形式的方程组:
令采用matlab中的非线性优化函数fsolve进行求解得到ωy再将ωy代入方程(12)中得到线性方程组,对线性方程组求解得到B1,B2,…,B6的值,从而得到y向的振动模态函数的表达式。
步骤三:利用最终得到的两个方向的振动模态进行叠加得到矩形板振动模态。
具体实例选择铝合金矩形板作为计算对象,具体的参数如下:
板长a=1m,板宽b=1m,板厚h=0.1m,板材料密度ρ=2700kg/m3,杨氏模量E=70×109N/m2,泊松比v=0.3,调节系数k=5/6,D=6.4103×106,G=2.6923×1010。
边界条件设定为y向均为简支边界条件,x向也为简支边界条件。已知y向第一阶模态函数及其导数,需要求解x向一阶模态函数。
本发明的矩形板振动模态计算方法,具体包括如下步骤:
步骤一:建立如图2坐标系,采用变分法和变量分离方法计算得到如下的两个方程组:
步骤二:确定矩形板两个方向的边界条件,在预设定的模态阶数下采用不同方向相互迭代来计算振动频率和振型,直至两个方向计算得到的振动频率误差在某个指定范围时停止。具体步骤如下:
(Ⅰ)设定y(或者x)向为第m阶模态,根据y(或者x)向边界条件选择相应自由振动梁的模态函数作为y(或者x)向的模态函数;
如图所示,图3中模态函数为先确定矩形板y向的边界条件为简支时得到的板带y向振动第一阶模态函数及其前两阶导数。
(Ⅱ)利用给定的y(或者x)向模态函数计算得到x(或者y)向的常微分方程组(Ⅰ)(或者(Ⅱ)),根据常微分方程特征值的情况确定x(或者y)向模态函数的形式,根据边界条件建立以模态函数未知系数为变量的线性方程组,对线性方程组6×6系数矩阵行列式值为0的非线性方程进行优化求解,得到频率值wx,再利用wx代入线性方程组进行求解得到模态函数系数,从而得到x(或者y)向的模态函数。
则在优化时输入初始值20hz时,得到如下的公式计算。
(D6+d1D4+d2D2+d3)ψ=0
其中ψ=wx,or
d1=a3+b3+c3-b1c1-a1c2;
d2=a3b3+a3c3+b3c3+a2b1c2+a1b2c1-a3b1c1-a2b2-a1b3c2;
d3=a3b3c3-a2b2c3。
其中d1=-1027.4,d2=18284.5,d3=281500;
用eλ代入(D6+d1D4+d2D2+d3)ψ=0产生如下的方程:
λ6+d1λ4+d2λ2+d3=0
令μ=λ2,则
μ3+d1μ2+d2μ+d3=0
则判别式Δ=18d1d2d3-4d1 3d3+d1 2d2 2-4d2 3-27d3 2>0,
此时特征根分别为:
μ1=1009,μ2=28.2672,μ3=-9.8696
由于判别式大于0,μ1>0:实根,μ2>0:实根,μ3<0:实根
其中模态函数表达式为:
其中δ1=0,δ2=0.1780,δ3=0.3174,γ1=10.1111,γ2=0.0989,γ3=0.8951。
将得到的表达式代入边界条件,整理以后得到如下形式的方程组:
令
采用matlab中的非线性优化函数进行求解得到ωx。
再将ωx代入式(3)到线性方程组,对线性方程组求解得到B1,B2,…,B6的值,从而得到x向的振动模态函数的表达式。具体的归一化后的x向幅值、绕x轴转角、绕y轴转角的第一阶模态如图4-6所示。
(Ⅲ)利用步骤(Ⅱ)中计算得到的x(或者y)的模态函数,计算得到y(或者x)向的常微分方程组(Ⅱ)(或者(Ⅰ)),根据常微分方程特征值的情况确定y(或者x)向模态函数的形式,根据边界条件建立以模态函数未知系数为变量的线性方程组,利用系数矩阵的行列式为0进行优化求解得到频率值wy,再利用wy代入线性方程组进行求解得到模态函数系数,从而得到y(或者x)向的模态函数。
按照如(Ⅱ)中相同步骤进行计算得到y向一阶模态函数如图7-9所示。
(Ⅳ)比较步骤(Ⅱ)和(Ⅲ)得到的频率值的大小,如果满足|wx-wy|≤ε(ε为误差值,取0.0001),则退出迭代。如果不满足条件,则将步骤(Ⅲ)中得到的y(或者x)向的模态函数作为给定的步骤(Ⅱ)模态函数进行迭代计算。
步骤三:利用最终得到的两个方向的一阶振动模态进行叠加得到矩形板振动模态,如图10-12所示。
Claims (6)
1.一种基于谱有限元的矩形板振动模态计算方法,包括以下步骤:
步骤一:基于能量函数变分原理得到板振动的边界值问题模型,再基于变量分离方法得到两个方向的自由振动方程;
步骤二:确定矩形板两个方向的边界条件,在预设定的模态阶数下采用不同方向相互迭代来计算振动频率和振型,直至两个方向计算得到的振动频率误差在某个指定范围时停止;
步骤三:利用最终得到的两个方向的振动模态进行叠加得到矩形板振动模态。
2.根据权利要求1所述的基于谱有限元的矩形板振动模态计算方法,其特征在于:所述步骤一具体步骤为
建立矩形板直角坐标系(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)=w0(x,y,t);
其中w0(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)ejωt,
其中二元变量W(x,y)表示板(x,y)位置的振动幅值,二元变量表示板(x,y)位置绕x轴的旋转角,二元变量表示板(x,y)位置绕y轴的旋转角;ω:圆频率(rad/s);
求Π的变分且令变分等于0,即
对变量W(x,y),进行二元变量分离有:
W(x,y)=Wx(x)Wy(y)
其中Wx(x),Wy(y)分别表示W(x,y)在x,y方向的分离模态函数,分别表示在x,y方向的分离模态函数,分别表示在x,y方向的分离模态函数
假定已知y向振动模态函数Wy,利用如下变分公式:
δW(x,y)=WyδWx
得到关于x向三个变量的常微分方程组Ⅰ:
其中
假定已知x向振动模态函数Wx,则利用如下变分公式:
δW(x,y)=WxδWy
得到关于y向三个变量的常微分方程组Ⅱ:
其中
3.根据权利要求1所述的基于谱有限元的矩形板振动模态计算方法,其特征在于:所述步骤二中,
针对x方向的边界条件定义为三种:
简支边界:Wx=0,Mxx=0;
夹紧边界:Wx=0,
自由边界:Vxz=0,Mxx=0,Mxy=0;
其中力矩
经过计算后得到如下的计算用边界条件:
简支边界:Wx=0,
夹紧边界:Wx=0,
自由边界:
针对y方向的边界条件定义为三种:
简支边界:Wy=0,Myy=0;
夹紧边界:Wy=0,
自由边界:Vyz=0,Myy=0,Mxy=0;
其中力矩
经过计算后得到如下的计算用边界条件:
简支边界:Wy=0,
夹紧边界:Wy=0,
自由边界:
4.根据权利要求3所述的基于谱有限元的矩形板振动模态计算方法,其特征在于:所述步骤二中,采用不同方向相互迭代计算振动频率和振型函数的具体步骤如下
2-1)设定y向为第m阶模态,根据y向边界条件选择相应自由振动梁的模态函数作为y向的模态函数;
2-2)利用给定的y向模态函数计算得到x向的常微分方程组Ⅰ,根据常微分方程特征值的情况确定x向模态函数的形式,根据边界条件建立以模态函数未知系数为变量的线性方程组,对线性方程组6×6系数矩阵行列式值为0的非线性方程进行优化求解,得到频率值wx,再利用wx代入线性方程组进行求解得到模态函数系数,从而得到x向的模态函数;
2-3)利用步骤2-2中计算得到的x向的模态函数,计算得到y向的常微分方程组Ⅱ,根据常微分方程特征值的情况确定y向模态函数的形式,根据边界条件建立以模态函数未知系数为变量的线性方程组,利用系数矩阵的行列式为0进行优化求解得到频率值wy,再利用wy代入线性方程组进行求解得到模态函数系数,从而得到y向的模态函数;
2-4)比较步骤2-2和2-3得到的频率值的大小,如果满足|wx-wy|≤ε,则退出迭代,ε为误差值,取值为0.0001;如果不满足条件,则将步骤2-3中得到的y向的模态函数作为给定的步骤2-2中模态函数进行迭代计算。
5.根据权利要求4所述的基于谱有限元的矩形板振动模态计算方法,其特征在于:所述步骤2-2中,基于已知y向模态函数计算得到频率值wx和x向的模态函数方法如下:
根据常微分方程组Ⅰ,令表示微分算子d/dx,则有如下的式子:
展开式(3)中的行列式得到如下的方程:
其中ψ=wx,
d1=a3+b3+c3-b1c1-a1c2;
d2=a3b3+a3c3+b3c3+a2b1c2+a1b2c1-a3b1c1-a2b2-a1b3c2;
d3=a3b3c3-a2b2c3;
将试验解ψ=eλ代入方程(4),由于表示微分算子,可产生如下的方程:
λ6+d1λ4+d2λ2+d3=0 (5)
令μ=λ2,则
μ3+d1μ2+d2μ+d3=0 (6)
令判别式Δ=18d1d2d3-4d1 3d3+d1 2d2 2-4d2 3-27d3 2,
定义wx(x)与函数表达式中的系数关系数值δ1,δ2,δ3,与函数表达式中的系数关系数值γ1,γ2,γ3如下
其中i=1,2,3,k为调节系数,为展开式用到的正弦函数频率值,根据展开式不同的项来进行取整数值,m=1,2,…,∞;ri为根据式(6)的解计算得到的值;根据判别式的符号以及式(6)中解μ1,μ2,μ3的情形得到如下10种x向振动模态函数形式:
(1)判别式小于0,μ1>0:实根,μ2,μ3:共轭复根
r2,r3为共轭复数,Re表示实部,Im表示虚部;
wx(x)=-δ1B2cosh(r1x)-δ1B1sinh(r1x)-δ2B4cosh(Re(r2)x)-δ2B3sinh(Re(r2)x)
-δ3B6cos(Im(r2)x)-δ3B5sin(Im(r2)x)
(2)判别式小于0,μ1<0:实根,μ2,μ3:共轭复根
令r2,r3为共轭复数,Re表示实部,Im表示虚部;
wx(x)=-δ1B2cos(r1x)-δ1B1sin(r1x)-δ2B4cosh(Re(r2)x)-δ2B3sinh(Re(r2)x)
-δ3B6cos(Im(r2)x)-δ3B5sin(Im(r2)x)
(3)判别式大于0,μ1>0:实根,μ2>0:实根,μ3>0:实根
令
(4)判别式大于0,μ1>0:实根,μ2>0:实根,μ3<0:实根
(5)判别式大于0,μ1>0:实根,μ2<0:实根,μ3>0:实根
(6)判别式大于0,μ1<0:实根,μ2>0:实根,μ3>0:实根
(7)判别式大于0,μ1>0:实根,μ2<0:实根,μ3<0:实根
(8)判别式大于0,μ1<0:实根,μ2>0:实根,μ3<0:实根
(9)判别式大于0,μ1<0:实根,μ2<0:实根,μ3>0:实根
(10)判别式大于0,μ1<0:实根,μ2<0:实根,μ3<0:实根
将得到的模态函数表达式代入边界条件,整理以后得到如下形式的方程组:
其中表示线性方程组系数矩阵,A11~A66表示对应位置的系数值;令采用matlab中的非线性优化函数fsolve进行求解得到ωx,再将ωx代入式(7)中得到线性方程组,对线性方程组求解得到B1,B2,…,B6的值,从而得到x向的振动模态函数的表达式。
6.根据权利要求5所述的基于谱有限元的矩形板振动模态计算方法,其特征在于:所述步骤2-3中基于已知x向模态函数计算得到频率值wy和y向模态函数方法如下:
根据常微分方程组Ⅱ,令表示微分算子d/dy,则有如下的式子:
展开式(8)中的行列式得到如下的方程:
其中
j1=e3+f3+g3-e1f2-g1f1;
j2=g3f3+e3g3+e3f3+e1f1g2+e2f2g1-e3g1f1-e2g2
j3=e3g3f3-e2g2f3-e1g3f2;
将试验解ψ=eλ代入方程(9),由于表示微分算子,则产生如下的方程:
λ6+j1λ4+j2λ2+j3=0 (10)
令μ=λ2,则
μ3+j1μ2+j2μ+j3=0 (11)
令判别式Δ=18j1j2j3-4j1 3j3+j1 2j2 2-4j2 3-27j3 2,
定义wy(y)与函数表达式中的系数关系数值仍表示为δ1,δ2,δ3,与函数表达式中的系数关系数值仍表示为γ1,γ2,γ3如下
其中i=1,2,3,k为调节系数,为展开式用到的正弦函数频率值,根据展开式不同的项来进行取整数值,m=1,2,…,∞;ri为根据式(11)的解计算得到的值;根据判别式的符号以及式(11)中解μ1,μ2,μ3的情形得到如下10种y向振动模态函数形式:
(11)判别式小于0,μ1>0:实根,μ2,μ3:共轭复根
r2,r3为共轭复数,Re表示实部,Im表示虚部;
wy(y)=-δ1B2cosh(r1y)-δ1B1sinh(r1y)-δ2B4cosh(Re(r2)y)-δ2B3sinh(Re(r2)y)
-δ3B6cos(Im(r2)y)-δ3B5sin(Im(r2)y)
(12)判别式小于0,μ1<0:实根,μ2,μ3:共轭复根
令r2,r3为共轭复数,Re表示实部,Im表示虚部;
wy(y)=-δ1B2cos(r1y)-δ1B1sin(r1y)-δ2B4cosh(Re(r2)y)-δ2B3sinh(Re(r2)y)
-δ3B6cos(Im(r2)y)-δ3B5sin(Im(r2)y)
(13)判别式大于0,μ1>0:实根,μ2>0:实根,μ3>0:实根
令
(14)判别式大于0,μ1>0:实根,μ2>0:实根,μ3<0:实根
(15)判别式大于0,μ1>0:实根,μ2<0:实根,μ3>0:实根
(16)判别式大于0,μ1<0:实根,μ2>0:实根,μ3>0:实根
(17)判别式大于0,μ1>0:实根,μ2<0:实根,μ3<0:实根
(18)判别式大于0,μ1<0:实根,μ2>0:实根,μ3<0:实根
(19)判别式大于0,μ1<0:实根,μ2<0:实根,μ3>0:实根
(20)判别式大于0,μ1<0:实根,μ2<0:实根,μ3<0:实根
将得到的表达式代入边界条件,整理以后得到如下形式的方程组:
令系数矩阵的行列式采用matlab中的非线性优化函数fsolve进行求解得到ωy,再将ωy代入方程(12)中得到线性方程组,对线性方程组求解得到B1,B2,…,B6的值,从而得到y向的振动模态函数的表达式。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710059235.1A CN106777825B (zh) | 2017-01-24 | 2017-01-24 | 一种基于谱有限元的矩形板振动模态计算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710059235.1A CN106777825B (zh) | 2017-01-24 | 2017-01-24 | 一种基于谱有限元的矩形板振动模态计算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106777825A true CN106777825A (zh) | 2017-05-31 |
CN106777825B CN106777825B (zh) | 2020-03-27 |
Family
ID=58943131
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710059235.1A Active CN106777825B (zh) | 2017-01-24 | 2017-01-24 | 一种基于谱有限元的矩形板振动模态计算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN106777825B (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109948180A (zh) * | 2019-01-25 | 2019-06-28 | 北京航空航天大学 | 一种正交各向异性对边简支矩形薄板振动分析方法 |
CN113358308A (zh) * | 2021-06-03 | 2021-09-07 | 哈尔滨工业大学 | 基于有限测点和全局模态的组合结构横向位移确定方法 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20140379314A1 (en) * | 2011-10-31 | 2014-12-25 | Osaka University | Analyzer, analysis method, and analysis program |
CN104778377A (zh) * | 2015-05-04 | 2015-07-15 | 中国矿业大学 | 一种组合梁弯曲振动的固有频率分析方法 |
CN105740541A (zh) * | 2016-01-29 | 2016-07-06 | 厦门大学 | 一种基于结构动力学模型修正的预应力识别方法 |
CN106326501A (zh) * | 2015-06-15 | 2017-01-11 | 上海东浩兰生国际服务贸易(集团)有限公司 | 建筑物结构动力分析用自振频率和振型的计算方法 |
-
2017
- 2017-01-24 CN CN201710059235.1A patent/CN106777825B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20140379314A1 (en) * | 2011-10-31 | 2014-12-25 | Osaka University | Analyzer, analysis method, and analysis program |
CN104778377A (zh) * | 2015-05-04 | 2015-07-15 | 中国矿业大学 | 一种组合梁弯曲振动的固有频率分析方法 |
CN106326501A (zh) * | 2015-06-15 | 2017-01-11 | 上海东浩兰生国际服务贸易(集团)有限公司 | 建筑物结构动力分析用自振频率和振型的计算方法 |
CN105740541A (zh) * | 2016-01-29 | 2016-07-06 | 厦门大学 | 一种基于结构动力学模型修正的预应力识别方法 |
Non-Patent Citations (1)
Title |
---|
曾军才 等: "正交各向异性矩形板的自由振动特性分析", 《振动与冲击》 * |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109948180A (zh) * | 2019-01-25 | 2019-06-28 | 北京航空航天大学 | 一种正交各向异性对边简支矩形薄板振动分析方法 |
CN113358308A (zh) * | 2021-06-03 | 2021-09-07 | 哈尔滨工业大学 | 基于有限测点和全局模态的组合结构横向位移确定方法 |
Also Published As
Publication number | Publication date |
---|---|
CN106777825B (zh) | 2020-03-27 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Lewicka et al. | The Föppl-von Kármán equations for plates with incompatible strains | |
Aref et al. | Vortex dynamics of the two-dimensional turbulent shear layer | |
Sung et al. | Process identification and PID control | |
Reddy Gorla | Heat transfer in an axisymmetric stagnation flow on a cylinder | |
Dawe et al. | Buckling of rectangular Mindlin plates | |
Gupta et al. | An assessment of a non-polynomial based higher order shear and normal deformation theory for vibration response of gradient plates with initial geometric imperfections | |
Janno et al. | Solitary waves in nonlinear microstructured materials | |
CN105205035B (zh) | 一种非均匀弹性约束边界条件矩形板结构面内振动分析方法 | |
Wang et al. | 3D free vibration analysis of multi-directional FGM parallelepipeds using the quadrature element method | |
CN106777825A (zh) | 一种基于谱有限元的矩形板振动模态计算方法 | |
CN110162826B (zh) | 薄壁结构热气动弹性动响应分析方法 | |
CN103838141B (zh) | 一种面向控制的大型天线建模方法 | |
CN110717216A (zh) | 不规则波下带柔性气囊直升机横摇响应预报方法 | |
Spescha | Framework for efficient and accurate simulation of the dynamics of machine tools | |
Zhou et al. | 3-D vibration analysis of skew thick plates using Chebyshev–Ritz method | |
CN104091003A (zh) | 一种基础运动时柔性壳结构大变形响应的有限元建模方法 | |
CN104636556B (zh) | 成任意角度连接的有限尺寸板结构振动响应计算方法 | |
CN109902350A (zh) | 对变截面梁的截面惯性矩进行模型修正中克服模态交换的方法 | |
Smith | Numerical study of two-dimensional stratified turbulence | |
CN107766686A (zh) | 基于matlab计算fgm薄板刚柔耦合动力学响应的仿真方法 | |
Yang | A duality theorem for plastic plates | |
Zhang et al. | Bifurcation behavior and chaotic dynamics of a three-degree-of-freedom aeroelastic system | |
CN108334688B (zh) | 基于matlab的旋转功能梯度厚板动力学行为计算方法 | |
Hassan et al. | Finite elements for the one variable version of Mindlin-Reissner plate | |
Alikakos et al. | Analysis of a corner layer problem in anisotropic interfaces |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |