load handel.mat
hfile = 'handel.wav'
wavwrite(y, Fs, hfile)
% 读入wav音频文件
[x,fs,bits]=wavread('handel.wav')
figure
% 显示输入信号样本
subplot(3,1,1)
plot(x)xlabel('number of sample')title('Inputput signals')[N col]=size(x)
Px=zeros(size(x))
Py=zeros(size(x))
y=zeros(size(x))
g=zeros(size(x))
% 参数设置
a=0.2
mu=0.5
Pref=0.4
Pmin=0.000000001
Pmax=0.5g(1,1)=1
y(1,1)=x(1,1)
Py(1,1)=(1-a)*x(1,1)*x(1,1)
Px(1,1)=(1-a)*x(1,1)*x(1,1)for n=2:N
% 计算公式
Px(n,1)=a*Px(n-1,1)+(1-a)*x(n,1)*x(n,1)
g(n,1)=g(n-1,1)*(1+mu*Px(n,1)*(Pref-Py(n-1,1)))
y(n,1)=g(n,1)*x(n,1)
if y(n,1)>2
y(n,1)=2
end
if y(n,1)<-2
y(n,1)=-2
end
Py(n,1)=g(n,1)*g(n,1)*Px(n,1)
end
% 显示输出信号样本
subplot(3,1,2)
plot(y)xlabel('number of sample')
title('Output signals')
subplot(3,1,3)
plot(g)xlabel('number of sample')title('gain')legend('gain')
figureplot(y,'r')hold onplot(x)
title('Output signals &Input signals ')xlabel('number of sample')
legend('Output signals','Input signals')
% 播放输出信号
sound(y)
你提出的问题我之前刚好做过,使用有限元方法来进行桁架结构分析。
Matlab编程实现平面杆单元分析
首先,明确Matlab程序要实现的5个重要模块分别为:单元刚度矩阵的求解、单元组装、节点位移的求解、单元应力的求解、节点力的求解。下面给出这5个模块的实现。
1. 单元刚度矩阵求解
定义函数Bar2D2Node_Stiffness,该函数计算单元的刚度矩阵,输入d性模量E,横截面积A,两个节点坐标输出单元刚度矩阵k(4X4)。具体代码如下:
function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2)
L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
x=acos((x2-x1)/L)
C=cos(x)
S=sin(x)
k=E*A/L*[C*C C*S -C*C -C*S C*S S*S -C*S -S*S
-C*C -C*S C*C C*S -C*S -S*S C*S S*S]
2. 单元组装
定义函数Bar2D2Node_Assembly,该函数进行单元刚度矩阵的组装,输入单元刚度矩阵k,单元的节点编号i、j。输出整体刚度矩阵KK,具体代码如下:
function z = Bar2D2Node_Assembly(KK,k,i,j)
DOF(1)=2*i-1
DOF(2)=2*i
DOF(3)=2*j-1
DOF(4)=2*j
for n1=1:4
for n2=1:4
KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2)
end
end
z=KK
3. 节点位移的求解
定义函数Bar2D2Node_Disp(KK,num,p),该函数输入KK为总体刚度矩阵;num为活动自由度编号数组;p为活动自由度方向上的节点力;输出节点位移列阵。具体代码如下:
function u = Bar2D2Node_Disp(KK,num,p)
k=KK(num,num)
u=k\p
4. 单元应力的求解
定义函数函数Bar2D2Node_Stress(E,x1,y1,x2,y2,u),该函数计算单元的应力输入d性模量E,第一个节点坐标(x1,y1),第二个节点坐标(x2,y2)单元节点位移矢量u,返回单元应力标量stress 。具体代码如下:
L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
x=acos((x2-x1)/L)
C=cos(x)
S=sin(x)
stress=E/L*[-C -S C S]*u
5. 计算节点力
定义函数Bar2D2Node_Forces(KK,q),该函数用于计算节点力,KK为刚度矩阵,q为节点位移阵列
function P= Bar2D2Node_Forces(KK,q)
q=zeros(8,1)
q(num)=u
P=KK*q
至此,基于Matlab的杆单元有限元分析的程序设计已经完成,遇到实际问题时可以直接调用这些函数就可以解决问。
经典算例
如图所示的结构,各个杆的d性模量和横截面积都为 , 。试基于MATLAB平台求解该结构的节点位移、单元应力以及支反力。
四杆桁架结构
解答:对该问题进行有限元分析的过程如下
(1) 结构的离散化与编号
对该结构进行自然离散,节点编号和单元编号如上图所示
(2)计算各单元的刚度矩阵(基于国际标准单位)
输入d性模量E、横截面积A,各点坐标。然后分别针对单元1,2,3和4,调用4次Bar2D2Node_Stiffness,就可以得到单元的刚度矩阵。
对应的主程序中代码:
E=2.95e11A=0.0001x1=0y1=0x2=0.4y2=0x3=0.4y3=0.3x4=0y4=0.3
k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2)
k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3)
k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3)
k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3)
(3) 建立整体刚度方程
由于该结构共有4个节点,因此,设置结构总的刚度矩阵为KK(8×8),先对KK清零,然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装。相关主程序代码为:
KK=zeros(8,8)
KK=Bar2D2Node_Assembly (KK,k1,1,2)
KK=Bar2D2Node_Assembly (KK,k2,2,3)
KK=Bar2D2Node_Assembly (KK,k3,1,3)
KK=Bar2D2Node_Assembly (KK,k4,4,3)
(4)边界条件的处理及刚度方程的求解
由图可以看出,被约束的自由度有:节点1的x,y方向自由度,节点2的y方向自由度,4节点的x、y方向两个自由度。则活动自由度编号为3,5,6.活动自由度对应的节点载荷F3=20000N,F5=0N,F6=25000N,采用高斯消去法进行求解,对应的代码为:
num=[3,5,6]%可活动的自由度编号
p=[200000-25000]
u=Bar2D2Node_Disp(KK,num,p)
(5)支反力的计算
在得到整个结构的节点位移后,由原整体刚度方程就可以计算出对应的支反力。这部分对应的主程序的代码如下:
q=zeros(8,1)
q(num)=u%节点位移阵列
P=Bar2D2Node_Forces(KK,q)
(6)单元应力的计算
先从整体位移列阵q中提取出单元的位移列阵,然后,调用计算单元应力的函数Bar2D2Node_ElementStress,就可以得到各个单元的应力分量。
u1=[q(1)q(2)q(3)q(4)]
stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1)
u2=[q(3)q(4)q(5)q(6)]
stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2)
u3=[q(1)q(2)q(5)q(6)]
stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3)
u4=[q(7)q(8)q(5)q(6)]
stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)
(7)计算结果的整理
通过主程序的运行得计算结果。
主程序
%计算各单元的刚度矩阵(以国际标准单位)
E=2.95e11
A=0.0001
x1=0
y1=0
x2=0.4
y2=0
x3=0.4
y3=0.3
x4=0
y4=0.3
k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2)
k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3)
k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3)
k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3)
%建立整体刚度方程
%由于该结构共有4个节点,因此,结构总的刚度矩阵为KK(8×8),先对K清零,然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装。
KK=zeros(8,8)
KK=Bar2D2Node_Assembly (KK,k1,1,2)
KK=Bar2D2Node_Assembly (KK,k2,2,3)
KK=Bar2D2Node_Assembly (KK,k3,1,3)
KK=Bar2D2Node_Assembly (KK,k4,4,3)
%边界条件的处理及刚度方程求解
num=[3,5,6]%可活动的自由度编号
p=[200000-25000]
u=Bar2D2Node_Disp(KK,num,p)
%支反力的计算
q=zeros(8,1)
q(num)=u%节点位移阵列
P=Bar2D2Node_Forces(KK,q)
%各单元的应力计算
u1=[q(1)q(2)q(3)q(4)]
stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1)
u2=[q(3)q(4)q(5)q(6)]
stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2)
u3=[q(1)q(2)q(5)q(6)]
stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3)
u4=[q(7)q(8)q(5)q(6)]
stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)
说的可能有些罗嗦,注意其中有5个function,和最后一个主程序,计算的时候直接运行主程序就可以了。希望能帮助到你。
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)