十杆桁架结构优化设计

1、题目：十杆桁架结构优化设计日期：2013.09.16 目录1设计题目12设计过程22.1一、运用abaqus求解各杆轴力应力22.1.1abaqus计算流程22.1.2结果32.2二、利用材料力学知识求解42.2.1基本思路42.2.2解题过程42.2.3结果52.3三、编写有限元程序求解62.3.1程序基本步骤62.3.2vs2012 中重要的程序段62.3.3程序输出文件92.3.4材料力学、有限元程序、abaqus结果比较102.4四、装配应力计算112.4.1处理技巧112.4.2abaqus处理技巧112.4.3不加外力（p1,p2,p3） 时 材力，ansys与abaqus结果12

2、2.4.4不加外力（p1,p2,p3） 时 材力，ansys与abaqus误差分析122.4.5加外力（p1,p2,p3）时 ansys与abaqus结果122.5五、优化设计142.5.1设计中变量的概念142.5.2优化步骤运用vs2012编写复合形法进行约束优化。142.5.3vs2012优化程序162.5.4优化结果192.5.5结果说明193设计感想204备注204.1参考书目204.2说明201 设计题目十字桁架结构优化设计现有十字桁架结构见图1，材料泊松比为0.3，e=2.1e11，密度为7.8103kg/m3, 许用应力为160mpa，p1=600k n ，p2=900k n

3、，p3=600k n，杆1-6面积为a1=0.03m2,杆7-10面积为a2=0.02m。1、利用计算各杆的应力；2、利用材料力学的知识求解，并与1计算出的结果做比较；3、编写有限元程序求解，与1和2计算结果进行比较；4、若杆5制作时短了0.001m，试求各杆的应力；5、若令2节点的位移小于0.005m，a1、a2为0.0050.05m2，试对结构进行优化，使其重量最小。 （同材料力学优化结果比较）。图 1十杆桁架2 设计过程2.1 一、运用abaqus求解各杆轴力应力利用abaqus求解，十字桁架结构可用2dtruss单元模拟。单元参数为：弹性模量，1-6杆截面面积，7-10杆截面面积。加载

4、求解输出各杆应力，结点位移。2.1.1 abaqus计算流程part:创建part trussmain,part45l,part45r,part trussmain,包涵除8,10杆外的所有杆，part45l包涵8杆，part45r包涵10杆。property：create material:elastic: 弹性模量,泊松比0.3 create section：beamtruss：section a1,截面面积30000。section a2,截面面积20000。并给各杆赋材料属性。assembly：组装part trussmain,part45l,part45r。step:创建一个分析步，

5、step1。interaction： 用tie把part trussmain,part45l,part45r,绑定。load： create load：5,6点加铰接约束，固定x,y方向位移。 create boundary condition：2,4点加相应力。mesh：划分网格，一个杆为一个单元。element tape，选trussjob：创建一个job，write input,data check,submit,通过result来查看应力云图。2.1.2 结果 图 2 abaqus各杆应力云图2.2 二、利用材料力学知识求解2.2.1 基本思路显然题目中的十字桁架结构是两次静不定问题。

6、对于一次静不定问题，材料力学给出了两类解法：去掉约束加力，找位移协调关系解题；力法正则方程求解。对于多次静不定，特别是上述桁架问题，找出其协调关系基本上是不可能的，而力法正则方程更适合于解这种结构。如图3所示，去掉多余约束，建立力法正则方程：图 3 去多余约束2.2.2 解题过程分别求出外力作用下各杆内力和单位力作用下的各杆内力，为计算方便，将其结果列入下表1中。应用莫尔积分定理有：表 1外力作用下各杆内力和单位力作用下的各杆内力杆号lp1p2p3fi1fi21a001-12a000013a1-14a0015a000106a000017008009000100000杆号轴力f(n)应力s(mp

7、a)1158056952.672347433.111.583-1119431-37.314474331.585128002.14.276347433.111.5871158850.457.948-962469.9-48.129781447.639.0710-491344.6-24.572.2.3 结果表 2 材料力学各杆应力结果2.3 三、编写有限元程序求解2.3.1 程序基本步骤 计算单元刚度矩阵单元坐标系下刚度矩阵：yxeye35x图 4 单元坐标系和结构坐标系结构坐标下刚度矩阵：组装总的刚度矩阵边界条件处理（固定约束，直接去掉约束对应的行和列）计算位移向量计算单元应力2.3.2 vs20

8、12 中重要的程序段 计算结点位移 private sub button1_click(sender as object, e as eventargs) handles button1.click 桁架结点位移计算 形成总刚度矩阵 dim tk(12, 12) as double 总体刚度矩阵 tk = matrix.stiffsofalltk() dim tkh(11, 11) as double 去除0行0列 for i = 1 to 12 for j = 1 to 12 tkh(i - 1, j - 1) = tk(i, j) next j next i 输入结点载荷p(i) dim

9、p(12) as double p = data.nodeloaddata() dim ph(11) as double for i = 1 to 12 ph(i - 1) = p(i) next i 边界条件处理 for i = 8 to 11 for j = 1 to 11 tkh(i, j) = 0.0 next j next i for j = 8 to 11 tkh(j, j) = 1.0 next j for j = 8 to 11 ph(j) = 0.0 next j 计算结点位移 dim z(11) as double 结点位移 dim tkht(11, 11) as doub

10、le 去除0行0列 tkht = matrix.inversionofmatrix(tkh) z = matrix.matrixmultiplevector(tkht, ph) 输出结点位移 io.output(z) end sub private sub button2_click(sender as object, e as eventargs) handles button2.click 桁架单元内力 计算总体坐标架单元新节点位移xnew(6，2) dim xnew(6, 2) as double dim znew(6, 2) as double dim x(6, 2) as doubl

11、e x = data.positiondata() dim z(12) as double 结点位移 z = matrix.displacement() for i = 1 to 6 znew(i, 1) = z(2 * i - 2) znew(i, 2) = z(2 * i - 1) next i xnew = matrix.add(x, znew) 计算变形后杆长 dim ddeltax(10) as double dim newddeltax(10) as double dim d(10) as double dim a(10) as double dim e1 as integer a

12、 = data.areadata e1 = data.edata dim n(10) as double 单元内力 for i = 1 to 10 x = data.positiondata() 单元结点编号 dim nx(2, 10) as double nx = data.nodedata() newddeltax(i) = math.sqrt(xnew(nx(1, i), 1) - xnew(nx(2, i), 1) 2 + (xnew(nx(1, i), 2) - xnew(nx(2, i), 2) 2) ddeltax(i) = matrix.elementlongger(i) d(

13、i) = newddeltax(i) - ddeltax(i) 计算单元内力 n(i) = d(i) * e1 * a(i) / ddeltax(i) next i 输出单元内力 n io.output(n) end sub 计算单元应力过程 private sub button3_click(sender as object, e as eventargs) handles button3.click dim s(10) as double dim f(10) as double dim a(10) as double a = data.areadata f = matrix.force f

14、or i = 1 to 10 s(i) = f(i) / a(i) next i 输出单元应力 s io.output(s) end sub2.3.3 程序输出文件图 5 有限元位移结果 图 6有限元应力结果图 7有限元轴力结果2.3.4 材料力学、有限元程序、abaqus结果比较表 3 材力、有限元、abaqus计算结果比较表 4 材力、有限元、abaqus计算误差分析通过误差图显示，最大误差在4%，在误差允许范围内。由此可见，有限元程序和abaqus计算是正确的。2.4 四、装配应力计算2.4.1 处理技巧abaqus与ansys提供了几种由于装配产生的应力的处理方法。耦合，当迫使某节点处

15、多个自由度取得相同的（未知的）某个值时，常用耦合处理，通常用于铰链、销接、外向节等连接处的处理；约束方程，提供了更为通用的联系自由度的方法，使得在某一节点处的自由度满足某个方程（而不是取得相同的值）；当然，对于特殊情况，可用加位移约束实现装配应力的处理。题目给出的十字桁架结构，由于5杆制造时短了一截，建立模型时将3点处建立两个节点（1、2杆对应的是3节点，5杆对应的是4节点，4节点在3节点下方处），则有，其中。不加力时，可以通过材力力法正则方程求得。2.4.2 abaqus处理技巧2.4.3 不加外力（p1,p2,p3）时 材力，ansys与abaqus结果表 5 不加力材力、有限元、abaq

16、us计算结果2.4.4 不加外力（p1,p2,p3）时 材力，ansys与abaqus误差分析表 6 不加力材力、有限元、abaqus计算结果误差分析通过不加力运算结果对比可知，abaqus与ansys中安装应力的处理是正确的。2.4.5 加外力（p1,p2,p3）时 ansys与abaqus结果（1）ansys结果（2）abaqus结果（3）加外力（p1,p2,p3）时 ansys与abaqus结果表 7加力材力、有限元、abaqus计算结果表 8加力材力、有限元、abaqus计算结果误差分析通过加力abaqus与ansys运算结果对比可知，计算应力结果是正确的。2.5 五、优化设计2.5.

17、1 设计中变量的概念 设计变量（dv）：1-10杆面积 i=1,210. 状态变量（sv）：各杆内最大应力max_s小于许用应力，2节点位移小于许用位移 目标函数（obj）：结构杆的总重量最小2.5.2 优化步骤运用vs2012编写复合形法进行约束优化。 复合形法优化原理求解最优化问题的一种算法。该法较为适合解决有约束优化问题。使用该法仅需比较目标函数值即可决定搜索方向，算法较简单，对目标函数的要求不苛刻。复合形是多个单纯形合并成的超多面体，顶点个数ge n+1（n维空间）。复合形法与单纯形法极为相似，却也有不同：1）复合形法不限制顶点个数为n+1，复合形法的顶点个数k取值范围为n+1le k

18、le2n；2）复合形法需要检查顶点的可行性，即是否满足约束。复合形法是由n+1个以上的顶点组合而成的多面体。他的基本思路是：在可行域内构造一初始复合型，然后通过比较各顶点目标函数值，在可行域中找一目标函数值有所改善的新点，并用其替换目标函数值较差的顶点，构成新的复合形。不断重复上述过程，复合形不断变形、转移、缩小，逐渐地逼近最优点。当复合形各顶点目标函数值相差不大或者各顶点相距很近时，则目标函数值最小的顶点即可作为最优点。复合形点点数目k一般取值(n+1)k2n，n是设计变量的个数。为了减小计算变量，复合形法在寻优过程中一般只以在可行域内的反射作为基本搜索策略。复合形法寻优方法主要工作是生成初

19、始复合形和更新复合形。综合来说复合型法的算法思路清晰，容易掌握；不需求导数，不需作一维搜索，对函数性态没有特殊要求；程序结构简单，计算量不大；对初始点要求低，能较快地找到最优解，算法较为可靠。求解时需给出变量取值区间及初始复合形；随着变量维数增多计算效率明显降低；对约束条件较多的非凸问题，常出现多次想形心收缩，使收敛速度减慢。 复合形法优化流程图 vs2012编写复合形法来做约束优化问题执行优化。2.5.3 vs2012优化程序(1) 目标函数（obj）：结构杆的总重量最小public shared function fitness(byval a() as double) as double

20、 dim n as integer n = a.getupperbound(0) 截面面积a dim m as double dim denserty as double denserty = 7.8 / 1000000 dim m1 as double dim l(n) as double dim lt(n, 1) as double l = data.ldata m1 = matrix.vectormultiplevector(a, l) m = denserty * m1 return m end function(2) 主程序，运用复合形法优化全局优化，质量最小 public shar

21、ed function opt(byval ll as double) as double()给定k,a,eps,ndim n as double 点的维度ndim k as double 顶点的数目kn = optimizisiondata.ndatak = optimizisiondata.ddatadim x(k, n) as double 复合形顶点dim xp(n) as double 最优解dim val as doubledim eps, q as doubleeps = 0.000001dim a as doublea = 0.5dim xr(n) as doubledim f

22、r as doubledim t, t1, t2, t3 as doubledim i, j as doubledim f(k) as double 顶点函数值dim u(12), s(n) as doubledim u2 as doubledim ss1, ss2, ss3, ss4 as double初始复合形顶点dim b1, b2 as doubleb1 = 5000b2 = 50000ss4 = 1while ss4 = 1randomize()for i = 1 to kfor j = 1 to nx(i, j) = int(b2 - b1) * rnd() + b1随机初始化位置

23、next jnext i计算顶点函数值()ss1 = 1while ss1 = 1dim xi(n) as doublefor i = 1 to kfor j = 0 to nxi(j) = x(i, j)nextf(i) = optimizision.fitness(xi)next计算好点和坏点 顶点函数值排序dim m(k) as doublefor i = 1 to km(i) = f(i)nextfor j = 1 to kfor i = 1 to k - 1if m(i) m(i + 1) thenelseq = m(i + 1)m(i + 1) = m(i)m(i) = qend

24、ifnext inext jdim l, h, sh as double好点和坏点及次坏点for i = 1 to kif f(i) = m(k) thenl = iend ifif f(i) = m(1) thenh = iend ifif f(i) = m(2) thensh = iend ifnext i是否满足终止条件dim sm as doublesm = 0for j = 1 to ksm = sm + (f(j) - f(l) 2next jdim sm1 as doublesm1 = (sm / k) 0.5 终止条件是否满足终止条件if sm1 b1 xr(i) b2 thent1 = t1 + 1elset1 = 0end ifif s(i) 160 thent3 = t3 + 1elset3 = 0end ifnextif u2 5 thent2 = 1elset2 = 0end ifif t1 = 10 and t2 = 1 and t3 = 10 thent =