十杆桁架结构优化设计

设计过程一、运用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,截面面积20000。并给各杆赋材料属性。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单元坐标系和结构坐标系结构坐标下刚度矩阵：②组装总的刚度矩阵③边界条件处理（固定约束，直接去掉约束对应的行和列）④计算位移向量⑤计算单元应力Vs2012中重要的程序段'计算结点位移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）：结构杆的总重量最小优化步骤运用VS2012编写复合形法进行约束优化。复合形法优化原理求解最优化问题的一种算法。该法较为适合解决有约束优化问题。使用该法仅需比较目标函数值即可决定搜索方向，算法较简单，对目标函数的要求不苛刻。复合形是多个单纯形合并成的超多面体，顶点个数\gen+1（n维空间）。复合形法与单纯形法极为相似，却也有不同：1）复合形法不限制顶点个数为n+1，复合形法的顶点个数k取值范围为n+1\lek\le2n；2）复合形法需要检查顶点的可行性，即是否满足约束。复合形法是由n+1个以上的顶点组合而成的多面体。他的基本思路是：在可行域内构造一初始复合型，然后通过比较各顶点目标函数值，在可行域中找一目标函数值有所改善的新点，并用其替换目标函数值较差的顶点，构成新的复合形。不断重复上述过程，复合形不断变形、转移、缩小，逐渐地逼近最优点。当复合形各顶点目标函数值相差不大或者各顶点相距很近时，则目标函数值最小的顶点即可作为最优点。复合形点点数目k一般取值(n+1)≤k≤2n，n是设计变量的个数。为了减小计算变量，复合形法在寻优过程中一般只以在可行域内的反射作为基本搜索策略。复合形法寻优方法主要工作是生成初始复合形和更新复合形。综合来说复合型法的算法思路清晰，容易掌握；不需求导数，不需作一维搜索，对函数性态没有特殊要求；程序结构简单，计算量不大；对初始点要求低，能较快地找到最优解，算法较为可靠。求解时需给出变量取值区间及初始复合形；随着变量维数增多计算效率明显降低；对约束条件较多的非凸问题，常出现多次想形心收缩，使收敛速度减慢。复合形法优化流程图vs2012编写复合形法来做约束优化问题执行优化。VS2012优化程序(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=1ToNIfXR(I)>B1<>XR(I)<B2ThenT1=T1+1ElseT1=0EndIfIfS(I)<160ThenT3=T3+1ElseT3=0E