计算机辅助设计制造实验指导书 精品.docx
《计算机辅助设计制造实验指导书 精品.docx》由会员分享,可在线阅读,更多相关《计算机辅助设计制造实验指导书 精品.docx(64页珍藏版)》请在冰豆网上搜索。
计算机辅助设计制造实验指导书精品
Hewlett-Packard
计算机辅助设计与制造
实验指导书
姜志宏
2012-6
该指导书为工业工程专业学习设计,旨在增强学生工业背景、计算机应用能力。
配套教材:
刘极峰.计算机辅助设计与制造.高等教育出版社
实验一AutoCAD绘制三维实体的方法和技巧
一、实验目的
1.熟悉创建实体的命令及功能;
2.掌握实体的编辑及三维操作;
3.绘制如图7-1所示的三维实体。
图7-1组合体
二、实验设备
1.586以上主机一台
2.显示器
3.鼠标器
4.键盘
5.打印机或绘图仪
三、实验内容:
创建图7-1所示立体的三维实体模型。
绘图步骤:
1.分析图形:
该立体可分为由左右对称的两个耳板和中间开有矩形槽的空心圆柱组成。
2.绘制平面图形,并把俯视图上的左右两个耳板创建成面域。
绘图(下拉式菜单)→边界→边界创建→面域→拾取点→拾取内部点↙
3.单击(西南等轴测)图标,视图显示方式为轴测图。
单击(建模工具栏)→拾取俯视图中左右耳板、耳板上的圆→指定高度6↙
单击(建模工具栏)→拾取俯视图中间的两个圆→指定高度20↙
拉伸后如图7-2所示。
4.在俯视图上绘制10×40的矩形,并定义成面域。
单击(建模工具栏)→拾取矩形→指定高度6↙
以长方体上表面中心为基准,将其移动到圆柱顶面;使长方体的上表面中心与圆柱顶面
图7-2在俯视图上拉伸出立体
的圆心对齐,如图7-3所示。
5.单击(差集)→拾取大圆柱、左右两个耳板↙→小圆柱、长方体、耳板上的两个小圆柱↙;
6.单击(并集)→窗选所有实体↙→完成实体绘制如图7-4所示。
7.单击(概念视觉样式)↙,效果图见图7-5所示。
图7-3移动长方体到圆柱上端面
图7-4布尔运算结果
图7-5着色后效果
实验三梁平面框架结构的有限元分析
实验目的:
了解MATLAB有限元分析功能;
了解有限元分析的过程;
实验内容:
如图3-19所示的框架结构,其顶端受均布力作用,结构中各个截面的参数都为:
,
,
。
试基于MATLAB平台求解该结构的节点位移以及支反力。
图3-19框架结构受一均布力作用
步骤:
对该问题进行有限元分析的过程如下。
(1)结构的离散化与编号
将该结构离散为3个单元,节点位移及单元编号如图3-20所示,有关节点和单元的信息见表3-5。
(a)节点位移及单元编号
(b)等效在节点上的外力
图3-20单元划分、节点位移及节点上的外载
(2)各个单元的描述
首先在MATLAB环境下,输入弹性模量E、横截面积A、惯性矩I和长度L,然后针对单元1,单元2和单元3,分别二次调用函数Beam2D2Node_ElementStiffness,就可以得到单元的刚度矩阵k1(6×6)和k2(6×6),且单元2和单元3的刚度矩阵相同。
>>E=3E11;
>>I=6.5E-7;
>>A=6.8E-4;
>>L1=1.44;
>>L2=0.96;
>>k1=Beam2D2Node_Stiffness(E,I,A,L1);
>>k2=Beam2D2Node_Stiffness(E,I,A,L2);
(3)建立整体刚度方程
将单元2和单元3的刚度矩阵转换成整体坐标下的形式。
由于该结构共有4个节点,则总共的自由度数为12,因此,结构总的刚度矩阵为KK(12×12),对KK清零,然后两次调用函数Beam2D2Node_Assemble进行刚度矩阵的组装。
>>T=[0,1,0,0,0,0;-1,0,0,0,0,0;0,0,1,0,0,0;0,0,0,0,1,0;0,0,0,-1,0,0;0,0,0,0,0,1];
>>k3=T'*k2*T;
>>KK=zeros(12,12);
>>KK=Beam2D2Node_Assemble(KK,k1,1,2);
>>KK=Beam2D2Node_Assemble(KK,k3,3,1);
>>KK=Beam2D2Node_Assemble(KK,k3,4,2)
KK=1.0e+008*
1.443100.0127-1.416700-0.026400.0127000
02.13280.00560-0.00780.00560-2.12500000
0.01270.00560.01350-0.00560.0027-0.012700.0041000
-1.4167001.443100.0127000-0.026400.0127
0-0.0078-0.005602.1328-0.00560000-2.12500
00.00560.00270.0127-0.00560.0135000-0.012700.0041
-0.02640-0.01270000.02640-0.0127000
0-2.1250000002.12500000
0.012700.0041000-0.012700.0081000
000-0.02640-0.01270000.02640-0.0127
0000-2.1250000002.12500
0000.012700.0041000-0.012700.0081
(4)边界条件的处理及刚度方程求解
该问题的位移边界条件为
。
因此,将针对节点1和节点2的位移进行求解,节点1和节点2的位移将对应KK矩阵中的前6行和前6列,则需从KK(12×12)中提出,置给k,然后生成对应的载荷列阵p,再采用高斯消去法进行求解。
注意:
MATLAB中的反斜线符号“\”就是采用高斯消去法。
>>k=KK(1:
6,1:
6);
>>p=[3000;-3000;-720;0;-3000;720];
>>u=k\p
u=0.0009-0.0000-0.00140.0009-0.0000-0.0000[将列排成了行]
(5)支反力的计算
在得到整个结构的节点位移后,由原整体刚度方程就可以计算出对应的支反力;先将上面得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系),可以得到整体的位移列阵U(12×1),再代回原整体刚度方程,计算出所有的节点力P(12×1),按式(3-180)的对应关系就可以找到对应的支反力。
>>U=[u;0;0;0;0;0;0]
U=0.0009-0.0000-0.00140.0009-0.0000-0.0000[将列排成了行]
000000[将列排成了行]
>>P=KK*U
P=1.0e+003*
3.0000-3.0000-0.72000.0000-3.00000.7200[将列排成了行]
-0.66582.20120.6014-2.33423.79881.1283[将列排成了行]
附件:
创建三个函数:
Beam2D2Node_Stiffness(E,I,A,L)
以上函数计算单元的刚度矩阵,输入弹性模量E,横截面积A,惯性矩I,长度L,输出单元刚度矩阵k(6
×6)。
Beam2D2Node_Assemble(KK,k,i,j)
以上函数进行单元刚度矩阵的组装,输入单元刚度矩阵k,单元的节点编号i、j,输出整体刚度矩阵KK。
Beam2D2Node_Forces(k,u)
以上函数计算单元的应力,输入单元刚度矩阵k,节点位移u,输出单元节点力forces。
每个函数的MATLAB程序如下。
%%%%%%%%%%Beam2D2Node%%begin%%%%%%%%%%%%
functionk=Beam2D2Node_Stiffness(E,I,A,L)
%以上函数计算单元的刚度矩阵
%输入弹性模量E,横截面积A,惯性矩I,长度L
%输出单元刚度矩阵k(6×6)
%-----------------------------------------
k=[E*A/L,0,0,-E*A/L,0,0;0,12*E*I/(L^3),6*E*I/(L^2),0,-12*E*I/(L^3),6*E*I/(L^2);0,6*E*I/(L^2),4*E*I/L,0,-6*E*I/
(L^2),2*E*I/L;-E*A/L,0,0,E*A/L,0,0;0,-12*E*I/(L^3),-6*E*I/(L^2),0,12*E*I/(L^3),-6*E*I/(L^2);0,6*E*I/(L^2),2*
E*I/L,0,-6*E*I/(L^2),4*E*I/L]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functiony=Beam2D2Node_Assemble(KK,k,i,j)
%以上函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k,单元的节点编号i、j
%输出整体刚度矩阵KK
%-----------------------------------------
DOF
(1)=3*i-2;
DOF
(2)=3*i-1;
DOF(3)=3*i;
DOF(4)=3*j-2;
DOF(5)=3*j-1;
DOF(6)=3*j;
forn1=1:
6
forn2=1:
6
KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
end
end
y=KK;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionforces=Beam2D2Node_Forces(k,u)
%以上函数计算单元的应力,输入单元刚度矩阵k,节点位移u,
%输出单元节点力forces
%-----------------------------------------
forces=k*u;
%%%%%%%%%%Beam2D2Node%%end%%%%%%%%%%%%
实验四基于MATLAB语言的机械优化设计
机械优化设计是在现代机械设计理论发展基础上产生的一种新的设计方法,在连杆、凸轮、齿轮、涡轮、轴承、传动轴、机床等产品的机械设计的应用中取得了良好的效果。
机械优化设计作用是在进行某种机械产品设计时,可根据规定的约束条件,优选设计参数,使某项或几项设计指标获得最优值。
实际中对一个优化问题的处理,主要有两个步骤:
首先要把工程设计问题转化成数学模型,即用数学表达式描述工程实际问题;然后按照数学模型的特点选择优化方法及其计算程序,作必要的简化和加工,用计算机求得最优设计方案。
目前,已经有很多成熟的优化程序可供选择,但由于机械优化设计的复杂性和程序算法的不稳定性,实际中往往会因为方法选择不当或初始值设置不合理而导致优化失败或只能得到局部优化结果。
MATLAB优化工具箱
美国MathWorks公司推出的MATLAB(MatrixLaboratory)是一套功能强大的工程计算软件,集数值计算、符号运算、可视化建模、仿真和图形处理等多种功能于一体,被广泛应用于机械设计、自动控制和数理统计等工程领域。
工程人员通过使用MATLAB提供的工具箱,可以高效的求解复杂的工程问题。
MATLAB优化工具箱提供了对各种问题的一个完整的解决方案,广泛应用于线性规划、二次规划、非线性优化、最小二乘问题、非线性方程求解、多目标决策、最小最大问题以及半无限问题。
其函数表达简洁,优化算法选择任意,参数设置自由,相比于其它很多成熟的优化程序具有明显的优越性。
大多数机械优化设计问题属于非线性多变量约束最小化问题,约束条件有等式和不等式约束,目标函数和约束函数中有一个或多个为非线性函数。
MATLAB通过调用fmincon函数来实现这一问题的优化。
优化问题可作下列描述:
minf(X)X∈Rn
约束条件:
AX≤b(线性不等式约束)
AeqX=beq(线性等式约束)
C(X)≤0(非线性不等式约束)
Ceq(X)=0(非线性等式约束)
Lb≤X≤Ub(边界约束)
fmincon函数相应的调用格式为[2]:
[x,f,ex2itflag,output,lambda,grad,hessian]=fmincon(‘fun’,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2,⋯)
其中x为最优设计点,f为目标函数在最优点x的值,exitflag是返回算法的终止标志,exitflag≤0,目标函数不收敛;exitflag>0,目标函数收敛。
output负责返回优化结果信息,主要有迭代次数、最后一次迭代计算次数、步长、算法、一阶导数等,lambda返回Lagrange乘子值,grad返回目标函数在最优点x的梯度值,hessian返回目标函数在最优点x的hessian矩阵。
优化实例
欲设计如图1所示某型号空心传动轴,D和d分别为空心轴的外径和内径,轴长L=4m。
轴的材料密度ρ=7.8×103kg/m3,剪切弹性模量G=80GPa,许用剪应力[τ]=40MPa,单位长度许用扭转角[φ]=1°/m,轴所传递的功率为P=5.5kW,转速n=200r/min。
要求在满足使用条件和结构尺寸限制的前提下使其质量最小。
3.1设计变量和目标函数
该传动轴的力学模型是一个受扭转的圆柱筒轴。
其外径D和内径d是决定圆轴的重要独立参数,故将其作为设计变量。
写成向量形式为:
X=[x1 x2]T=[D d]T
取质量最小为优化目标。
则目标函数空心圆轴的质量可按下式计算:
M=π/4*ρL(D2-d2)
可见这是一个合理选择D和d而使质量M最小的优化问题。
3.2约束条件
应满足的使用条件和结构尺寸限制是:
(1)扭转强度
根据扭转强度,要求扭转剪应力须满足:
τmax=TWt≤[τ]
式中,T是圆轴所受扭矩,T=9549P/n,Wt是抗扭截面模量,Wt=π(D4-d4)/16D。
(2)扭转刚度
为了确保传动轴正常工作,除满足强度条件外,还要限制轴的变形。
限制变形条件即为刚度条件。
通常要求单位长度的最大扭转角不超过规定的许用值,即:
φ=TGIp≤[φ]
式中,φ是单位长度扭转角,G是剪切模量,Ip是
极惯性矩。
(3)结构尺寸
由结构尺寸要求决定的约束条件为:
d≥0
D≥d
3.3优化模型
将所有函数表达式规范化并代入已知数据,可得传动轴优化设计的数学模型为:
minf(X)=24.504×10-6×(x12-x22)
满足约束条件:
g1(X)=33435×x1x14-x24-1≤0
g2(X)=1.9157×106x14-x24-1≤0
g3(X)=-x2≤0
g4(X)=-x2-x1≤0
3.4MATLAB优化程序实现
此例为非线性约束优化问题。
利用文件编辑器为目标函数编写M文件(fun.m):
functionf=fun(x)
f=24.504e-63(x
(1)^2-x
(2)^2);
编写约束函数的M文件(nonlcon.m):
function[c,ceq]=nonlcon(x)
c=[33435e43x
(1)/(x
(1)^4-x
(2)^4)-1;1.9157e6/(x
(1)^4-x
(2)^4)-1];
ceq=[];
在命令窗口编写主程序(调用函数):
x0=[2010];
A=[-11];
b=0;
lb=[00];
[x,f,exitflag,output]=fmincon(‘fun’,x0,A,b,[],[],lb,[],‘nonlcon’)
得到运行结果如下:
x=37.2301 8.6223
f=0.0321
exitflag=1
output=
iterations:
7
funcCount:
31
stepsize:
1
algorithm:
‘medium-scale:
SQP,
Quasi-Newton,line-search’
firstorderopt:
3.9983e-004
cgiterations:
[]
实验五用Matlab进行数据拟合
1.多项式曲线拟合:
polyfit
p=polyfit(x,y,m)
其中,x,y为已知数据点向量,分别表示横,纵坐标,m为拟合多项式的次数,结果返回m次拟合多项式系数,从高次到低次存放在向量p中.
y0=polyval(p,x0)
可求得多项式在x0处的值y0.
例1已知观测数据点如表所示
分别用3次和6次多项式曲线拟合这些数据点.
编写Matlab程序如下:
x=0:
0.1:
1
y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2]
plot(x,y,'k.','markersize',25)
axis([01.3-216])
p3=polyfit(x,y,3)
p6=polyfit(x,y,6)
x=0:
0.1:
1
y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2]
plot(x,y,'k.','markersize',25)
axis([01.3-216])
p3=polyfit(x,y,3)
p6=polyfit(x,y,6)
t=0:
0.1:
1.2
s=polyval(p3,t)
s1=polyval(p6,t)
holdon
plot(t,s,'r-','linewidth',2)
plot(t,s,'b--','linewidth',2)
grid
注意观察输出,并对各行添加注释。
例2用切削机床进行金属品加工时,为了适当地调整机床,需要测定刀具的磨损速度.在一定的时间测量刀具的厚度,得数据如表所示:
解:
描出散点图,在命令窗口输入:
t=[0:
1:
16]
y=[30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0]
plot(t,y,'*')
描出散点图,在命令窗口输入:
t=[0:
1:
16]
y=[30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0]
plot(t,y,'*')
a=polyfit(t,y,1)
holdon
y1=-0.3012*t+29.3804
plot(t,y1),holdoff
例3一个15.4cm×30.48cm的混凝土柱在加压实验中的应力-应变关系测试点的数据如表所示
1.55
2.47
2.93
3.03
2.89
500e-6
1000e-6
1500e-6
2000e-6
2375e-6
3.103e3
2.465e3
1.953e3
1.517e3
1.219e3
已知应力-应变关系可以用一条指数曲线来描述,即假设
式中,
表示应力,单位是N/m2;
表示应变.
解选取指数函数作拟合时,在拟合前需作变量代换
化为k1,k2的线性函数
在命令窗口输入:
x=[500*1.0e-61000*1.0e-61500*1.0e-62000*1.0e-62375*1.0e-6]
y=[3.103*1.0e+32.465*1.0e+31.953*1.0e+31.517*1.0e+31.219*1.0e+3]
z=log(y)
a=polyfit(x,z,1)
k1=exp(8.3009)
w=[1.552.472.933.032.89]
plot(x,w,'*')
y1=exp(8.3009)*x.*exp(-494.5209*x)
plot(x,w,'*',x,y1,'r-')
求出a数组中的值即为k1,k2
在实际应用中常见的拟合曲线有:
直线、多项式、双曲线、指数曲线。
2.非线性曲线拟合:
nlinfit.
x=lsqcurvefit(fun,x0,xdata,ydata)
[x,resnorm]=lsqcurvefit(fun,x0,xdata,ydata)
功能:
根据给定的数据xdata,ydata(对应点的横,纵坐标),按函数文件fun给定的函数,以x0为初值作最小二乘拟合,返回函数fun中的系数向量x和残差的平方和resnorm.
例4已知观测数据点如表所示
求三个参数a,b,c的值,使得曲线f(x)=aex+bx2+cx3与已知数据点在最小二乘意义上充分接近.
编写下面的程序调用拟合函数
x=0:
0.1:
1;
y=[3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17];
x0=[0,0,0];
[beta,r,J]=nlinfit(x',y','nihehanshu',x0);
求三个参数a,b,c的值,使得曲线f(x)=aex+bx2+cx3与已知数据点在最小二乘意义上充分接近.
说明:
最小二乘意义上的最佳拟合函数为f(x)=3ex+4.03x2+0.94x3.
此时的残差是:
0.0912.
拟合函数为:
f(x)=3ex+4.03x2+0.94x3.
实验作业:
1.已知观测数据点如表所示
求用三次多项式进行拟合的曲线方程.
2.已知观测数据点如表所示
求a,b,c的值,使得曲线f(x)=aex+bsinx+clnx与已知数据点在最小二乘意义上充分接近.
实验六用MATLAB实现拉格朗日插值、分段线性插值
一、实验目的:
1)学会使用MATLAB软件;
2)会使用MATLAB软件进行拉格朗日插值算法和分段线性差值算法;
二、实验内容:
1用MATLAB实现y=1./(x.^2+1);(-1<=x<=1)的拉格朗日插值、分段线性
2.选择以下函