模态综合法实例:含ansys命令流和matlab源程序

  • A+
所属分类:学霸天地

解析解说明

三边自由一边固支的悬臂板可以利用单向解析行数带入变分方程降为另一方向的常微分方程求解。

解析解求解过程:具体可以参考曹志远的《板壳振动理论》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直接求解结果

模态综合法实例:含ansys命令流和matlab源程序

模态综合法实例:含ansys命令流和matlab源程序

以上参考了徐斌,高跃飞,余龙 的《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直接求解结果

 

模态综合法实例:含ansys命令流和matlab源程序

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

采用模态综合法求解

模态综合法的基本原理可以参考刘明杰发表的《固定界面模态综合法的动力学原理》一文,该文发表于《浙江大学学报(工学版)》1985年06期。

 

模态综合法的核心是两次坐标变换,论文讲得比较清楚:

第一次坐标变换:模态矩阵对每个子结构的刚度阵和质量阵进行缩减;

第二次坐标变换:界面处节点的位移必须相等,实现子结构间的结合。

ANSYS模态综合法求解

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

1.创建超单元

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

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

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

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

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

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

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

还是一楼的例子,接着分析:

 

需要先建立三个集合,两个单元集合以及一个interface节点集合;写成apdl命令流比较麻烦用GUI操作方便些。

 

子结构A:

/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

!求解并保存

模态综合法实例:含ansys命令流和matlab源程序

子结构B:

/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

!求解并保存

模态综合法实例:含ansys命令流和matlab源程序

 

 

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

 

/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模态综合法求解结果

模态综合法实例:含ansys命令流和matlab源程序

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

主模态分析Matlab程序:

子结构的划分与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子结构求解结果

模态综合法实例:含ansys命令流和matlab源程序

 

约束模态分析程序如下:

 

子结构一(只列出修改部分):

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;

%master_mode为主模态保留阶数

add_num=jdx*3;

%add_num为界面上自由度数

Tol_num=length(k1);

%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模态综合法求解结果

模态综合法实例:含ansys命令流和matlab源程序

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

始发于微信公众号:声振之家

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: