模态综合法求解悬臂板模态

wei_x

解析解求解：
三边自由一边固支的悬臂板可以利用单向解析行数带入变分方程降为另一方向的常微分方程求解。
具体可以参考曹志远的《板壳振动理论》47-49页
ansys求解：

以长和宽2的正方形悬臂板为例：材料参数：弹性模量E=2.1e11 泊松比g=0.3密度r=7.3e3 厚度t=0.05

网格划分：在长度和宽度方向上均划分20个单元（为了与后面程序一致）

/PREP7    BLC4, , ,2,2    ！前处理 画2×2矩形

MP,DENS,1,7.3e3           ！密度

MP,EX,1,2.1e11            ！弹性模量

MP,PRXY,1,0.3             ！泊松比

ET,1,SHELL63              ！单元类型

R,1,0.05,0.05,0.05,0.05,, ,！实常数，厚度

TYPE,1 MAT,1REAL,1         ！材料及单元类型编号

ESIZE,0.1,0,              ！变长0.1的单元

AMESH,all                 ！画网格

DL,4,,ALL,                ！施加固定约束

FINISH                    ！直接求解命令流 !整体求解

/SOL     ANTYPE,2   MODOPT,LANB,10     EQSLV,SPAR      MODOPT,LANB,20,0,99999999, ,OFF   SOLVE

表一 ANSYS直接求解结果

阶数	1	2	3	4	5	6	7	8	9	10
频率	11.2	27.4	68.7	87.8	100.1	175.2	197.8	207.3	229.5	281.3

wei_x

参考了徐斌，高跃飞，余龙　的《Matlab有限元动力学分析与工程应用》

有限元方法求解悬臂板模态MATLAB程序：

没有包括计算单元刚度和质量的子函数

clear all
clc%清空变量

tic%记时起点

% 定义材料

E=2.1e11;%弹性模量

poisson =0.3;% 泊松比

density=7.3e3;%密度

t=0.05;%板的厚度

lx=2;%x方向长度

ly=2;%y方向长度

jdx=21;%x向节点数

jdy=21;%y向节点数

%初始化

%刚度和质量阵

Tol_dof=3*jdx*（jdy-1）;%在y=0边界上施加固定约束，去掉21个节点

k(1:Tol_dof,1:Tol_dof)=0;%总刚阵

m(1:Tol_dof,1:Tol_dof)=0;%总质量阵

Tol_element=(jdx-1)*(jdy-1); % 总单元数

%节点矩阵

en(1:Tol_element,1:4)=0;%每个单元四个节点，记录节点编号

for ni=1:jdx-1
for nj=1:jdy-1
en(ni+(nj-1)*(jdx-1),1)=ni+(nj-1)*jdx;
en(ni+(nj-1)*(jdx-1),2)=ni+1+(nj-1)*jdx;
en(ni+(nj-1)*(jdx-1),4)=ni+nj*jdx;
en(ni+(nj-1)*(jdx-1),3)=ni+1+nj*jdx;
end

end

%节点位移，约束以及自由度坐标

disp(1:jdx*jdy,1:3)=1;% 位移阵

constraints=1:1:jdx;% y=0上被约束节点

disp(constraints,:)=0;dof=0;%自由度标记

for ni=1:jdx*jdy
for nj=1:3
if disp(ni,nj)~=0%有约束不编号
dof=dof+1;
disp(ni,nj)=dof;
end
end

end

%单元长度

el=lx/(jdx-1);

eh=ly/(jdy-1);

%四节点12自由度 单元刚度和质量矩阵

% 输入材料参数和节点集中长度

[ek,dm]=km(el/2,eh/2,poisson,t,E,density);
%整体刚度矩阵和质量矩阵
tg(1:12)=0; % 标记，单元矩阵维数

forloopi=1:(jdx-1)*(jdy-1) %外循环，单元
for zi=1:4
%每个节点三个自由度
tg((zi-1)*3+1)=disp(en(loopi,zi),1);
tg((zi-1)*3+2)=disp(en(loopi,zi),2);
tg((zi-1)*3+3)=disp(en(loopi,zi),3);
end
for jx=1:12
for jy=1:12
if(tg(jx)*tg(jy)~=0)
k(tg(jx),tg(jy))=k(tg(jx),tg(jy))+ek(jx,jy);
m(tg(jx),tg(jy))=m(tg(jx),tg(jy))+dm(jx,jy);
end
end

end
%总装刚度阵和质量阵

end

%求解特征值问题

[v,d] = eig(k,m);
%eig函数求解特征值问题

tempd=diag(d);
%记录特征值

[nd,sortindex]=sort(tempd);
%对特征值排序

v=v(:,sortindex);
%特征向量排序

mode_number=1:15;
%求解阶数

frequency(mode_number)=sqrt(nd(mode_number))/(2*pi);%记录角频率

toc%记时终点

表三 MATLAB直接求解结果

阶数	1	2	3	4	5	6	7	8	9	10
频率	11.208	27.469	68.764	87.797	99.961	174.77	197.92	207.16	229.18	299.24
相对误差（%）	0.07	0.25	0.09	0.03	0.13	0.24	0.06	0.06	0.13	6

第十阶频率有很大差别，取出ANSYS分析结果的振型分析，可知第十节频率为平行于固支边的弯曲振型，在薄板弯曲单元中不包括这个方向上的自由度，因而取ANSYS分析结果中的第十一阶频率，其大小为300.1Hz，与MATLAB结果对比，误差为0.29%

wei_x

模态综合法介绍：
刘明杰 《固定界面模态综合法的动力学原理》
模态综合法主要内容就是两个坐标变换，论文讲得比较清楚：
第一次坐标变换：模态矩阵对每个子结构的刚度阵和质量阵进行缩减；
第二次坐标变换：界面处节点的位移必须相等，实现子结构间的结合。

wei_x

ANSYS模态综合法求解

ANSYS模态综合法是子结构在动力分析中的应用，其基本过程包括三方面：生成过程、使用过程、扩展过程。ANSYS提供了友好的模态综合法（ComponentMode Synthesis）CMS向导（Preprocessor>Modeling>CMS），可以方便的定义超单元和交界面，而且可以对模态综合法分析生成的文件进行管理和组织。

1.创建超单元

2.选择CMS求解方法(CMSOPT)：固定界面法或自由界面法。

3.命名超单元矩阵文件(SEOPT)。

4.对于考虑阻尼时，输入集中质量矩阵公式。

5.定义主自由度：在超单元的交界面定义主自由度。

6.保持数据库：这是必须做的，因为在扩展模态时需要有相同的数据库文件。

7.求解生成超单元矩阵文件。

CMS的使用和扩展部分与子结构基本相同，但是CMS只支持模态分析。在子结构分析中，需要对整体结构中的每一个子结构进行生成和扩展，然后将各个子结构集合成整个模型。在自由界面模态综合法中，使用MODOPT扩展模态数，应小于每个超单元求解的模态数。若需要扩展更多的模态，需要CMS的超单元重新求解更多的阶数。

wei_x

还是一楼的例子，接着分析：

需要先建立三个集合，两个单元集合以及一个interface节点集合

写成apdl命令流比较麻烦用GUI操作方便些。

/filnam,part1             ！文件名 part1 子结构A

/solu  antype,substr         ！分析类型 子结构

seopt,part1,2            ！子结构A

cmsopt,fix,20            ！固定界面法求前20阶频率

cmsel,s,part1             ！选定part1单元集合

cmsel,s,interface          ！interface节点集合

m,all,all                   ！选择所有

nsle solve finish save       ！求解并保存

wei_x

/filnam,part2                ！文件名 part2 子结构B

/solu     antype,substr       ！分析类型 子结构

seopt,part2,2               ！子结构B

cmsopt,fix,20               ！固定界面法求前20阶频率

cmsel,s,part2                ！选定part2单元集合

cmsel,s,interface             ！interface节点集合

m,all,all                      ！选择所有

nsle solve finish save          ！求解并保存

wei_x

固定界面模态综合法求解命令流：

/clear,nostart      ！清空变量，启动一个新分析         
/filnam,use       ！文件名 use 整体体结构

/prep7            

et,1,matrix50    ！定义超单元     type,1      se,part1    ！选择子结构一       se,part2   ！选择子结构一

finish

/solu
antype,modal    ！分析类型        modopt,lanb,10 ！求解前10阶频率     solve

finish

save

表二 ANSYS模态综合法求解结果

子结构A
阶数	1	2	3	4	5	6	7	8	9	10
频率	287.5	302.3	355.3	458.9	624.8	792.3	812.7	857.9	879.8	995.8
子结构B
阶数	1	2	3	4	5	6	7	8	9	10
频率	45.1	69.2	131.8	247.1	282.1	318.7	407.1	439.6	557.5	682.2
整体结构
阶数	1	2	3	4	5	6	7	8	9	10
频率	11.2	27.5	68.7	87.8	100.1	175.2	197.8	207.3	229.5	281.5

与表一中直接求解结构对比，可以看出本例中，ANSYS模态综合法具有较高的精度（保留一位小数），前10阶频率基本上是一致的。在计算大型复杂结构，模态综合法可以节省大量的计算时间和计算机资源，提高效率；同时可以灵活修改大系统的子系统设计。因此，对于复杂大型结构，如飞机、车辆、船舶、高层建筑等结构，采用ANSYS模态综合法来对结构进行模态分析，可以在精度和计算速度上得到较好的解决方案。

wei_x

子结构的划分与ANSYS中处理一致，主模态分析程序如下：

子结构A（只列出修改部分）：

lx=2;
%x方向长度 ly=1; %y方向长度

jdx=21; %x向节点数
jdy=11; %y向节点数

%y=0和1沿x轴去掉42个点，得到总自由度

Tol_dof=3*jdx*(jdy-2);

cons_Ori=1:1:jdx;
%原有约束节点

% 固定界面法人为施加约束

cons_Arti=jdx*(jdy-1)+1:1:jdx*jdy;

%总约束

constraints=[cons_Ori,cons_Arti];

master_mode=15;add_num=jdx*3;

% 取前master_mode阶振型作为模态综合

fi_unres1=cat(1,v(1:master_mode,:)',zeros(add_num,master_mode));

子结构B（只列出修改部分）：

lx=2;
%x方向长度
ly=1;
%y方向长度

jdx=21; %x向节点数
jdy=11; %y向节点数

%总自由度

Tol_dof=3*jdx*(jdy-1);

%固定界面法人为施加约束

cons_Arti=1:1:jdx;

% 单边约束节点编号

constraints=cons_Arti; %总约束

master_mode=15;add_num=jdx*3;

fi_unres2=cat(1,zeros(add_num,master_mode),v(1:master_mode,:)');

% 主模态中，认为施加固定约束的节点，其位移为零，补充完整加入模态综合

表四 MATLAB子结构求解结果

子结构A
阶数	1	2	3	4	5	6	7	8	9	10
频率	287.8	302.1	353.9	455.9	620.3	794.8	813.7	852.3	877.7	988.8
子结构B
阶数	1	2	3	4	5	6	7	8	9	10
频率	45.1	69.1	131.5	246.4	282.3	318.5	405.4	438.8	553.6	681.2

wei_x

约束模态分析程序如下：

子结构一（只列出修改部分）：

lx=2;
%x方向长度 ly=1; %y方向长度

jdx=21; %x向节点数
jdy=11; %y向节点数

Tol_dof=3*jdx*(jdy-1);

cons_Ori=1:1:jdx;
%原有约束节点

constraints=cons_Ori; %总约束

% 总刚度矩阵中按是否在界面上分块

add_num=jdx*3;master_mode=15;

kii=k(1:end-add_num,1:end-add_num);

kij=k(1:end-add_num,end-add_num+1:end);

%提取约束模态

fi_res1=cat(1,-inv(kii)*kij,eye(add_num));

%合并总模态

fi_1=cat(2,fi_unres1,fi_res1);

%用模态对刚阵和质量阵做第一次坐标变换

k1=fi_1'*k*fi_1; m1=fi_1'*m*fi_1;

子结构二（只列出修改部分）：

lx=2;
%x方向长度
ly=1;
%y方向长度

jdx=21; %x向节点数
jdy=11; %y向节点数

%总自由度

Tol_dof=3*jdx*jdy; %子结构没有约束

% 总刚度矩阵中按是否在界面上分块

add_num=jdx*3;master_mode=15;

kii=k(add_num+1:end,add_num+1:end);

kij=k(add_num+1:end,1:add_num);

%提取约束模态

fi_res2=cat(1,eye(add_num),-inv(kii)*kij);

%合并总模态

fi_2=cat(2,fi_res2,fi_unres2);

%用模态对刚阵和质量阵做第一次坐标变换

k2=fi_2'*k*fi_2;
m2=fi_2'*m*fi_2;

模态综合程序如下：

master_mode=15;add_num=jdx*3;Tol_num=length(k1);

%master_mode为主模态保留阶数，add_num为界面上自由度数，Tol_num为变换后刚度矩阵维数

K=zeros(2*Tol_num);M=zeros(2*Tol_num); %初始化模态综合矩阵

K(1:Tol_num,1:Tol_num)=k1;K(Tol_num+1:end,Tol_num+1:end)=k2;

M(1:Tol_num,1:Tol_num)=m1;M(Tol_num+1:end,Tol_num+1:end)=m2;

%非耦合模态综合矩阵，直接写成矩阵形式

beta=zeros(2*Tol_num,2*master_mode+add_num);

beta(1:master_mode,1:master_mode)=eye(master_mode);

beta(master_mode+1:master_mode+add_num,master_mode+1:master_mode+add_num)=eye(add_num);

beta(master_mode+1+add_num:master_mode+2*add_num,master_mode+1:master_mode+add_num)=eye(add_num);

beta(end-master_mode+1:end,end-master_mode+1:end)=eye(master_mode);

%第二次坐标变换矩阵beta，由位移列车很容易求的

KK=beta'*K*beta;MM=beta'*M*beta;%做第二次坐标变换

[v,d] = eig(KK,MM);tempd=diag(d);[nd,sortindex]=sort(tempd);v=v(:,sortindex);

%求解特征值问题

mode_number=1:5;frequency_sum(mode_number)=sqrt(nd(mode_number))/(2*pi);

%得到最终模态综合方法求出频率

MATLAB模态综合法求解结果

阶数	1	2	3	4	5
频率	11.5	29.1	87.7	106.8	152.5

wei_x

matlab求解的结果不是很好，只有基频还过得去，不明白是因为固定界面模态综合法本身的原因还是程序的原因。
终于写完了，嘎嘎。

Rainyboy

回复 wei_x 的帖子

强烈建议楼主整理一下楼层，合并成一个大楼，我个人觉得很有精华的潜质哦~

16443

wei_x 发表于 2010-10-10 20:22

                               
登录/注册后可看大图

参考了徐斌，高跃飞，余龙　的《Matlab有限元动力学分析与工程应用》
有限元方法求解悬臂板模态MATLAB程序 ...

有限元第10阶与解析的第10阶结果不一致这个正常，属于离散化带来的自由度问题。

16443

wei_x 发表于 2010-10-10 20:28

                               
登录/注册后可看大图

模态综合法介绍：
刘明杰 《固定界面模态综合法的动力学原理》
模态综合法主要内容就是两个坐标变换，论文 ...

在耦合的时候只保证了位移平衡，是否考虑过力平衡 ？

wei_x

16443 发表于 2010-10-24 14:48

                               
登录/注册后可看大图

有限元第10阶与解析的第10阶结果不一致这个正常，属于离散化带来的自由度问题。

级数表示的解析解，要用到很多系数，偷懒就改用更方便的ansys建模计算了。
第10阶频率有较大误差是因为我用的板单元没有考虑到z方向的自由度，用薄板理论简化了。
ansys的壳单元自由度比我用的板单元多了z方向的自由度

wei_x

16443 发表于 2010-10-24 14:50

                               
登录/注册后可看大图

在耦合的时候只保证了位移平衡，是否考虑过力平衡 ？

界面是人为强加给结构的，实际情况下并不存在，故耦合的时候只考虑了位移平衡，
界面上应力的是否平衡没有验证过，这可能就是固定界模态综合法算出来的结果只能保证前几阶频率正确的原因，谢谢你的建议。
ANSYS中模态综合法用的是超单元法，没找到相关的资料，论坛中有哪位大虾用过么？