matlab 材料力学计算程序

1．1 基本方法和模型的建立

在小变形条件下，根据力的独立作用原理，无

论载荷多么复杂，都可以将其分解为若干简单载

1 复杂载荷作用F的简支梁

然后应用叠加法得到复杂载荷下梁内的弯矩：如

图1所示为受任意载荷的简支梁。在集中力偶、

集中力和分布载荷单独作用下的弯矩方程分别为

MM( )= RMA-丁+M( 一＆)

R ：一 M (1)

Mp( )= R^P —P( —b)

= 一 P

MfJ( )=R 一 1 q( —c)生二【】： ± 2 L (3)

则在集中力偶M 、集中力P和分布载荷q共同作

用下的弯矩方程为

M ( )= MM( )+MfJ( )+M ( ) (4)

以上各式中应用了跳跃函数，其意义如下

f 0 (I『≤ )

一L『 1( 。， ) ( >)

1．2 计算机分析的实现过程

根据上面介绍的计算模型可应用Matlab编

制如下的计算程序

clear；

L=input( L(In)= )．

M=input( M(KNm)= )．

a input( a(In)= )．

P=input( P(KNm)= )．

b=input( b(In)= )．

q=input( q(KN／m)= )；

c=input( c(In)= )．

d=input( d(In)= )．

nd= 3000；

nf=nd+l：

x=linspace(0，L，nf)；

dx=L／nd；

％ * * * * * * * * * * * * * * * * *

RMA= M／L；nl=a／dx+l；

MM1=RMA X(1：n1)：

MⅣI2=RMA X(nl+l：nf)+M ：

MM=[MMl，MM2 J；

％ * * * * * * * * * * * * * * * * * *

nl=b／dx+1；bb=L—b；

RPA=bb／L*P：

M[Pl=RPA*X(1：n1)；

M口f)2=RPA X(nl+l：nf)

P*(x(nl+l：nf)一b)；

MP=[MP1，MP2 J；

％ * * * * * * * * * * * * * * * * * *

nl=c／dx+l：

n2=d／dx+1；

RqA (L一0．5*(c+d))／L q*(d—c)；

Mql=RqA*X(1：n1)；

Mq2=RqA X(nl+l：n2)一

0．5 q (X(nl+l：n2)一禅祥c)． 2；

Mq3=RqA X(n2+l：nf)一

0．5 q*(X(n2+l：nf) c)．陆世 2

+0．5 q*(X(n2+l：nf)一d)． 2；

Mq=[Mql，Mq2，Mq3]；

‘J／n* * * * * * * * * * * * * * * * *

M = MM +MP+Mq：

subplot(2，l，1)；Mmax=max(M)，Mmin=

rain(M)

plot(X，M)，

title(’复杂早袭肢载荷作用下的弯矩图’)

grid

当L =3 In、M =15 kNm、q=30 kN／m、“

= 0．5 nl、6= l In、f=1．5 In、d =2．5 In时，运

行程序时，得到如图2所示的弯矩图，最大和最小

弯矩分别为Mmax=33．333 kN、Mmin=0=

代码的解释（供参考，下列是根据枣大经验大概的凳袭竖判断，如不对请谅解）

%自定义函数 spacemodel（）作用是当flag某值，执行相当于的程序命令

function [sys,x0,str,ts] = spacemodel(t,x,u,flag)

switch flag,

case 0,

[sys,x0,str,ts]=mdlInitializeSizes

case 1,

sys=mdlDerivatives(t,x,u)

case 3,

sys=mdlOutputs(t,x,u)

case {2,4,9}

sys=[]

otherwise

error(['Unhandled flag = ',num2str(flag)])

end

%自定义函数 mdlInitializeSizes（）作用是当flag=0，执行禅做相当于的程序命令

function [sys,x0,str,ts]=mdlInitializeSizes

sizes = simsizes

sizes.NumContStates = 4

sizes.NumDiscStates = 0圆盘

sizes.NumOutputs = 4输出

sizes.NumInputs = 2输入

sizes.DirFeedthrough = 0反馈

sizes.NumSampleTimes = 1采样时间

sys = simsizes(sizes)大小

x0 = [3001] 初值

str = []空矩阵

ts = [0 0] 初值

%自定义函数 mdlDerivatives（）作用是（flag=1）根据t,x,u的值，将t,x,u变量值赋值给相对应的sys变量

function sys=mdlDerivatives(t,x,u)

%a=1

a=1000*rands(1)

d1=a*0.3*sin(3*pi*t)随机作用力

d2=a*0.1*(1-exp(-t))随机作用力

dq1=x(1)%Angle1 speed:dq1 角速度

dq2=x(2)%Angle2 speed:dq2 角速度

q1=x(3)%Angle1:q1 角度

q2=x(4)%Angle2:q2 角度

tol1=u(1)

tol2=u(2)

%%%%%%%%%%%%%%%%%%%%%%%%% From aslpd_con_s.m

m1=10m2=5 %质量

l1=1l2=0.5%长度

r1=0.5r2=0.25%半径

i1=0.83+m1*r1^2+m2*l1^2%惯性半径

i2=0.3+m2*r2^2% %惯性半径

g=9.8%重力加速度

D=[i1+i2+2*m2*r2*l1*cos(q2) i2+m2*r2*l1*cos(q2)i2+m2*r2*l1*cos(q2) i2]质量

C=[-m2*r2*l1*dq2*sin(q2),-m2*r2*l1*(dq1+dq2)*sin(q2)m2*r2*l1*dq1*sin(q2),0]阻尼

G=[(m1*r1+m2*l1)*g*cos(q1)+m2*r2*g*cos(q1+q2)m2*r2*g*cos(q1+q2)]重量

%%%%%%%%%%%%%%%%%%%%%%%%%%

D2=inv(D)D的逆矩阵

Ta=[d1d2]

A=-D2*C振幅

Z=-D2*G相位角

sys(1)=A(1,1)*x(1)+A(1,2)*x(2)+Z(1)+D2(1,1)*(-Ta(1)+tol1)+D2(1,2)*(-Ta(2)+tol2)动载荷

sys(2)=A(2,1)*x(1)+A(2,2)*x(2)+Z(2)+D2(2,1)*(-Ta(1)+tol1)+D2(2,2)*(-Ta(2)+tol2)动载荷

sys(3)=x(1) 角速度

sys(4)=x(2) 角速度

%自定义函数 mdlOutputs（）作用是（flag=3）根据t,x,u的值，将x变量值赋值给相对应的sys变量

function sys=mdlOutputs(t,x,u)

sys(1)=x(1)%Angle1 speed:dq1 角速度

sys(2)=x(2)%Angle2 speed:dq2 角速度

sys(3)=x(3)%Angle1:q1 角度

sys(4)=x(4)%Angle2:q2 角度

---