现代设计方法及鲍威尔法课程设计

绪论◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆

一.有限元方法和优化方法的选择◆◆◆◆

二.大坝内应力分析◆◆◆◆◆◆◆◆◆◆◆◆ 2.1工程问题◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆ 2.2解析法求解◆◆◆◆◆◆◆◆◆◆◆◆◆◆ 2.3有限元法求解◆◆◆◆◆◆◆◆◆◆◆◆◆ 2.4结果分析与比较◆◆◆◆◆◆◆◆◆◆

三.鲍威尔法（powell）◆◆◆◆◆◆◆◆◆◆◆◆ 3.1鲍威尔法简介◆◆◆◆◆◆◆◆◆ 3.2鲍威尔法计算框图◆◆◆◆◆◆ 3.3问题与结果◆◆◆◆◆◆◆◆◆◆◆

附录◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆

参考文献◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆

绪论

有限元法是一种高能效，常用的计算方法。有限元法在早期是以变

分原理为基础发展起来的，所以它广泛的应用于拉普拉斯方程和泊松方程所描述的各种物理场中。自从1969年以来，，某些学者在流体力学中应用加权余数法中的迦辽金法或最小二乘法等同样获得了有限元方程，因为有限元法可以应用于任何微分方程所描述的各种物理场中，而不再要求这类物理场和泛函数的极值问题有所联系。基本思想：有解给定的泊松方程化为求解泛函数的极值问题。

它的原理：将连续的求解或离散的的一组单元的组合体，用在每个单元内假设的近似函数来分片的表示求解或尙待求的未知场函数来表达。从而使一个连续的无限自由度问题变成离散的有限自由度问题 经过多年的发展，有限元法已经成为现代结构分析的有效方法和主要手段，它的应用非常广泛，如：对拱坝、涡轮叶片、飞机和船体等复杂结构进行应力分析。由平衡问题扩展到稳定问题与动力问题。如：对结构在地震力与波浪力作用下的动力反应进行分析，由弹性力学问题扩展到弹塑性与粘弹性问题，如：土力学与岩石力学问题、疲劳力学与脆性断裂问题；由固体力学扩展到流体力学、渗流与固结理论，热传导与热应力问题以及建筑声学与噪声问题。由工程力学扩展到力学的其他领域，如：冰川与地质力学，血管与眼球力学等。由结构计算问题扩展到结构优化设计问题和可靠性问题。

有限单元由假定的应变方程式导出，有些单元可假设应变是常量，而另一些可采用更高阶的函数，利用给定单元的这些方程和实际几何体，则可以写出外力和节点位移之间的平衡方程。对于单元的每个节点来说，每个自由度就有一个方程，这些方程被十分便利的写成矩阵的形式以用于计算机的演算中。

近年来，机械优化设计的应用俞来愈广，但还面临许多问题需要解决。例如，机械产品设计中零部件的通用化、系列化和标准化，整机优化模型及方法的研究，机械设计中离散变量优化方法的研究，更为有效的优化设计方法的发掘等一系列问题，都需做比较大的努力才能适应机械工业发展的需要。

近年来，在计算机辅助设计中，应用优化方法后，使得在设计过程中即能不断选择设计参数并选出最优方案，又能加快设计速度，缩短设计周期。把优化设计与计算机辅助设计结合起来，使设计过程完全自动化，已成为设计方法的一个发展趋势。

一. 有限元方法和优化方法的选择

1.有限元问题：大坝利用ANSYS软件进行分析。

2.优化问题：利用鲍威尔法求解最优解。

二．大坝内应力分析

2.1工程问题

图示为一水坝示意图，水面高度为25m，坝体材料弹性模量为50GPa，泊松比为0.3，试对坝体进行应力分析。

P=5KN 2.2 解析法求解

此问题属于线性静力学问题，由于大坝的跨度远大于其他尺寸的方向，因此分析过程中按平面应变问题求解。应用有限元求解方法，，

粗略划分有限元进行求解，同后面软件计算结果进行比较。将均布载荷处的10*25划分成两个三角形单元，分别计算出两个单元的刚度，受力等等，再整体求解。经计算，单元1的计算数值是

m

I j bi=-25 bj =25 bm=0 ci=-10 cj=0 cm=10 应变矩阵

-5 0 5 0 0 0 [B]= 1/50 0 -2 0 0 0 2 -2 -5 0 5 2 0 单元载荷

[p]=125t[1/3 0 0 0 1/6 0] 刚度矩阵

经计算i、j,以及ij边的位移均为零，在i节点处的应力最大，为1.6*e6

2.3 有限元法求解

（一）建立工作文件名和工作标题

（1）选择Utility Menu/File/ Change Jobname命令，出现Jobname对话框。在【FILNAM】Enter nee jobname输入栏中输入工作文件名xiti,单击OK按钮，关闭该对话框。

（2）Utility Menu /File/ change Tile命令，出现Change Tile对话框，在输入栏中输入A DAM UNDER THE PREESURE OF WATER,单击OK按钮关闭该对话框。

（二）定义单元类型 （1）选择

main menue/preprocessor/element type

|add/edit/delete命令，出现element types对话框，单击add按钮，出现library of element types对话框，在library of element types列表中选择structural solid ,quad 8node 82，在element tyoe reference number输入栏中输入1，单击ok 按钮关闭该对话框。

（2）单击element type 对话框上的options按钮，出现plane82 element type option对话框，在element behavior K3下拉选项中选

择plane strain,单击OK按钮关闭该对话框。

（3）单击element type对话框上的close按钮，关闭该对话框。

（三）定义材料性能参数

（1）选择main menu|preprocessor|material props|material models命令，出现define material mode behavior对话框。

（2）在material model available一栏中依次双击structural 、linear、elastic、isotroptiac选项，出现linear isotroptiac properties for material numeber1对话框，在EX输入栏中输入5E10，在PRXY输入栏中输入0.3，单击OK按钮挂壁该对话框。

（3）在define material model behavior对话框中选择material|Exit命令，关闭该对话框。 （四）创建几何模型、划分网格

（1）选择main menu|preprocessor|modeling|creat|keypoint|in active CS命令，出现create keypoint in active coordinate systerm对话框。在NPT keypoint nember输入栏中输入1，在X,Y,Z location an active CS输入栏中分别输入0，0，0，如图示，单击apply按钮。

（2）参照上一步的操作步骤，依次创建一下关键点及其编号： 2（20,0,0）；3（10，30,0）；4（0,30,0）；5（0,25,0）

（3）单击create keypoint in ctive coordinate system对话框上的OK按钮关闭该对话框。

（4）选择 main menu | preprocessor | modeling | create | areas |arbitrary|through kps命令，出现create area thru拾取菜单，在输入栏中输入1,2,3,4,5，单击OK按钮关闭该菜单。

（5）选择main menu|preprocessor|meshing |size cntrls| manualsize|lines|picked lines命令，出现elemennt size on picked lines拾取菜单，在输入栏中输入1,3，单击OK按钮，出现element size on picked lines 对话框，在NDIV No.of element divisions输入栏中输入30，单击OK按钮关闭该对话框。如下图示。

提示：设置线段等份数。

（6）选择main menu|preprocessor|meshing| size cncrls |manualsize |lines |picked lines命令，出现elemennt size on picked lines拾取菜单，在输入栏中输入2, 单击OK按钮，出现element size on picked lines 对话框，在NDIV No.of element divisions输入栏中输入60，单击OK按钮关闭该对话框。

（7）选择main menu|preprocessor|meshing| size cncrls |manualsize |lines |picked lines命令，出现elemennt size on picked lines拾取菜单，在输入栏中输入4, 单击OK按钮，出现element size on picked lines 对话框，在NDIV No.of element divisions输入栏中输入10，单击OK按钮关闭该对话框。

（8）选择main menu|preprocessor|meshing| size cncrls |manualsize |lines |picked lines命令，出现elemennt size on picked lines拾取菜单，在输入栏中输入5, 单击OK按钮，出现element size on picked lines 对话框，在NDIV No.of element divisions输入栏中输入50，单击OK按钮关闭该对话框。

（9）选择main menu|preprocessor|meshing| mesh| areas| mapped |by corners命令。出现map mesh area by拾取菜单，在输入栏中输入1，单击OK按钮，再次出现map mesh areas by拾取菜单，，在输入栏中输入，4,3,2,1，单击OK 按钮关闭菜单。

（10）选择utility menu|select |everything命令。

（11）选择utility menu| plot |elements 命令，ansys显示窗口显示网格划分结果。，如图示，单击OK按钮关闭该对话框。

（12）选择utility menu |file |save as命令，出现save database对话框，在save database to输入栏中输入xiti.db,保存操作过程，单击OK按钮关闭该对话框。 （五）加载求解

（1）选择main menu| solution |analysis type|new analysis命令，出现new analysis 对话框，选择分析类型为static，单击OK按钮关闭该对话框。

（2）选择选择

main menu| solution |define loads

|apply|functions|define/edit命令，出现function editor对话框，参照图进行设置，选择对话框上的FILE|save命令，出现“另存为”对话框，将所做函数取名为FUNC保存在ANSYS工作目录下，单击：“保存”按钮关闭该对话框，在function editor对话框上选择file|close按

钮关闭该对话框。

函数编辑器对话框

（3）选择main menu| solution |define loans |apply |functions

|read file命令，出现“打开”对话框，选择FUNC.func，单击“打开”按钮，出现 function loader对话框，在table parameter name输入栏中输入PRES1，单击OK按钮 关闭该对话框。

（4）选择main menu |solution |define loads |apply | structural| on

lines |命令，出现apply PRES on lines 拾取菜单，在输入栏中输

入5，单击OK按钮，出现apply PRES on lines对话框，在【SFL】 apply PRES on lines as a下拉选框中选择existing table，单击OK按钮，出现apply PRES on lines对话框，在列表框中选择PRES1，如图所示

（5）选择main menu| solution |define loads| apply |structural |

displacement |on lines命令，出现apply U,ROT on L拾取菜单，

在输入栏中输入1，单击OK按钮，出现applyU,ROT on lines对话框，参照下图设置，单击OK 按钮关闭该对话框。

（6）选择utility menu |file | save as 命令，出现save database to

输入栏中输入xiti.db,保存操作过程，单击OK按钮关闭该对话框。

（7）选择main menu|solution |sove|current LS命令，单击sove

current load step对话框上的OK按钮，ansys开始求解计算。求解结束后，ansys显示窗口出现NOTE提示框，单击close按钮关闭该对话框。

（8）选择utility menu |file| save as命令，出现save database to 输

入栏中输入xiti.db,保存求解结果，单击OK按钮关闭该对话框。

（六）查看求解结果

（1）选择main menu |general postproc| contour plot |nodal solu 命

令，出现contour nodal solution data对话框。选择nodal solution} displacement vector sum,单击OK按钮，ANSYS显示窗口显示合位移等值线图。

（2）选择main menu |general postproc| contour plot |nodal solu

命令，出现contour nodal solution data对话框，选择nodal solution |stess| von mises stress,单击OK按钮，ANSYS显示窗口显示等效应力场等值线图

。 合位移等值线图

等效应力场等效线图

（3）选择utility menu|list |nodes,出现各个节点的信息。

（4）选择utility menu| file|exit命令，出现EXITfrom ansys对话框，

选择Quit-No Save!,单击OK按钮，关闭ANSYS. 2.4 结果分析与比较

由于计算量很大，在此只计算比较了节点I处的位移与应力，通过

理论计算的结果与ansys软件计算的结果相对比，位移结果相同。应力结果有所偏差，是在计算过程中的有效数字取舍导致，可认为结果一致。

三 鲍威尔法（ Powell法）

3.1鲍威尔法简介

• 两次平行搜索产生一个共轭方向，Powell法也是一种共轭方向法，能在有限步长内极小化一个二次函数，是直接搜索方法中使用效果最佳的一种方法。

• 对于维数n<20的目标函数求最优化问题，此法可获得满意效果。 Ⅰ、鲍威尔法基本原理、迭代格式

• 原始的Powell法是沿着逐步产生的共轭方向进行一维搜索的。 • 现以二维二次目标函数为例来说明。 如下图所示，选定初始点X0(1)，初始方向： S1(1)=e1=[1,0]T S2(1)=e2=[0,1]T

• 由图可知点X0(2) 、X2(2)是先后两次沿S(1)方向一维搜索的极小点。 • 由共轭性质知：连接X0(2) ，X2(2)构成的矢量S(2) 与S(1)对H共轭。 • 从理论上讲，二维二次正定函数经过这组共轭方向的一维搜索，迭代点已达到函数的极小点X* 。

• 将此结构推广至n维二次正定函数，即依次沿n个（S(1) ，S(2)，…，S(n)）共轭方向一维搜索就能达到极小点。 Ⅱ、鲍威尔法缺陷

• 当某一循环方向组中的矢量系出现线性相关的情况（退化、病态）时，搜索过程在降维的空间进行，致使计算不能收敛而失败。

• 为了避免鲍威尔法缺陷，提出了修正算法。

• 和原始Powell法的主要区别在于：在构成第k+1次循环方向组时，不用淘汰前一循环中的第一个方向S1(k)的办法，而是计算函数值并根据是否满足条件计算： f1=f(Xk(0))

f2=f(Xk(n)) f3=f(Xk(n+2))

• 找出前一轮迭代法中函数值下降最多的方向m及下降量△m，即：

△m=max｛[f(Xk(i))－f(Xk(i+1))]（i=0,1,…,n-1）｝ = f(Xk(m-1))－f(Xk(m)) • 可以证明：若 f3 22 (f1 -2f2+f3)(f1-f2- △m)< 0.5△m(f1-f3)

同时成立

• 表明方向Sk(n)与原方向组成线性无关，可以用来替换对象△m所对应的方向Sk(m)。否则仍用原方向组进行第k+1轮搜索。

3.2鲍威尔法计算框图

3.3问题与结果

附录（源程序）

#include \"stdio.h\" #include \"stdlib.h\" #include \"math.h\" double objf(double x[]) {double ff;

ff=x[0]*x[0]+2*x[1]*x[1]-4*x[0]-2*x[0]*x[1]; return(ff); }

void jtf(double x0[],double h0,double s[],int n,double a[],double b[]) {int i;

double *x[3],h,f1,f2,f3; for(i=0;i<3;i++)

x[i]=(double *)malloc(n*sizeof(double)); h=h0;

for(i=0;i*(x[1]+i)=*(x[0]+i)+h*s[i]; f2=objf(x[1]); if(f2>=f1) {h=-h0;

for(i=0;ifor(i=0;ifor(i=0;i*(x[2]+i)=*(x[1]+i)+h*s[i]; f3=objf(x[2]); if(f2{ for(i=0;i{*(x[0]+i)=*(x[1]+i); *(x[1]+i)=*(x[2]+i); } f1=f2; f2=f3; } }

if(h<0)

for(i=0;ifor(i=0;ifor(i=0;i<3;i++) free(x[i]); }

double gold(double a[],double b[],double eps,int n,double xx[]) {int i;

double f1,f2,*x[2],ff,q,w; for(i=0;i<2;i++)

x[i]=(double *)malloc(n*sizeof(double)); for(i=0;i{*(x[0]+i)=a[i]+0.618*(b[i]-a[i]); *(x[1]+i)=a[i]+0.382*(b[i]-a[i]); }

f1=objf(x[0]); f2=objf(x[1]); do

{if(f1>f2)

{for(i=0;i*(x[0]+i)=*(x[1]+i); } f1=f2;

for(i=0;i*(x[1]+i)=a[i]+0.382*(b[i]-a[i]); f2=objf(x[1]);

} else

{ for(i=0;i*(x[1]+i)=*(x[0]+i);} f2=f1;

for(i=0;i*(x[0]+i)=a[i]+0.618*(b[i]-a[i]); f1=objf(x[0]); } q=0;

for(i=0;iq=q+(b[i]-a[i])*(b[i]-a[i]); w=sqrt(q); }while(w>eps); for(i=0;ixx[i]=0.5*(a[i]+b[i]); ff=objf(xx); for(i=0;i<2;i++) free(x[i]); return(ff); }

double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[]) {double *a,*b,ff;

a=(double *)malloc(n*sizeof(double)); b=(double *)malloc(n*sizeof(double)); jtf(x0,h0,s,n,a,b); ff=gold(a,b,epsg,n,x); free(a); free(b); return (ff); }

double powell(double p[],double h0,double eps,double epsg,int n,double x[]) {int i,j,m;

double *xx[4],*ss,*s;

double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;

ss=(double *)malloc(n*(n+1)*sizeof(double)); s=(double *)malloc(n*sizeof(double)); for(i=0;i*(ss+i*(n+1)+i)=1; }

for(i=0;i<4;i++)

xx[i]=(double *)malloc(n*sizeof(double)); for(i=0;i{for(i=0;i{*(xx[1]+i)=*(xx[0]+i); x[i]=*(xx[1]+i); }

f0=f1=objf(x); dlt=-1;

for(j=0;j*(s+i)=*(ss+i*(n+1)+j); }

f=oneoptim(xx[0],s,h0,epsg,n,x); df=f0-f; if(df>dlt) {dlt=df; m=j; } }

sdx=0;

for(i=0;isdx=sdx+fabs(x[i]-(*(xx[1]+i))); if(sdxfor(i=0;i<4;i++) free(xx[i]); return(f); }

for(i=0;ifor(i=0;i{*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i))); x[i]=*(xx[3]+i); }

fx=objf(x); f3=fx;

q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt); d=0.5*dlt*(f1-f3)*(f1-f3); if((f3for(i=0;i*(xx[0]+i)=*(xx[2]+i); else

for(i=0;i*(xx[0]+i)=*(xx[3]+i); } else

{for(i=0;i{*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i)); *(s+i)=*(ss+(i+1)*(n+1)); }

f=oneoptim(xx[0],s,h0,epsg,n,x); for(i=0;i*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j); } }

void main()

{double p[]={1,2}; double ff,x[2];

ff=powell(p,0.3,0.001,0.0001,2,x);

printf(\"x[0]=%f,x[1]=%f,ff=%f\\n\getchar(); }

运行结果：

参考文献

1.《现代机械设计方法》 化学工业出版社 倪洪启、谷耀新主编 2.《材料力学》 高等教育出版社 刘鸿文主编

3.《有限元分析理论与工程应用》 电子工业出版社 张朝晖、李树奎

主编

4.《C语言程序设计》 清华大学出版社 谭浩强主编 5.《ANSYS>工程应用实例解析》 机械工业出版社 龚曙光主编