十杆桁架结构优化设计

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

2、4.3 不加外力( P1,P2,P3 )时材力， Ansys与Abaqus 结果 122.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.3 VS2012优化程序162.5.4 优化结果 202.5.5 结果说明 203 设计感想 204 备注 204.1 参考书目 20说明 20十杆桁架结构优化设计1设计题目十字桁架结构优化设计现有十

3、字桁架结构见图 1,材料泊松比为0.3 , E=2.1e11，密度为7.8 X 103kg/m3,许用应力为 160Mpa P仁600k N , P2=900k N , P3=600k N,杆 1-6 面积为A1=0.03m2杆7-10面积为A2=0.02m1、利用Abaqus计算各杆的应力；2、利用材料力学的知识求解，并与1计算出的结果做比较；3、编写有限元程序求解，与1和2计算结果进行比较；4、若杆5制作时短了 0.001m,试求各杆的应力；5、若令2节点的位移小于0.005m, A1、A2为0.0050.05m2，试对结构进行优化，使其重量最小。（同材料力学优化结果比较）。1十杆桁架结构

4、优化设计#十杆桁架结构优化设计2设计过程2.1 一、运用Abaqus求解各杆轴力应力利用Abaqus求解，十字桁架结构可用2Dtruss单元模拟。单元参数为：弹性模量E =2.1 1011Pa， 1-6杆截面面积A =0. Om2， 7-10杆截面面积A =o. om2。加载求解输出各杆应力，结点位移2.1.1Abaqus计算流程Part:仓U建 part trussmain,part45l,part45r,parttrussmain, 包涵除 8,10杆外的所有杆，part45l包涵8杆，part45r包涵10杆。弹性模量E =2.1 1011Pa ,泊松比0.3Property : cre

5、ate Material:Elastic:Mattfial Be-haviorsDataElasticYoung r£Modulus1 210000Poissan'fRatioD,3Create Section:beam'Truss : Section A1,截面面积 30000。Section A2,截面面积20000。并给各杆赋材料属性。=Section ManagerNameAlA2Cnnv.Assembly:组装 part trussmain,part45l,part45rStep :创建一个分析步，step1。3十杆桁架结构优化设计Interaction*

6、Model： I Model-1 Stef绑疋Interaction: 用 Tie 扌巴 part trussmain,part45l,part45r,#十杆桁架结构优化设计#十杆桁架结构优化设计Load: create load : 5,6点加铰接约束，固定 x,y方向位移。Create boundary Condition: 2,4 点加相应力。Mesh划分网格，一个杆为一个单元。Element tape，选trussElement Type#十杆桁架结构优化设计#十杆桁架结构优化设计FamilyPiezoelectricPipeThermal ElectricTrussS, S11 (A

7、vg: 7 5%)+ 5.794e+01 + 4.910e+01 + 4.026e+01 + 3.143e+01 + 2.259e+01 r +1.37 5e+01 k +4.910e+00 k -3.929e+00 k -1.277e+01 k -2.161e+0L k -3.045e+01 L -3.928e+01 L -4.812e+01Element ILibraryo Standard ExplicitGeo metric OrderQ Linear ! QuadraticJob:创建一个 job，Write Input,Data Check,Submit, 通过 Result 来查

8、看应力云图。结果#十杆桁架结构优化设计#十杆桁架结构优化设计图2 Abaqus各杆应力云图#十杆桁架结构优化设计5十杆桁架结构优化设计2.2二、利用材料力学知识求解2.2.1 基本思路显然题目中的十字桁架结构是两次静不定问题。对于一次静不定问题，材料力学给出了两类解法：去掉约束加力，找位移协调关系解题；力法正则方程 求解。对于多次静不定，特别是上述桁架冋题，找出其协调关系基本上是不可能 的，而力法正则方程更适合于解这种结构。如图 3所示，去掉多余约束，建立力 法正则方程：、Xi| "qX 2 ： 1P1 ： 1P2* ： 1P3 0、21 X1、- 22X 2 ' '

9、 2P1 ' ' 2P2 ' ' 2P3 _02.2.2 解题过程分别求出外力作用下各杆内力和单位力作用下的各杆内力，为计算方便，将其结果列入下表1中。应用莫尔积分定理有：表1外力作用下各杆内力和单位力作用下的各杆内力杆号LP1P2P3Fi1Fi21a0-2P201-12a000013aP1P2P31-14a0P2P3015a0 t00106a000017辰-72P100-近7280V2P20429>/2a0-V2P200-42100000-巫eal込4忑Q = 6.563 X 10"EA2EA1+ 4 = -5.9915x10-EA 1 EA

10、2込工=工10】尸110=z10=E1EA,F丽EA：F :帀EAjEAlPKEA1 JP2=£ F：S2P210=2 1P3EAt=£F :g2 P3EA.£ F而EA.2半P f = J9974EAl'忑P $ 二 1.7974EA22/2 pi 2.69618EA1 EA22P尹=十EAlEA2= 0.342S4EA1杆号轴力F(N)应力S(MPa)1158056952.672347433.111.583-1119431-37.31447433,.585128002.14.276347433.1,1.5871158850.457.948-962469

11、.9-48.129781447.639.0710-491344.6-24.57结果223表2材料力学各杆应力结果2.3三、编写有限兀程序求解2.3.1程序基本步骤计算单元刚度矩阵-10-101- EA0000Ke-l-10100000 一结构坐标下刚度矩阵:Ke KeTcossi n：0cos000cos：sin 0-sin ：cos：0-sin :0 组装总的刚度矩阵 边界条件处理（固定约束，直接去掉约束对应的行和列） 计算位移向量 计算单元应力2.3.2 Vs2012中重要的程序段As Object , e As Eve ntArgs) Han dies总体刚度矩阵'计算结点位移P

12、rivate Sub Butto n1_Click(se nder Butto n1.Click'桁架结点位移计算'形成总刚度矩阵Dim TK(12, 12) As DoubleTK = Matrix .STIFFSOfAIITK()Dim TKH(11, 11) As Double '去除 0行0列For I = 1 To 12For J = 1 To 12TKH(I - 1, J - 1) = TK(I, J)Next JNext I'输入结点载荷P(I)Dim P(12) As DoubleP = Data.NodeLoadData()Dim PH(11)

13、 As DoubleFor I = 1 To 12PH(I - 1) = P(I)Next I'边界条件处理For I = 8 To 11For J = 1 To 11TKH(I, J) = 0.0Next JNext IFor J = 8 To 11TKH(J, J) = 1.0Next JFor J = 8 To 11PH(J) = 0.0Next J'计算结点位移Dim Z(11) As Double'结点位移Dim TKHT(11, 11) As Double '去除 0行0列TKHT = MatrixnversionOfMatrix(TKH)Z = M

14、atrix .MatrixMultipleVector(TKHT, PH)'输出结点位移IO.Output(Z)End SubPrivate Sub Button2_Click(sender As Object , e As EventArgs) Handles Butto n2.Click'桁架单元内力'计算总体坐标架单元新节点位移 XNEW(6 2)Dim XNEW(6, 2) As DoubleDim ZNEW(6, 2) As DoubleDim X(6, 2) As DoubleX =Data.Positio nData()Dim Z(12) As Doubl

15、e'结点位移Z =Matrix .Displacement()For I = 1 To6ZNEW(I, 1) = Z(2 * I - 2)ZNEW(I, 2) = Z(2 * I - 1)Next IXNEW = Matrix .Add(X, ZNEW)'计算变形后杆长Dim DDELTAX(IO) As DoubleDim NEWDDELTAX(10As DoubleDim D(10) As DoubleDim A(10) As DoubleDim E1 As IntegerA =Data.AreaDataE1 =Data.EDataDim N(10) As Double&#

16、39;单元内力For I = 1 To 10X =Data.Positio nData()'单元结点编号Dim NX(2, 10) As DoubleNX =Data.NodeData()NEWDDELTAX(I) = Math.Sqrt(XNEW(NX(1, I),1) - XNEW(NX(2,I),1) A 2 + (XNEW(NX(1, I), 2) - XNEW(NX(2, I), 2) A 2)DDELTAX(I) =Matrix .ElementLongger(l)D(I) = NEWDDELTAX(I) - DDELTAX(I)'计算单元内力N(I) = D(I)

17、 * E1 * A(I) / DDELTAX(I)Next I'输出单元内力NIO.Output(N)End Sub'计算单元应力过程Private Sub Butt on 3_Click(se nderAs Object , e As Eve ntArgs) Ha ndlesButto n3.ClickDim S(10) As DoubleDim F(10) As DoubleDim A(10) As DoubleA = Data.AreaDataF = Matrix .FORCEFor I = 1 To 10 S(I) = F(I) / A(I)Next I'输出单

18、元应力SIO.Output(S)End Sub233 程序输出文件髓吃羿00L102 -OO5J1O -000.6B-005*409000.903-002.553-000* 640-002.626000.000000.000000.000000.000图5有限元位移结果1轴比”记亭本丈件旧 虽盲 摘式g）_至看辿_wh（h） 000,0001582153.030349149. 291 -1117754.194 49315.135 128581.338 348147. 949 1159282. 297-961985. 978782220.608 -490587. 585勘黑釧F髓勺醱7靄VI豔I

19、H：052.738OIL 638 -037.258 00L644001286OIL 605051964-048.099 min-024.529图6有限元应力结果图7有限元轴力结果234 材料力学、有限元程序、Abaqus结果比较ABAQUSrVSp材力计算结果对比80III III11十杆桁架结构优化设计-6012345绘囲区i 7&g10 A 5|MPa>52.691158-37.311.S84.2757J94-4B.1239.07M.57 VS$(MPaS2.73Bll.W-3Z2581.6444.21C.11.60557.964-4S.D3939.111-24.529材力S

20、(MPa52-57115B-37.311.584.2711.5857.M4B-123SUD7-245? ASIMPak VS 5 MPa) 材力 S(|VIPa表3材力、有限元、Abaqus计算结果比较abaqus. V&2012与材力误差对比0.0450.04D50632S-0.0003797330-0.00139(1723 0.D4D5Q63W0I30JJ 5*895Q .O.0PMH222 D-0.0004364W-0 QQ16&57QZ1Q屿村力邊差 VS与村力溟差o.oosooasC.D0215B8950.000414222-0.00043S4® 0.001

21、049399-0001666702#十杆桁架结构优化设计#十杆桁架结构优化设计屿时力误差晡与村力邊差表4材力、有限元、Abaqus计算误差分析通过误差图显示，最大误差在4%在误差允许范围内。由此可见，有限元#十杆桁架结构优化设计#十杆桁架结构优化设计程序和Abaqus计算是正确的。#十杆桁架结构优化设计2.4四、装配应力计算2.4.1 处理技巧Abaqus与Ansys提供了几种由于装配产生的应力的处理方法。耦合，当迫 使某节点处多个自由度取得相同的（未知的）某个值时，常用耦合处理，通常用 于铰链、销接、外向节等连接处的处理；约束方程，提供了更为通用的联系自由 度的方法，使得在某一节点处的自由度

22、满足某个方程（而不是取得相同的值）；当然，对于特殊情况，可用加位移约束实现装配应力的处理。题目给出的十字桁架结构，由于 5杆制造时短了一截，建立模型时将3 点处建立两个节点（1、2杆对应的是3节点，5杆对应的是4节点，4节点在3 节点下方处），则有-UY3 +UY4 =也，其中A=1。不加力时，可以通过材力力法正则厲八Mi，仁乂2：1方程2 二-2iX O22X2 =0求得。f Comtreint 临叫y EUlLWmmiLe2.4.2Abaqus处理技巧Name2 Constraint-?Const naint-57 ConstraintJlV Constaii価V Constnaint-3

23、7 Constraint*Name： Constraiit-7Typsi EqjaborevidentSet NanifDOFQ Conslraint'711 1141Name* 伽¥ CN5trairt-2lP Cwshd祁/级畑iirt-"Enter cne mw : dab 喻 each lem in lhe 呷痂亍Ckkrnous? bLtton J fcr table :ptor5.CSVS ID伽岫(u obal；#十杆桁架结构优化设计13十杆桁架结构优化设计不加力时abaqu$T材力结果对比O#十杆桁架结构优化设计#十杆桁架结构优化设计4512a456

24、7g1村力帥网4.7054.1954.7D54.29594.295>-9.9B-9.9S-9411-9 111 Abaqus SMpa|i4.70464.29494.70464.2949B.99954 294 S-9.981-g.112 ANSYS 5 （Mpd）470524.7W24 2K5BOW?4 2W5-9.9812-9W12-9.11221122科力 Abuqu 5 ANSVS S （ Mp»表5不加力材力、有限元、Abaqus计算结果244 不加外力（P1,P2,P3） 时 材力，Ansys与Abaqus误差分析村力与辰误差-S.5CU5E-ffi -2.32S29

25、E-ffi -S.50159f-EK -2.a2829E-ffi -5.55556E4S -2.32E291E-ffi 0.D001D02O.aDOJOOZ Q-CQaiD9757 D.DDD109757科力与Aft巧嗾差 4J.25DK-CB O.OOOUB414 4.25OSE-D5 0.D0011fefll4 7.7777BE-05 0.000116414 0.00012024 0_00012024 0_«H13；17D9 D.DD013171M枕tl与硼砂溟差村力2酗林溟雀表6不加力材力、有限元、Abaqus计算结果误差分析 通过不加力运算结果对比可知，abaqus与ansy

26、s中安装应力的处理是正确2.4.5 加外力（P1,P2,P3）时 Ansys 与 Abaqus结果（1）Ansys 结果Elie 注iec* Lisi E"i PEmC好h 呂口匕曲1 世士 FagiueieTfl亡皿凸rM He-Lt-* Tj 圃” 3Sj It EllAKY5 TaolbErSA1EJ：创 隆创1_睡qjT| PafflSRPH,AKYS .3 leroS TrcferenctB： E rr-Hprwunswr Bl 5ra| ill i nnB Ci=3iirj-ia>l P-uh I pruc HUnln il TileB Etad RcdliI tf

27、i FmlsrB Criteria日 Plot flriful tv ffl Dn Friran'-d Shagair 曰匚IICll out PH u E ZUHEU! Q Eluant S3u Qtle> Table QLioe Elm Rce3 Tnclnr Flol3 Pst Fml h I EiRh3 EuciEr-i: Cc Fl E313 ThLidFElMB List Jk：Gul tfi Que ft Result口 Oplfcr (hi4pS Rrnil tv rats FGR 卩和. e Vadfci £.al“=IleKnc TAbLe H Pa

28、i h Over at Lons Surf-Bee- Opernt imic i Load «de=Clwck Elv» E Vlri 1 e Knni E h=14(Ml DpCFHl XEDUC-JjKOU. fiJUHTOH51XF-2 两-LTDT-JUW (IMS)Kn-jEX 4.4M5T see记-“SEP £ MUL4l4IM¥3. £15.«J71751.921523.035B74.25022吓刮o5创F-中IIcsyECHET1)I Pick a xrtj Ltea «* enter anSm躬世T(2)

29、 Abaqus 结果S, Sil (Avg： 75%)+ 5.744e+01+ 4.782O+0 丄+ 3.819A+01+ 2.8 56e+01 + 1.894e+01 + 9.312e+00 -3.143e-01-Q.QdOp+OO-1.957e+01 -2.919e+01-3.882e+01 -4844吕十 0i -5,807e+01(3) 加外力(P1,P2,P3 )时 Ansys 与 Abaqus 结果ABCAnsys 力口力 S (MPa)Abaqus 应力呂(KPa>Abaqus Ansys 结果误差57. 39157. 4416-0. 000880895115. S771

30、5. 9365-0,003733568-32. 609-32.56070.001433383一5.8767S.9369-6 010139972)13 26813 2895-0. 001617819L15. 87715. 8991-0. 00139001647. 96147.9607& 25512E-06-58. 105-53.07010.00060099329. 9629 9945-0. 001150211-33. 68-33 64220.00112353915十杆桁架结构优化设计-B012345678910Msystfl力 SIMPfl)5739115S77證图区k朗鬥13 2631

31、5.87747,961-58.10529,96-33.68Abaqus应力5 (MPa)57.441615.9365-32-5607 5.93691389515.89914Z9607-58.070129.9945-33.6422血那加力S(MPa) AMquS应力$ ( MPa)80604020020-4060表7加力材力、有限元、Abaqus计算结果Abaqus与Airys结果误差-D.006D.0DO6O(m&.&.2S512E4J6. 1-0>00aSECS95JLJL 一患-D.OO 咖O.Q014435B5-或DOS-0.01-0 010139972M2 123

32、45679!帥勒3 勻尿移 结果误差DQOBNS95-£JK3735«3 a.0ClX4E 3-3fi3<-O.01D139972-D.DD1617319-O CC 13>90016 6 2551J E-Dg D DDKD099B-£.0ai 150211 Q.0011235B9表8加力材力、有限元、Abaqus计算结果误差分析通过加力abaqus与ansys运算结果对比可知，计算应力结果是正确的2.5五、优化设计2.5.1 设计中变量的概念 设计变量（DV）: 1-10杆面积A i=1,210. 状态变量（SV）：各杆内最大应力max_s小于许用应力

33、，2节点位移小于 许用位移 目标函数（OBJ:结构杆的总重量最小2.5.2 优化步骤运用VS2012编写复合形法进行约束优化。复合形法优化原理求解最优化问题的一种算法。该法较为适合解决有约束优化问题。使用该法仅需比较目标函数值即可决定搜索方向，算法较简单，对目标 函数的要求不苛刻。复合形是多个单纯形合并成的超多面体，顶点个数 ge n+1（ n 维空 间）。复合形法与单纯形法极为相似，却也有不同：1）复合形法不限制顶点个数为n+1,复合形法的顶点个数k取值范 围为 n+1le kle2n ；2）复合形法需要检查顶点的可行性,即是否满足约束。 复合形法是由 n+1 个以上的顶点组合而成的多面体。

34、他的基本思路是：在可行域内构造一初始复合型,然后通过比较各顶点目标函数值, 在可行域中找一目标函数值有所改善的新点,并用其替换目标函数值较 差的顶点,构成新的复合形。不断重复上述过程,复合形不断变形、转 移、缩小,逐渐地逼近最优点。当复合形各顶点目标函数值相差不大或 者各顶点相距很近时,则目标函数值最小的顶点即可作为最优点。复合 形点点数目k 一般取值（n+1） <k<2n, n是设计变量的个数。为了减小 计算变量,复合形法在寻优过程中一般只以在可行域内的反射作为基本 搜索策略。复合形法寻优方法主要工作是生成初始复合形和更新复合形。 综合来说复合型法的算法思路清晰,容易掌握；不需求

35、导数,不需作一 维搜索,对函数性态没有特殊要求；程序结构简单,计算量不大；对初 始点要求低,能较快地找到最优解,算法较为可靠。求解时需给出变量 取值区间及初始复合形；随着变量维数增多计算效率明显降低；对约束 条件较多的非凸问题,常出现多次想形心收缩,使收敛速度减慢。17十杆桁架结构优化设计复合形法优化流程图2.5.3 VS2012优化程序 目标函数(OBJ:结构杆的总重量最小Public Shared Function Fitness( ByVai A() As Double) As DoubleDim n As Integern = A.GetUpperBou nd(0)'截面面积A

36、Dim M As DoubleDim DENSERTYKs DoubleDENSERTY = 7.8 / 1000000Dim M1 As DoubleDim L(n) As DoubleDim LT(n, 1) As DoubleL = Data.LDataM1 = Matrix .VectorMultipleVector(A, L)M = DENSERTY * M1Retur n MEnd Function(2)主程序，运用复合形法优化 全局优化，质量最小19十杆桁架结构优化设计#十杆桁架结构优化设计Public Shared Function OPT(ByValII As Double)

37、 As Double()'给定 K,a,Eps,NDim N As Double'点的维度NDim K As Double'顶点的数目KN = Optimizisi on data.NDataK = Optimizisi on data.DDataDim X(K, N) As Double '复合形顶点Dim Xp(N) As Double '最优解Dim VAL As DoubleI Dim Eps, Q As DoubleEps = 0.000001Dim A As DoubleA = 0.5Dim XR(N) As DoubleDim FR As

38、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 = 5000_!L_B2 = 50000SS4 = 1While SS4 = 1Ran domize()For I = 1 To KFor J = 1 To NX(I, J) = In t(B2 - B1) * Rn d()

39、+B1'随机初始化位置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) = Optimizisio n .Fit ness(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 =

40、M(I + 1)M(I + 1) = M(IM(I) = QEnd IfNext INext JDim L, H, SH As Double'好点和坏点及次坏点End IfNext I'是否满足终止条件Dim SM As DoubleSM = 0For J = 1 ToKSM = SM + (F(J) - F(L) A 2Next JDim SM1 As DoubleSM1 = (SM / K) a 0.5'终止条件I'是否满足终止条件If SM1 <= Eps ThenFor I = 1 ToNXp(l) = X(L, I)NextVAL = F(L)I

41、O.Output(Xp)'IO.Output(F)If II = 1ThenReturn XpElseIf II = 0ThenRetur n FIEnd IfEndElseSS2= 1While SS2 = 1'计算XcDim S1(N) As DoubleDim S2(N) As DoubleDim S3(N) As DoubleDim Xc(N) As DoubleFor J = 1 To K '计算 XcFor I = 1 To NS2(I) = X(J, I)NextS1 = Matrix .Add(S1, S2)NextFor I = 1 To NS3(I) = -X(H, I)NextXc = Matrix .Add(S1, S3)For I = 1 To NXc(I) = Xc(I) / (K - 1)NextSS3= 1While SS3 = 1'计算XRFor I = 1 To NXR(I) = Xc(I) + A * (Xc(I) - X(H,I) INext IFR = Optimizis