题目:十杆桁架结构优化设计日期：2013。
09.16目录1 设计题目 (1)2 设计过程 (2)2.1 一、运用Abaqus求解各杆轴力应力 (2)2.1。
1 Abaqus计算流程 (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。
2 Vs2012 中重要的程序段 (6)2。
3.3 程序输出文件 (9)2.3.4 材料力学、有限元程序、Abaqus结果比较 (10)2.4 四、装配应力计算 (11)2。
4。
1 处理技巧 (11)2.4。
2 Abaqus处理技巧 (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。
3 VS2012优化程序 (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 0 2- 10 2a0 0 0 0 2- 11111221P11P21P322112222P12P22P3=0=0X X X X δδδδ∆=++∆+∆+∆⎧⎨∆=++∆+∆+∆⎩图 3 去多余约束2.2.3结果杆号轴力F(N)应力S(MPa）152.67 2347433.111.583-1119431—37。
31 447433 1.58 5128002。
1 4.27 6347433.111.58 71158850.457.94 8—962469.9-48.12 9781447。
639.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％,在误差允许范围内。