十杆桁架结构优化设计

题目：十杆桁架结构优化设计日期：2023.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 Vs2023中重要的程序段 62.3.3 程序输出文献 92.3.4 材料力学、有限元程序、Abaqus结果比较 102.4 四、装配应力计算 112.4.1 解决技巧 112.4.2 Abaqus解决技巧 112.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 优化环节运用VS2023编写复合形法进行约束优化。 142.5.3 VS2023优化程序 162.5.4 优化结果 192.5.5 结果说明 193 设计感想 204 备注 204.1 参考书目 204.2 说明 20设计题目十字桁架结构优化设计现有十字桁架结构见图1，材料泊松比为0.3，E=2.1e11，密度为7.8×103kg/m3,许用应力为160Mpa，P1=600kN，P2=900kN，P3=600kN，杆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，试对结构进行优化，使其重量最小。（同材料力学优化结果比较）。图图SEQ图\*ARABIC1十杆桁架设计过程一、运用Abaqus求解各杆轴力应力运用Abaqus求解，十字桁架结构可用2Dtruss单元模拟。单元参数为：弹性模量，1-6杆截面面积，7-10杆截面面积。加载求解输出各杆应力，结点位移。Abaqus计算流程Part:创建parttrussmain,part45l,part45r,parttrussmain,包涵除8,10杆外的所有杆，part45l包涵8杆，part45r包涵10杆。Property：createMaterial:Elastic:弹性模量,泊松比0.3CreateSection：beam\Truss：SectionA1,截面面积30000。SectionA2,截面面积20230。并给各杆赋材料属性。Assembly：组装parttrussmain,part45l,part45r。Step:创建一个分析步，step1。Interaction：用Tie把parttrussmain,part45l,part45r,绑定。Load：createload：5,6点加铰接约束，固定x,y方向位移。CreateboundaryCondition：2,4点加相应力。Mesh：划分网格，一个杆为一个单元。Elementtape，选trussJob：创建一个job，WriteInput,DataCheck,Submit,通过Result来查看应力云图。结果图图SEQ图\*ARABIC2Abaqus各杆应力云图二、运用材料力学知识求解基本思绪显然题目中的十字桁架结构是两次静不定问题。对于一次静不定问题，材料力学给出了两类解法：①去掉约束加力，找位移协调关系解题；②力法正则方程求解。对于多次静不定，特别是上述桁架问题，找出其协调关系基本上是不也许的，而力法正则方程更适合于解这种结构。如图3所示，去掉多余约束，建立力法正则方程：图图SEQ图\*ARABIC3去多余约束解题过程分别求出外力作用下各杆内力和单位力作用下的各杆内力，为计算方便，将其结果列入下表1中。应用莫尔积分定理有：表1外力作用下各杆内力和单位力作用下的各杆内力杆号LP1P2P3Fi1Fi21a001-12a000013a1-14a0015a000106a000017008009000100000杆号轴力F(N)应力S(MPa)1158056952.672347433.111.583-1119431-37.314474331.585128002.14.276347433.111.5871158850.457.948-962469.9-48.129781447.639.0710-491344.6-24.57 结果表2材料力学各杆应力结果三、编写有限元程序求解程序基本环节计算单元刚度矩阵单元坐标系下刚度矩阵：YYXeYe35xx图图SEQ图\*ARABIC4单元坐标系和结构坐标系结构坐标下刚度矩阵：②组装总的刚度矩阵③边界条件解决（固定约束，直接去掉约束相应的行和列）④计算位移向量⑤计算单元应力Vs2023中重要的程序段'计算结点位移PrivateSubButton1_Click(senderAsObject,eAsEventArgs)HandlesButton1.Click'桁架结点位移计算'形成总刚度矩阵DimTK(12,12)AsDouble'总体刚度矩阵TK=Matrix.STIFFSOfAllTK()DimTKH(11,11)AsDouble'去除0行0列ForI=1To12ForJ=1To12TKH(I-1,J-1)=TK(I,J)NextJNextI'输入结点载荷P(I)DimP(12)AsDoubleP=Data.NodeLoadData()DimPH(11)AsDoubleForI=1To12PH(I-1)=P(I)NextI'边界条件解决ForI=8To11ForJ=1To11TKH(I,J)=0.0NextJNextIForJ=8To11TKH(J,J)=1.0NextJForJ=8To11PH(J)=0.0NextJ'计算结点位移DimZ(11)AsDouble'结点位移DimTKHT(11,11)AsDouble'去除0行0列TKHT=Matrix.InversionOfMatrix(TKH)Z=Matrix.MatrixMultipleVector(TKHT,PH)'输出结点位移IO.Output(Z)EndSubPrivateSubButton2_Click(senderAsObject,eAsEventArgs)HandlesButton2.Click'桁架单元内力'计算总体坐标架单元新节点位移XNEW(6，2)DimXNEW(6,2)AsDoubleDimZNEW(6,2)AsDoubleDimX(6,2)AsDoubleX=Data.PositionData()DimZ(12)AsDouble'结点位移Z=Matrix.Displacement()ForI=1To6ZNEW(I,1)=Z(2*I-2)ZNEW(I,2)=Z(2*I-1)NextIXNEW=Matrix.Add(X,ZNEW)'计算变形后杆长DimDDELTAX(10)AsDoubleDimNEWDDELTAX(10)AsDoubleDimD(10)AsDoubleDimA(10)AsDoubleDimE1AsIntegerA=Data.AreaDataE1=Data.EDataDimN(10)AsDouble'单元内力ForI=1To10X=Data.PositionData()'单元结点编号DimNX(2,10)AsDoubleNX=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)NextI'输出单元内力NIO.Output(N)EndSub'计算单元应力过程PrivateSubButton3_Click(senderAsObject,eAsEventArgs)HandlesButton3.ClickDimS(10)AsDoubleDimF(10)AsDoubleDimA(10)AsDoubleA=Data.AreaDataF=Matrix.FORCEForI=1To10S(I)=F(I)/A(I)NextI'输出单元应力SIO.Output(S)EndSub程序输出文献图5有限元位移结果图6有限元应力结果图7有限元轴力结果材料力学、有限元程序、Abaqus结果比较表3材力、有限元、Abaqus计算结果比较表4材力、有限元、Abaqus计算误差分析通过误差图显示，最大误差在4%，在误差允许范围内。由此可见，有限元程序和Abaqus计算是对的的。四、装配应力计算解决技巧Abaqus与Ansys提供了几种由于装配产生的应力的解决方法。耦合，当迫使某节点处多个自由度取得相同的（未知的）某个值时，常用耦合解决，通常用于铰链、销接、外向节等连接处的解决；约束方程，提供了更为通用的联系自由度的方法，使得在某一节点处的自由度满足某个方程（而不是取得相同的值）；当然，对于特殊情况，可用加位移约束实现装配应力的解决。题目给出的十字桁架结构，由于5杆制造时短了一截Δ，建立模型时将3点处建立两个节点（1、2杆相应的是3节点，5杆相应的是4节点，4节点在3节点下方Δ处），则有，其中。不加力时，可以通过材力力法正则方程求得。Abaqus解决技巧不加外力（P1,P2,P3） 时材力，Ansys与Abaqus结果表5不加力材力、有限元、Abaqus计算结果不加外力（P1,P2,P3） 时材力，Ansys与Abaqus误差分析表6不加力材力、有限元、Abaqus计算结果误差分析通过不加力运算结果对比可知，abaqus与ansys中安装应力的解决是对的的。加外力（P1,P2,P3）时Ansys与Abaqus结果（1）Ansys结果（2）Abaqus结果（3）加外力（P1,P2,P3）时Ansys与Abaqus结果表7加力材力、有限元、Abaqus计算结果表8加力材力、有限元、Abaqus计算结果误差分析通过加力abaqus与ansys运算结果对比可知，计算应力结果是对的的。五、优化设计设计中变量的概念设计变量（DV）：1-10杆面积i=1,210.状态变量（SV）：各杆内最大应力max_s小于许用应力，2节点位移小于许用位移目的函数（OBJ）：结构杆的总重量最小优化环节运用VS2023编写复合形法进行约束优化。复合形法优化原理求解最优化问题的一种算法。该法较为适合解决有约束优化问题。使用该法仅需比较目的函数值即可决定搜索方向，算法较简朴，对目的函数的规定不苛刻。复合形是多个单纯形合并成的超多面体，顶点个数\gen+1（n维空间）。复合形法与单纯形法极为相似，却也有不同：1）复合形法不限制顶点个数为n+1，复合形法的顶点个数k取值范围为n+1\lek\le2n；2）复合形法需要检查顶点的可行性，即是否满足约束。复合形法是由n+1个以上的顶点组合而成的多面体。他的基本思绪是：在可行域内构造一初始复合型，然后通过比较各顶点目的函数值，在可行域中找一目的函数值有所改善的新点，并用其替换目的函数值较差的顶点，构成新的复合形。不断反复上述过程，复合形不断变形、转移、缩小，逐渐地逼近最优点。当复合形各顶点目的函数值相差不大或者各顶点相距很近时，则目的函数值最小的顶点即可作为最优点。复合形点点数目k一般取值(n+1)≤k≤2n，n是设计变量的个数。为了减小计算变量，复合形法在寻优过程中一般只以在可行域内的反射作为基本搜索策略。复合形法寻优方法重要工作是生成初始复合形和更新复合形。综合来说复合型法的算法思绪清楚，容易掌握；不需求导数，不需作一维搜索，对函数性态没有特殊规定；程序结构简朴，计算量不大；对初始点规定低，能较快地找到最优解，算法较为可靠。求解时需给出变量取值区间及初始复合形；随着变量维数增多计算效率明显减少；对约束条件较多的非凸问题，常出现多次想形心收缩，使收敛速度减慢。复合形法优化流程图vs2023编写复合形法来做约束优化问题执行优化。VS2023优化程序(1)目的函数（OBJ）：结构杆的总重量最小PublicSharedFunctionFitness(ByValA()AsDouble)AsDoubleDimnAsIntegern=A.GetUpperBound(0)'截面面积ADimMAsDoubleDimDENSERTYAsDoubleDENSERTY=7.8/1000000DimM1AsDoubleDimL(n)AsDoubleDimLT(n,1)AsDoubleL=Data.LDataM1=Matrix.VectorMultipleVector(A,L)M=DENSERTY*M1ReturnMEndFunction(2)主程序，运用复合形法优化'全局优化，质量最小PublicSharedFunctionOPT(ByValllAsDouble)AsDouble()'给定K,a,Eps,NDimNAsDouble'点的维度NDimKAsDouble'顶点的数目KN=Optimizisiondata.NDataK=Optimizisiondata.DDataDimX(K,N)AsDouble'复合形顶点DimXp(N)AsDouble'最优解DimVALAsDoubleDimEps,QAsDoubleEps=0.000001DimAAsDoubleA=0.5DimXR(N)AsDoubleDimFRAsDoubleDimT,T1,T2,T3AsDoubleDimI,JAsDoubleDimF(K)AsDouble'顶点函数值DimU(12),S(N)AsDoubleDimU2AsDoubleDimSS1,SS2,SS3,SS4AsDouble'初始复合形顶点DimB1,B2AsDoubleB1=5000B2=50000SS4=1WhileSS4=1Randomize()ForI=1ToKForJ=1ToNX(I,J)=Int((B2-B1)*Rnd())+B1'随机初始化位置NextJNextI'计算顶点函数值()SS1=1WhileSS1=1DimXI(N)AsDoubleForI=1ToKForJ=0ToNXI(J)=X(I,J)NextF(I)=Optimizision.Fitness(XI)Next'计算好点和坏点'顶点函数值排序DimM(K)AsDoubleForI=1ToKM(I)=F(I)NextForJ=1ToKForI=1ToK-1IfM(I)>M(I+1)ThenElseQ=M(I+1)M(I+1)=M(I)M(I)=QEndIfNextINextJDimL,H,SHAsDouble'好点和坏点及次坏点ForI=1ToKIfF(I)=M(K)ThenL=IEndIfIfF(I)=M(1)ThenH=IEndIfIfF(I)=M(2)ThenSH=IEndIfNextI'是否满足终止条件DimSMAsDoubleSM=0ForJ=1ToKSM=SM+(F(J)-F(L))^2NextJDimSM1AsDoubleSM1=(SM/K)^0.5'终止条件'是否满足终止条件IfSM1<=EpsThenForI=1ToNXp(I)=X(L,I)NextVAL=F(L)'IO.Output(Xp)'IO.Output(F)Ifll=1ThenReturnXpElseIfll=0ThenReturnFEndIfEndElseSS2=1WhileSS2=1'计算XcDimS1(N)AsDoubleDimS2(N)AsDoubleDimS3(N)AsDoubleDimXc(N)AsDoubleForJ=1ToK'计算XcForI=1ToNS2(I)=X(J,I)NextS1=Matrix.Add(S1,S2)NextForI=1ToNS3(I)=-X(H,I)NextXc=Matrix.Add(S1,S3)ForI=1ToNXc(I)=Xc(I)/(K-1)NextSS3=1WhileSS3=1'计算XRForI=1ToNXR(I)=Xc(I)+A*(Xc(I)-X(H,I))NextIFR=Optimizision.Fitness(XR)'边界条件U=Optimizision.Displacement(XR)S=Optimizision.FORCE(XR)U2=(U(2)^2+U(3)^2)^0.5ForI=1To