题目：十杆桁架结构优化设计日期：2013.09.16目录1设计题目 (1)2设计过程 (2)2.1一、运用Abaqus求解各杆轴力应力 (2)2.1.1Abaqus计算流程 (2)2.1.2结果 (3)2.2二、利用材料力学知识求解 (4)2.2.1基本思路 (4)2.2.2解题过程 (4)2.2.3结果 (5)2.3三、编写有限元程序求解 (6)2.3.1程序基本步骤 (6)2.3.2Vs2012 中重要的程序段 (6)2.3.3程序输出文件 (9)2.3.4材料力学、有限元程序、Abaqus结果比较 (10)2.4四、装配应力计算 (11)2.4.1处理技巧 (11)2.4.2Abaqus处理技巧 (11)2.4.3不加外力（P1,P2,P3）时材力，Ansys与Abaqus结果 (12)2.4.4不加外力（P1,P2,P3）时材力，Ansys与Abaqus误差分析 (12)2.4.5加外力（P1,P2,P3）时 Ansys与Abaqus结果 (12)2.5五、优化设计 (14)2.5.1设计中变量的概念 (14)2.5.2优化步骤运用VS2012编写复合形法进行约束优化。

(14)2.5.3VS2012优化程序 (16)2.5.4优化结果 (20)2.5.5结果说明 (20)3设计感想 (20)4备注 (20)4.1参考书目 (20)4.2说明 (20)十字桁架结构优化设计现有十字桁架结构见图1，材料泊松比为0.3，E=2.1e11，密度为7.8×103kg/m3, 许用应力为160Mpa，P1=600k N ，P2=900k N ，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.005～0.05m2，试对结构进行优化，使其重量最小。

（同材料力学优化结果比较）。

图1十杆桁架2.1 一、运用Abaqus 求解各杆轴力应力利用Abaqus 求解，十字桁架结构可用2Dtruss 单元模拟。

单元参数为：弹性模量112.110E Pa =⨯，1-6杆截面面积210.03A m =，7-10杆截面面积220.02A m =。

加载求解输出各杆应力，结点位移。

2.1.1 Abaqus 计算流程Part:创建part trussmain,part45l,part45r,part trussmain,包涵除8,10杆外的所有杆，part45l 包涵8杆，part45r 包涵10杆。

Property ：create Material:Elastic: 弹性模量112.110E Pa =⨯,泊松比0.3Create Section ：beam\Truss ：Section A1,截面面积30000。

Section A2,截面面积20000。

并给各杆赋材料属性。

Assembly ：组装part trussmain,part45l,part45r 。

Step:创建一个分析步，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结果图 2Abaqus各杆应力云图2.2 二、利用材料力学知识求解 2.2.1 基本思路显然题目中的十字桁架结构是两次静不定问题。

对于一次静不定问题，材料力学给出了两类解法：①去掉约束加力，找位移协调关系解题；②力则方程求解。

对于多次静不定，特别是上述桁架问题，找出其协调关系基本上是不可能的，而力则方程更适合于解这种结构。

如图3所示，去掉多余约束，建立力则方程：2.2.2 解题过程分别求出外力作用下各杆力和单位力作用下的各杆力，为计算方便，将其结果列入下表1中。

应用莫尔积分定理有：杆号 L P1 P2 P3 Fi1 Fi21 a 0 -2P2 0 1 -12 a 0 0 0 0 13 a 1P P2 P3 1 -14 a 0 P2 P3 0 15 a 0 0 0 1 06 a 0 0 0 0 17 2a 21P - 0 02- 2 8 2a 0 2P2 0 2-2 9 2a 0 2P2- 0 02- 10 2a 0 0 0 02- 11111221P11P21P322112222P12P22P3=0=0X X X X δδδδ∆=++∆+∆+∆⎧⎨∆=++∆+∆+∆⎩图 3 去多余约束2.2.3结果杆号轴力F(N) 应力S(MPa)1 1580569 52.672 347433.1 11.583 -1119431 -37.314 47433 1.585 128002.1 4.276 347433.1 11.587 1158850.4 57.948 -962469.9 -48.129 781447.6 39.0710 -491344.6 -24.572.3 三、编写有限元程序求解 2.3.1 程序基本步骤① 计算单元刚度矩阵 单元坐标系下刚度矩阵：1 0 1 00 0 0 01 0 1 00 0 0 0e EA K l --⎡⎤⎢⎥⎢⎥=⎢⎥-⎢⎥⎣⎦结构坐标下刚度矩阵：Te e K T K T -=cos sin 0 0-sin cos 0 00 0 cos sin 0 0 -sin cos T αααααααα⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦②组装总的刚度矩阵③边界条件处理（固定约束，直接去掉约束对应的行和列） ④计算位移向量 ⑤计算单元应力2.3.2 Vs2012 中重要的程序段'计算结点位移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 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) 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 = Matrix.InversionOfMatrix(TKH)Z = Matrix.MatrixMultipleVector(TKHT, PH)'输出结点位移IO.Output(Z)End SubPrivate Sub Button2_Click(sender As Object, e As EventArgs) Handles Button2.Click'桁架单元力'计算总体坐标架单元新节点位移XNEW(6，2)Dim XNEW(6, 2) As DoubleDim ZNEW(6, 2) As DoubleDim X(6, 2) As DoubleX = Data.PositionData()Dim Z(12) As Double'结点位移Z = Matrix.Displacement()For I = 1 To 6ZNEW(I, 1) = Z(2 * I - 2)ZNEW(I, 2) = Z(2 * I - 1)Next IXNEW = Matrix.Add(X, ZNEW)'计算变形后杆长Dim DDELTAX(10) As DoubleDim NEWDDELTAX(10) As DoubleDim D(10) As DoubleDim A(10) As DoubleDim E1 As IntegerA = Data.AreaDataE1 = Data.EDataDim N(10) As Double'单元力For I = 1 To 10X = Data.PositionData()'单元结点编号Dim NX(2, 10) As DoubleNX = 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(I) = NEWDDELTAX(I) - DDELTAX(I)'计算单元力N(I) = D(I) * E1 * A(I) / DDELTAX(I)Next I' 输出单元力 NIO.Output(N)End Sub'计算单元应力过程Private Sub Button3_Click(sender As Object, e As EventArgs) Handles Button3.ClickDim S(10) As DoubleDim F(10) As DoubleDim A(10) As DoubleA = Data.AreaDataF = Matrix.FORCEFor I = 1 To 10S(I) = F(I) / A(I)Next I' 输出单元应力 SIO.Output(S)End Sub2.3.3程序输出文件图 5 有限元位移结果图 6有限元应力结果图 7有限元轴力结果2.3.4材料力学、有限元程序、Abaqus结果比较表 3 材力、有限元、Abaqus计算结果比较表 4 材力、有限元、Abaqus计算误差分析通过误差图显示，最大误差在4%，在误差允许围。