有一圆形的等截面销轴.doc

系统建模与仿真实验报告

姓名 孙阳
班级 流体11-1
学号 1107240117

一、MATLAB在优化设计中的应用
有一圆形的等截面销轴,一端固定在机架上,另一端受集中载荷P=10kN和转矩M=100N.m的作用,其简化模型如图1所示。由于结构的需要,轴的长度l󰀁8cm,已知销轴材料的许用弯曲应力[󰀁w]=120MPa;许用切应力[󰀁]=80MPa;允许挠度[f]=0.01cm;密度󰀁=7.8t/m3;弹性模量E=2󰀁105MPa。现要求设计这根销轴,在满足使用要求的前提下,使其质量最轻[5]。
1.1数学模型的建立 优化设计追求的指标为销轴质量Q最轻,Q的计算式为: Q(X)=Q(x1,x2)= 1 4 󰀁x21x2󰀁=14󰀁d2l󰀁=0.00613d2l式中:设计变量X=[x1,x2]T=[d,l]T。 显然,d、l越小,Q值越小。但二者的取值要受强度、刚度 及结构尺寸等各种使用要求的限制,将以上各限制条件整理后可得: (1)抗弯强度条件: 󰀁max= Pl 0.1d 3 󰀁[󰀁w]󰀁󰀁d3-8.33l󰀁0(2)抗扭强度条件: 󰀁max= M 0.2d3 󰀁[󰀁]󰀁󰀁d3-6.25󰀁0(3)刚度条件: fmax= Pl33EJ=64Pl3 3E󰀁d4 󰀁[f]󰀁󰀁d4-0.34l3󰀁0(4)结构尺寸限制: l󰀁lmin󰀁󰀁l-8󰀁0 据此,可写出优化设计的数学模型为: minQ(X)=0.00613d2l󰀁󰀁X=[x1,x2]=[d,l]Ts.t󰀁d3-8.33l󰀁0 d3-6.25󰀁0d4-0.34l3󰀁0l-8󰀁0 这是一个单目标非线性二维约束优化问题。 1.1.2󰀁求优结果 为了便于对采用不同手段得到的优化结果进行比较,以及确认MATLAB分析结果的可信度,特意选择了已经采用其它方法得到最优结果的例子(引自文献[5],并假定其结果正确),最优方案为:d*=4.309cm,l*=8cm,Q*=0.911kg。1.2󰀁采用MATLAB软件对优化问题示例的分析1.2.1󰀁约束函数的性态分析 应用MATLAB语言编程: [x,y]=meshgrid(linspace(0,17,19),linspace(0,11,19));󰀁%根据函数的定义划分网格区域 Q=0.00613*y.^2*x; mesh(x,y,Q);󰀁%通过三维网格模拟目标函数图形xlabel(󰀁l󰀁);ylabel(󰀁d󰀁); zlabel(󰀁Q󰀁);󰀁%标注三坐标轴 title(󰀁目标函数的图像󰀁)󰀁%标注图形名称 󰀁收稿日期:2003-09-25;修订日期:2004-01-30 作者简介:王春香(1962-),女,内蒙古包头人,副教授,硕士,主要研究方向:机械优化设计方法及CAD/CAM/CAI技术。

󰀁󰀁运行程序可画出目标函数的三维图形,如图2所示
同理,可作如下编程: ezplot(󰀁x^3-8.33*y󰀁,[0,10,0,20])󰀁%弯曲强度条件的图形holdony=0󰀁0.01󰀁20;x=6.25^(1/3); plot(x,y,󰀁k-󰀁)󰀁%扭转强度条件的图形holdon ezplot(󰀁x^4-0.34*y^3󰀁,[0,10,0,20])󰀁%刚度条件的图形holdonx=0󰀁0.01󰀁10;y=8; plot(x,y,󰀁k-󰀁)󰀁%长度边界条件的图形title(󰀁各约束函数的图像󰀁)text(6.5,11,󰀁可行域󰀁)holdoff 执行程序绘出各约束函数的图形,如图3所示。 从图2可以看出:目标函数的图像规则,即性态好,对于多数优化方法均适用;但由图3可行域的构成分析可知:d3 -6.25󰀁0 (即抗扭强度条件)不是起作用约束,
设计时可不予考虑。
应用MATLAB软件的优化工具箱对优化问题示例求解 a.取设计变量的初值为:x0=[4󰀁10]t 首先,编写目标函数的m文件:Objfun.m,返回x处的函数值f。 functionf=Objfun(x) f=0.00613*x(1)^2*x(2); b.因设计约束含3个非线性约束,故需编写一个描述非线性约束的m文件:NonLinConstr.m function[c,ceq]=NonLinConstr(x)c(1)=-x(1)^3+8.33*x(2);
c(2)
=-x(1)^3+6.25; c(3)=-x(1)^4+0.34*x(2)^3;ceq=[]; c.设置线性约束的系数: A=[0󰀁-1 0󰀁0];b=[-8󰀁0]; d.给定变量的初值,并且调用优化函数: x0=[4󰀁10];A=[0󰀁-1 0󰀁0];b=[-8󰀁0]; 1b=zeros(2,1); options=optimset(󰀁Display󰀁,󰀁iter󰀁,󰀁LargeScale󰀁,󰀁off󰀁); [x,fval,exitflag,output]=fmincon(󰀁Objfun󰀁,x0,A,b,[],[],lb,[],󰀁NonLinConstr󰀁,options) e.计算结果: x=󰀁󰀁󰀁󰀁󰀁󰀁󰀁%最后的优化结果:󰀁4.0543󰀁8.0000%d*=4.0543,l*=8cmfval=%优化后的最小质量0.8061%Qmin=0.8061kgexitflag=%算法退出处条件 󰀁󰀁󰀁󰀁1output= 󰀁󰀁󰀁iterations󰀁7󰀁󰀁󰀁%函数调用次数󰀁󰀁󰀁funcCount󰀁29%函数评价次数󰀁󰀁󰀁stepsize󰀁1%步长 󰀁󰀁󰀁Algorithm󰀁󰀂medium-scale:SQP,Quasi-Newton,line- search󰀁󰀂%算法 󰀁󰀁󰀁firstorderopt󰀁[]󰀁󰀁󰀁cgiterations󰀁[] 由上可知,最优方案为:d*=4.0543cm,l*=8cm,Q=0.8061kg 1.2.3󰀁优化结果的分析与比较 采用MATLAB优化工具箱的求优结果与原优化方案的比较,请参见表1。 比较2种优化结果表明:MATLAB的优化方案不仅是可信的,而且可使销轴的质量更小(即优化数据朝着对设计有利的方向变化,这与作者采用MATLAB所进行的其它已知优化案例验证分析所得到的结论趋势一致),减轻质量的比率󰀁为: 󰀁= 原方案质量-现优化质量现优化质量=0.911-0.8601 0.8601 =13% d/cm l/cmQ/kg原来的优化方案 4.30980.911MATLAB的优化结果 4.0543 8 0.8061 1.2.4󰀁MATLAB的优化求解步骤 归纳总结使用MATLAB软件求解如上示例(以及综合所做的其它大量的优化题目)的过程和经验分析后得知,求解时的一般步骤: (1)判断优化问题的类型。分析时要区分:󰀁单目标与多目标问题;󰀁线性与非线性问题;󰀁是否为线性规划问题等几种情况。 (2)根据优化问题的类型来选定优化函数。例如,本例属 53 2004年7月王春香,等:MATLAB软件在机械优化设计中的应用研究
于单目标多变量非线性约束优化问题,故选定fmincon优化模块。 (3)为优化模块fmincon提供输入参数。比如初值、等式约束、不等式约束、变量的上下限,编制目标函数的程序(也可以用inline函数直接定义,不编制函数模块),对于含有非线性约束的问题,需要编制非线性约束的函数模块。 (4)根据目标函数的性态,预设优化选项,即options的设置。 (5)在所有的输入参数定义后,调用优化函数进行优化程序调试。 (6)根据优化过程的具体提示信息,修改优化选项的设置,直到达到满足优化函数fmincon所需的优化条件为止。(7)对所得优化数据和结果进行分析和评价(图表的绘制、数据的曲线模拟及可信度分析等)。

二、MATLAB在控制系统中的应用
借助MATLAB及其控制系统工具箱和SIMULINK的强大功能,可以迅速解决大量的数字计算、特殊图形绘制、动态系统建模、系统仿真和分析󰀁下面以几个典型问题来说明其在控制系统中的应用󰀁 2󰀁1根轨迹的绘制 用手工绘制根轨迹是一件很麻烦的事,所绘图形是一个大致形状,而且不是很美观󰀁利用MAT-LAB绘制就可以克服这些缺点󰀁下面用MATLAB绘制根轨迹[2]: 已知传递函数:GH(s)=K(s+8)/[s(s+2)(s2 +8s+32)],绘制根轨迹并求下面两种情形下的K值:(i)两条分支进入右半平面(ii)两条分支从复数极点出发在实轴相交󰀁 解决以上问题很简单,只需几条语句即可完成,源程序如下: den=conv([120],[1832])󰀁%分母有四个极点G=tf([18],den)󰀁%定义分子并创建G(s)Gzp=zpk(G)󰀁%显示G(s)的零点、极点和增益rlocus(G)󰀁󰀁%计算根轨迹 axis([-155-1010])󰀁%调整绘制区域 [kk,clroots]=rlocfind(G)󰀁%计算增益值和极点 MATLAB绘制的根轨迹如图1所示画出了根轨迹后,我们可以交互的使用rlocfind命令来确
定用户点击鼠标所选的根轨迹上任意点所对应的K值,K值所对应的所有闭环极点值也可以使用形如[kk,clroots]=rlocfind(G)的命令来显示󰀁本例中所要求的K值利用以上方法很容易就可以求出:K󰀁50时,从实数极点上出发的两条分支穿入右半平面;k󰀁2070时,从复数极点出发的两条分支到达实轴
󰀁 2󰀁2󰀁阶跃性能分析 利用MATLAB来计算系统性能指标,了解控制器的参数对其性能的影响也是一件很容易的事情.考虑图2中的反馈系统[2],其中Gp(s)=1/[(s+0.01)(s+1)(s+20)],Gc(s)=50,H(s)=1,当d(t)=0时,求单位阶跃输入作用下系统的响应曲线,并求出MO,tp,tr,ts2和ess󰀁再对r(t)=0时由单位阶跃扰动作用引起的响应曲线,并求稳态值󰀁解决本题的源程序如下: denGp=conv(conv([10.01],[11],[120])󰀁󰀁󰀁󰀁󰀁%对象模型Gp=tf(1,denGp) Gcp=50*Gp󰀁%控制器和对象 Tr=feedback(Gcp,1.0)󰀁%单位反馈t=[0󰀁0.05󰀁20]󰀁%
时间向量 yr=step(Tr,t)󰀁%闭环阶跃响应 figure(1);plot(t,yr);grid󰀁%绘制响应曲线[M0,tp,tr,ts2,ess]=tstats(t,yr,1)󰀁%性能指󰀁󰀁󰀁󰀁󰀁标 [numTr,denTr]=tfdata(Tr,󰀁v󰀁)󰀁%提取分子󰀁󰀁󰀁󰀁󰀁分母多项式 [resS,polS,otherS]=residue(numTr,[denTr0])󰀁󰀁󰀁%部分分式展开式 Td=feedback(Gp,50)󰀁%计算Td(s)yd=step(Td,t)󰀁%扰动响应 figure(2);plot(t,yd);grid󰀁%扰动响应曲线运行以上程序后,我们可得响应曲线如图3和图4所示,对单位阶跃响应使用函数tstats,我们可求得下列性能指标:上升时间tr=0.83s;在tp=2󰀁1s时,最大超调量M0=39.5%;2%的调节时间ts2=9.0s;稳态误差为ess=0.4%󰀁扰动响应的稳态值为0.02
󰀁 2󰀁3󰀁动态系统的建模和仿真 MATLAB在工程技术方面另外一个令人振奋的应用是它的强大的系统建模和仿真能力.MAT-LAB提供的Simulink是一个用来对动态系统进行建模、仿真和分析的软件包,它支持线性和非线性系统,能够在连续时间域、离散时间域或者两者的混合时间域里进行建模,它同样支持具有多种采样速率的系统[4]󰀁下面以一个实例来说明它在这方面的应用󰀁 建立一个如图5所示的典型PID控制系统模 型,其中Gp(s)=(s3+7s2+24s+24)/(s4+10s3+35s2+50s+24),R(s)为单位阶跃输入,绘制出系 统的仿真曲线[1]
󰀁 利用Simulink解决以上问题变得很容易.
Simulink有一个空白的模型编辑窗口,允许用户输入自己的模型框图;也有现成的模块库,模块库里封装了一些常用的典型的模块,包含典型的输入信号、典型的环节以及输出信号显示模块等󰀁用户只需将自己所需要的模块从模块库拖入模型编辑窗口,然后用箭头将各模块连起来即可󰀁当然,模块库里没有的特殊模块需要用户自己定制,这项工作由MATLAB提供的S-函数来完成󰀁利用以上方法可建立本例的模型如图6所示,其中Step为单位阶跃输入信号;Gain0,Gain1和Gain2分别为比例,积分和微分增益;Integrator为积分环节;Derivative为微分环节;zeros(s)/poles(s)为被控对象GP(s)的传函的零极点表达式;Scope为模拟示波器,可以显示输出曲线的波形;outl为输出端口模块,可以将 输出值送到工作空间,这样用户可以利用绘图函数plot(tout,yout)绘制响应曲线的波形,利用绘图函数plot绘制的波形与Scope显示的波形是一样的,只是前者在使用时更灵活󰀁在启动仿真之前,用户可以任意修改Simulink默认的仿真控制参数,比如:仿真算法的选择、仿真范围的指定、仿真步长范围的指定以及仿真精度等󰀁另外,各模块的参数也是可以根据需要任意修改的,对如PID控制器来说,其参数的整定是一件很困难的事,在仿真的环境下,用户可以不断的调整其参数,一直到满足要求为止󰀁选择控制器参数KP=10,Ki=3,Kd=2,我们就可以启动仿真󰀁仿真结束后,我们可以得到系统的响应曲线如图7所示
三、 MATLAB在机械振动学中的应用
已知振动方程,求t=8秒时，位移、速度及加速度各位多少？并绘制振 动曲线和速度、加速度曲线。 编写MATLAB程序如下 t=8; %赋值 y=2*cos(3*t+pi/5); %确定y、t的涵数关系 v=-2*3*sin(3*t+pi/5); %确定v、t的涵数关系 a=-2*3*3*cos(3*t+pi/5) %确定a、t的涵数关系 敲回车键，得到结果如下 a = -15.7582 通过MATLAB绘制绘制振动曲线编写MATLAB程序如下 clc;%清除指令窗 t=0:0.001:10;%确定时间变量的变化范围及步长 y=2*cos(3*t+pi/5);%确定y、t的涵数关系 v=-2*3*sin(3*t+pi/5);%确定v、t的涵数关系，pi表示圆周率 a=-2*3*3*cos(3*t+pi/5);%确定a、t的涵数关系 plot(t,y,'-',t,v,'-.',t,a,':')%绘制图形 敲回车键，得到如图1

四、 MATLAB在液压系统仿真中的应用
所谓仿真,就是模仿真实事物,也就是用一个模型来模仿真实系统。基于这种定义,仿真分为物理仿真和数学仿真。数学仿真是通过数学模型来仿真动力学系统的幅频特性。 图2 工艺系统质量元件位移 (1)由图中曲线关系可以看出,工件和刀架振动发生在低频区,其频率f4=0.5kHz和f1=0.14kHz;刀杆和可换多面刀片发生在高频区,其频率f2=2.15kHz和f3=6.62kHz。 (2)用振动测量仪测量,在频率f=0～44.1时,记录加工数字信号,可以观测到所研究信号的时间 关系的不同形式的谱。 (3)利用数据总线接通计算机进行摸拟数字转换,由振动测量仪得到的信号转换成数字形式进行记录或加工,作为模数转换应用标准声图片,低频振动测量用霍尔式加速度传感器,包括振动测量仪的成份,而高频振动(f=1～7kHz)用DWY-3型加速器传感器。 (4)由上述用量工作时得到切断刀的幅频特性可知,切刀振动发生在高频和低频,刀杆属高频振动(f=2.05kHz),机械夹固可换多面刀片(f=6.53kHz),刀架和工件属工艺系统中低频(f=0.15～0.4kHz)振动。