有限元法大作业

合集下载

有限单元法课后习题全部答案_王勖成

有限单元法课后习题全部答案_王勖成


∂ 2φ ∂ 2φ ∂φ ∂φ k − ∫ k 2 + k 2 + Q δφ d Ω + ∫ δφ d Ω − ∫ αφ − q − k δφ d Γ Ω Γ − Γ Γ q q ∂y ∂n ∂n ∂x
欧拉方程: k
∂ 2φ ∂ 2φ + +Q = k 0 ∂x 2 ∂y 2
习题 1.2: 在用有限元法求解时,边界条件总是满足的,控制方程的不完全匹配,会产生误差。题中所 ,代入边 给出的近似函数: φ =a0 + a1 x + a2 x + a3 x ,应该满足边界条件,对于情况(1)
2 3
界条件可得 = a0 0, = a3
1 − a1 L − a2 L2 ,从而 L3 x3 x3 x3 2 ) + a ( x − ) + 2 L2 L L3

= =
∑{ A
m k =1 m
T
( N j ( xk )) [ A( N i ( xk )ai ) − f ( xk )]
m
}
( N j )A( N i )ai − ∑ AT ( N j ) f = k 1= k 1
T
∑A
= Ka-P
(写成矩阵形式)
因此, kij =
d 2 w dw d 3 w 0 dx 2 δ dx − dx3 δ w = 0
L
1.5 如有一问题的泛函为 = Π ( w)

L
0
EI d 2 w 2 kw2 + qwdx ,其中 E, I, k 是常数,q 2 + 2 dx 2

(完整版)有限元大作业matlab---课程设计例子

(完整版)有限元大作业matlab---课程设计例子

有限元大作业程序设计学校:天津大学院系:建筑工程与力学学院专业:01级工程力学姓名:刘秀学号:\\\\\\\\\\\指导老师:连续体平面问题的有限元程序分析[题目]:如图所示的正方形薄板四周受均匀载荷的作用,该结构在边界上受正向分布压力,m kNp 1=,同时在沿对角线y 轴上受一对集中压力,载荷为2KN ,若取板厚1=t ,泊松比0=v 。

[分析过程]:由于连续平板的对称性,只需要取其在第一象限的四分之一部分参加分析,然后人为作出一些辅助线将平板“分割”成若干部分,再为每个部分选择分析单元。

采用将此模型化分为4个全等的直角三角型单元。

利用其对称性,四分之一部分的边界约束,载荷可等效如图所示。

[程序原理及实现]:用FORTRAN程序的实现。

由节点信息文件NODE.IN和单元信息文件ELEMENT.IN,经过计算分析后输出一个一般性的文件DATA.OUT。

模型基本信息由文件为BASIC.IN生成。

该程序的特点如下:问题类型:可用于计算弹性力学平面问题和平面应变问题单元类型:采用常应变三角形单元位移模式:用用线性位移模式载荷类型:节点载荷,非节点载荷应先换算为等效节点载荷材料性质:弹性体由单一的均匀材料组成约束方式:为“0”位移固定约束,为保证无刚体位移,弹性体至少应有对三个自由度的独立约束方程求解:针对半带宽刚度方程的Gauss消元法输入文件:由手工生成节点信息文件NODE.IN,和单元信息文件ELEMENT.IN结果文件:输出一般的结果文件DATA.OUT程序的原理如框图:(1)主要变量:ID:问题类型码,ID=1时为平面应力问题,ID=2时为平面应变问题N_NODE:节点个数N_LOAD:节点载荷个数N_DOF:自由度,N_DOF=N_NODE*2(平面问题)N_ELE:单元个数N_BAND:矩阵半带宽N_BC:有约束的节点个数PE:弹性模量PR:泊松比PT:厚度LJK_ELE(I,3):单元节点编号数组,LJK_ELE(I,1),LJK_ELE(I,2),LJK_ELE(I,3)分别放单元I的三个节点的整体编号X(N_NODE), Y(N_NODE):节点坐标数组,X(I),Y(I)分别存放节点I的x,y 坐标值P_LJK(N_BC,3):节点载荷数组,P_LJK(I,1)表示第I个作用有节点载荷的节点的编号,P_LJK(I,2),P_LJK(I,3)分别为该节点沿x,y方向的节点载荷数值AK(N_DOF,N_BAND):整体刚度矩阵AKE(6,6):单元刚度矩阵BB(3,6):位移……应变转换矩阵(三节点单元的几何矩阵)DD(3,3):弹性矩阵SS(3,6);应力矩阵RESULT_N(N_NOF):节点载荷数组,存放节点载荷向量,解方程后该矩阵存放节点位移DISP_E(6)::单元的节点位移向量STS_ELE(N_ELE,3):单元的应力分量STS_ND(N_NODE,3):节点的应力分量(2)子程序说明:READ_IN:读入数据BAND_K:形成半带宽的整体刚度矩阵FORM_KE:计算单元刚度矩阵FORM_P:计算节点载荷CAL_AREA:计算单元面积DO_BC:处理边界条件CLA_DD:计算单元弹性矩阵SOLVE:计算节点位移CLA_BB:计算单元位移……应变关系矩阵CAL_STS:计算单元和节点应力(3)文件管理:源程序文件:chengxu.for程序需读入的数据文件:BASIC.IN,NODE.IN,ELEMENT.IN(需要手工生成)程序输出的数据文件:DATA.OUT(4)数据文件格式:需读入的模型基本信息文件BASIC.IN的格式如下表需读入的节点信息文件NODE.IN的格式如下表需读入的单元信息文件ELEMENT.IN的格式如下表输出结果文件DATA.OUT格式如下表[算例原始数据和程序分析]:(1)模型基本信息文件BASIC.IN的数据为1,4,6,5,31.,0.,1.1,1,0,2,1,0,4,1,1,5,0,1,6,0,11,-0.5,-1.5,3.,-1.,-1,6,-0.5,-0.5(2)手工准备的节点信息文件NODE.IN的数据为1 0.0 2.02 0.0 1.03 1.0 1.04 0. 0.5 1.0 0.6 2.0 0.(3)手工准备的单元信息文件ELEMENT.IN的数据为1 2 3 3 0 0 0 0 1 1 1 1 0 12 4 5 5 0 0 0 0 1 1 1 1 0 25 3 2 2 0 0 0 0 1 1 1 1 0 33 5 6 6 0 0 0 0 1 1 1 1 04 (4)源程序文件chengxu.for为:PROGRAM FEM2DDIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),&RESULT_N(500),AK(500,100)D IMENSION STS_ELE(500,3),STS_ND(500,3)OPEN(4,FILE='BASIC.IN')OPEN(5,FILE='NODE.IN')OPEN(6,FILE='ELEMENT.IN')OPEN(8,FILE='DATA.OUT')OPEN(9,FILE='FOR_POST.DAT')READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOADIF(ID.EQ.1)WRITE(8,20)IF(ID.EQ.2)WRITE(8,25)20 FORMAT(/5X,'=========PLANE STRESS PROBLEM========')25 FORMAT(/5X,'=========PLANE STRAIN PROBLEM========')CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT, & IJK_ELE,X,Y,IJK_U,P_IJK)CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,& IJK_ELE,X,Y,PE,PR,PT,AK)CALL FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)CALL SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, & STS_ELE,STS_ND)c to putout a data fileWRITE(9,70)REAL(N_NODE),REAL(N_ELE)70 FORMAT(2f9.4)WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I),& STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE)71 FORMA T(7F9.4)WRITE(9,72)(REAL(IJK_ELE(I,1)),REAL(IJK_ELE(I,2)),&REAL(IJK_ELE(I,3)),REAL(IJK_ELE(I,3)),&STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1, N_ELE)72 FORMAT(7f9.4)cCLOSE(4)CLOSE(5)CLOSE(6)CLOSE(8)CLOSE(9)E NDcc to get the original data in order to model the problemSUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR, &PT,IJK_ELE,X,Y,IJK_U,P_IJK)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3), & P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14)REAL ND_ANSYS(N_NODE,3)READ(4,*)PE,PR,PTREAD(4,*)((IJK_U(I,J),J=1,3),I=1,N_BC)READ(4,*)((P_IJK(I,J),J=1,3),I=1,N_LOAD)READ(5,*)((ND_ANSYS(I,J),J=1,3),I=1,N_NODE)READ(6,*)((NE_ANSYS(I,J),J=1,14),I=1,N_ELE)DO 10 I=1,N_NODEX(I)=ND_ANSYS(I,2)Y(I)=ND_ANSYS(I,3)10 CONTINUEDO 11 I=1,N_ELEDO 11 J=1,3IJK_ELE(I,J)=NE_ANSYS(I,J)11 CONTINUEN_BAND=0DO 20 IE=1,N_ELEDO 20 I=1,3DO 20 J=1,3IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J))IF(N_BAND.LT.IW)N_BAND=IW20 CONTINUEN_BAND=(N_BAND+1)*2IF(ID.EQ.1) THENELSEPE=PE/(1.0-PR*PR)PR=PR/(1.0-PR)END IFR ETURNENDcC to form the stiffness matrix of elementSUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3), & AKE(6,6), SS(6,6)CALL CAL_DD(PE,PR,DD)CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 10 I=1,3DO 10 J=1,6SS(I,J)=0.0DO 10 K=1,310 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 20 I=1,6DO 20 J=1,6AKE(I,J)=0.0DO 20 K=1,320 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PTRETURNENDcc to form banded global stiffness matrixSUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE, & PR,PT,AK)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100)N_DOF=2*N_NODEDO 40 I=1,N_DOFDO 40 J=1,N_BAND40 AK(I,J)=0DO 50 IE=1,N_ELECALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DO 50 I=1,3DO 50 II=1,2IH=2*(I-1)+IIIDH=2*(IJK_ELE(IE,I)-1)+IIDO 50 J=1,3DO 50 JJ=1,2IL=2*(J-1)+JJIZL=2*(IJK_ELE(IE,J)-1)+JJIDL=IZL-IDH+1IF(IDL.LE.0) THENELSEAK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL)END IF50 CONTINUERETURNENDcc to calculate the area of elementSUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=(XIJ*YIK-XIK*YIJ)/2.0RETURNENDcc to calculate the elastic matrix of elementSUBROUTINE CAL_DD(PE,PR,DD)DIMENSION DD(3,3)DO 10 I=1,3DO 10 J=1,310 DD(I,J)=0.0DD(1,1)=PE/(1.0-PR*PR)DD(1,2)=PE*PR/(1.0-PR*PR)DD(2,1)=DD(1,2)DD(2,2)=DD(1,1)DD(3,3)=PE/((1.0+PR)*2.0)RETURNENDcc to calculate the strain-displacement matrix of elementSUBROUTINE CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)DO 10 II=1,3DO 10 JJ=1,310 BB(II,JJ)=0.0BB(1,1)=Y(J)-Y(K)BB(1,3)=Y(K)-Y(I)BB(1,5)=Y(I)-Y(J)BB(2,2)=X(K)-X(J)BB(2,4)=X(I)-X(K)BB(2,6)=X(J)-X(I)BB(3,1)=BB(2,2)BB(3,2)=BB(1,1)BB(3,3)=BB(2,4)BB(3,4)=BB(1,3)BB(3,5)=BB(2,6)BB(3,6)=BB(1,5)CALL CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DO 20 I1=1,3DO 20 J1=1,620 BB(I1,J1)=BB(I1,J1)/(2.0*AE)RETURNENDcc to form the global load matrixSUBROUTINE FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3), & RESULT_N(N_DOF)DO 10 I=1,N_DOF10 RESULT_N(I)=0.0DO 20 I=1,N_LOADII=P_IJK(I,1)RESULT_N(2*II-1)=P_IJK(I,2)20 RESULT_N(2*II)=P_IJK(I,3)RETURNENDcc to deal with BC(u) (here only for fixed displacement) using "1-0" method SUBROUTINE DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100)DO 30 I=1,N_BCIR=IJK_U(I,1)DO 30 J=2,3IF(IJK_U(I,J).EQ.0)THENELSEII=2*IR+J-3AK(II,1)=1.0RESULT_N(II)=0.0DO 10 JJ=2,N_BAND10 AK(II,JJ)=0.0DO 20 JJ=2,II20 AK(II-JJ+1,JJ)=0.0END IF30 CONTINUERETURNENDcc to solve the banded FEM equation by GAUSS eliminationSUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),AK(500,100)DO 20 K=1,N_DOF-1IF(N_DOF.GT.K+N_BAND-1)IM=K+N_BAND-1IF(N_DOF.LE.K+N_BAND-1)IM=N_DOFDO 20 I=K+1,IML=I-K+1C=AK(K,L)/AK(K,1)IW=N_BAND-L+1DO 10 J=1,IWM=J+I-K10 AK(I,J)=AK(I,J)-C*AK(K,M)20 RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K)RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1)DO 40 I1=1,N_DOF-1I=N_DOF-I1IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1IF(N_BAND.LE.N_DOF-I-1)JQ=N_BANDDO 30 J=2,JQK=J+I-130 RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)40 RESULT_N(I)=RESULT_N(I)/AK(I,1)WRITE(8,50)50 FORMAT(/12X,'* * * * * RESULTS BY FEM2D * * * * *',//8X,&'--DISPLACEMENT OF NODE--'//5X,'NODE NO',8X,'X-DISP',8X,'Y-DISP') DO 60 I=1,N_NODE60 WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I)70 FORMAT(8X,I5,7X,2E15.6)RETURNENDcc calculate the stress components of element and nodeSUBROUTINECAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, &STS_ELE,STS_ND)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),DD(3,3),BB(3,6), &SS(3,6),RESULT_N(N_DOF),DISP_E(6)DIMENSION STS_ELE(500,3),STS_ND(500,3)WRITE(8,10)10 FORMAT(//8X,'--STRESSES OF ELEMENT--')CALL CAL_DD(PE,PR,DD)DO 50 IE=1,N_ELECALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 20 I=1,3DO 20 J=1,6SS(I,J)=0.0DO 20 K=1,320 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 30 I=1,3DO 30 J=1,2IH=2*(I-1)+JIW=2*(IJK_ELE(IE,I)-1)+J30 DISP_E(IH)=RESULT_N(IW)STX=0STY=0TXY=0DO 40 J=1,6STX=STX+SS(1,J)*DISP_E(J)STY=STY+SS(2,J)*DISP_E(J)40 TXY=TXY+SS(3,J)*DISP_E(J)STS_ELE(IE,1)=STXSTS_ELE(IE,2)=STYSTS_ELE(IE,3)=TXY50 WRITE(8,60)IE,STX,STY,TXY60 FORMAT(1X,'ELEMENT NO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)c the following part is to calculate stress components of nodeWRITE(8,55)55 FORMAT(//8X,'--STRESSES OF NODE--')DO 90 I=1,N_NODEA=0.B=0.C=0.II=0DO 70 K=1,N_ELEDO 70 J=1,3IF(IJK_ELE(K,J).EQ.I) THENII=II+1A=A+STS_ELE(K,1)B=B+STS_ELE(K,2)C=C+STS_ELE(K,3)END IF70 CONTINUESTS_ND(I,1)=A/IISTS_ND(I,2)=B/IISTS_ND(I,3)=C/IIWRITE(8,75)I,STS_ND(I,1),STS_ND(I,2),STS_ND(I,3)75 FORMAT(1X,'NODE NO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)90 CONTINUERETURNENDc FEM2D programm end[算例结果]:chengxu.for所输出的数据文件DATA.OUT数据内容如下:=========PLANE STRESS PROBLEM========* * * * * RESULTS BY FEM2D * * * * *--DISPLACEMENT OF NODE--NODE NO X-DISP Y-DISP1 .000000E+00 -.525275E+012 .000000E+00 -.225275E+013 -.108791E+01 -.137363E+014 .000000E+00 .000000E+005 -.824176E+00 .000000E+006 -.182418E+01 .000000E+00--STRESSES OF ELEMENT--ELEMENT NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00ELEMENT NO.= 2STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00ELEMENT NO.= 3STX=-.108791E+01 STY=-.137363E+01 TXY= .307692E+00ELEMENT NO.= 4STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00--STRESSES OF NODE--NODE NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00NODE NO.= 2STX=-.100000E+01 STY=-.220879E+01 TXY= .249084E+00NODE NO.= 3STX=-.105861E+01 STY=-.191575E+01 TXY= .205128E+00NODE NO.= 4STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00NODE NO.= 5STX=-.970696E+00 STY=-.166667E+01 TXY= .586081E-01NODE NO.= 6STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00[结论与体会]:通过本次的课程设计,我对有限元的概念有了更加深刻的理解,同时也弥补了平时学习是疏忽的地方,充实了有限元知识。

有限元分析大作业报告

有限元分析大作业报告

有限元分析大作业报告一、引言有限元分析是工程领域中常用的数值模拟方法,通过将连续的物理问题离散为有限个子区域,然后利用数学方法求解,最终得到数值解。

有限元分析的快速发展和广泛应用,为工程领域提供了一种强大的工具。

本报告将介绍在大作业中所进行的有限元分析工作及结果。

二、有限元模型建立本次大作业的研究对象是工程结构的应力分析。

首先,通过对结构进行几何建模,确定了结构的尺寸和形状。

然后,将结构离散为有限个单元,每个单元又可以看作一个小的子区域。

接下来,为了求解结构的应力分布,需要为每个单元确定适当的单元类型和单元属性。

最后,根据结构的边界条件,建立整个有限元模型。

三、材料属性和加载条件在建立有限元模型的过程中,需要为材料和加载条件确定适当的参数。

本次大作业中,通过实验获得了结构材料的弹性模量、泊松比等参数,并将其输入到有限元模型中。

对于加载条件,我们选取了其中一种常见的加载方式,并将其施加到有限元模型中。

四、数值计算和结果分析为了求解结构的应力分布,需要进行数值计算。

在本次大作业中,我们选用了一种常见的有限元求解器进行计算。

通过输入模型的几何形状、材料属性和加载条件,求解器可以根据有限元方法进行计算,并得到结构的应力分布。

最后,我们通过对计算结果进行分析,得出了结论。

五、结果讨论和改进方法根据计算结果,我们可以对结构的应力分布进行分析和讨论。

根据分析结果,我们可以得出结论是否满足设计要求以及结构的强度情况。

同时,根据分析结果,我们还可以提出改进方法,针对结构的特点和问题进行相应的优化设计。

六、结论通过对工程结构进行有限元分析,我们得到了结构的应力分布,并根据分析结果进行了讨论和改进方法的提出。

有限元分析为工程领域提供了一种有效的数值模拟方法,可以帮助工程师进行结构设计和分析工作,提高设计效率和设计质量。

【1】XXX,XXXX。

【2】XXX,XXXX。

以上是本次大作业的有限元分析报告,总结了在建立有限元模型、确定材料属性和加载条件、数值计算和结果分析等方面的工作,并对计算结果进行讨论和改进方法的提出。

有限元分析大作业

有限元分析大作业

有限元大作业一题目要求:图1所示为一悬臂梁,在端部承受载荷,材料弹性模量为E,泊松比为1/3,悬臂梁的厚度(板厚)为t,若该粱被划分为两个单元,单元和节点编号如图所示,试按平面应力问题计算各个节点位移计支反力。

一、单元划分1.计算简图及单元划分如下所示:2.进行节点及单元编号节点i j m单元① 2 3 4② 3 2 13.节点坐标值节点号1 2 3 4坐标值X 2 2 0 0Y 1 0 1 0二、计算单元刚度矩阵1、计算每个单元面积△以及i b ,i c (m j i i ,,=) ①②单元的面积相等,即12121=⨯⨯=∆ 单元①的i b ,i c⎩⎨⎧=--==-=0)(1m j i m j i y x c y y b ⎩⎨⎧=--==-=2)(0i m ji m j x x c y y b ⎩⎨⎧-=--=-=-=2)(1j i mj i m y x c y y b 对平面应力问题,其表达式为[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+-+-+-+∆-=s r s r sr s r s r s r s r s r b b uc c cb u b uc b c u c ub c c u b b u Et Krs 21212121)1(42 然后对单元①求解单元刚度子矩阵2==i r 2==i s []⎥⎦⎤⎢⎣⎡=3/1001329)1(22Et K 2==i r 3==j s []⎥⎦⎤⎢⎣⎡=03/23/20329)1(23Et K2==i r 4==m s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)1(24Et K 3==j r 3==j s []⎥⎦⎤⎢⎣⎡=4003/4329)1(33Et K 3==j r 2==i s []⎥⎦⎤⎢⎣⎡=03/23/20329)1(32Et K 3==j r 4==m s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)1(34Et K 4==m r 4==m s []⎥⎦⎤⎢⎣⎡=3/133/43/43/7329)1(44Et K 4==m r 2==i s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)1(42Et K 4==m r 3==j s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)1(43Et K由子矩阵[]e rs K 合成单元刚度矩阵[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------=3/133/443/23/13/23/43/73/23/43/2143/24003/23/23/403/43/203/13/203/23/103/213/2001329)1(Et K将单元①的单元刚度矩阵补零升阶变为单元刚度矩阵,其在总体刚度矩阵中的位置为:节点号→单元②的i b ,i c⎩⎨⎧=--=-=-=0)(1m j im j i y x c y y b ⎩⎨⎧-=--==-=2)(0i m ji m j x x c y y b ⎩⎨⎧=--==-=2)(1j i mj i m y x c y y b 然后对单元 求解单元刚度子矩阵:3==i r 3==i s []⎥⎦⎤⎢⎣⎡=3/1001329)2(33Et K 3==i r 2==j s []⎥⎦⎤⎢⎣⎡=03/23/20329)2(32Et K 3==i r 1==m s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)2(31Et K 1 2 3 412[])1(22K[])1(23K[])1(24K3[])1(32K[])1(33K[])1(34K4[])1(42K[])1(43K[])1(44K2==j r 2==j s []⎥⎦⎤⎢⎣⎡=4003/4329)2(22Et K 2==j r 3==i s []⎥⎦⎤⎢⎣⎡=03/23/20329)2(23Et K 2==j r 1==m s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)2(21Et K 1==m r 1==m s []⎥⎦⎤⎢⎣⎡=3/133/43/43/7329)2(11Et K 1==m r 3==i s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)2(13Et K 1==m r 2==j s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)2(12Et K 由子矩阵[]e rs K 合成单元刚度矩阵[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------=3/133/443/23/13/23/43/73/23/43/2143/24003/23/23/403/43/203/13/203/23/103/213/2001329)2(Et K将单元②的单元刚度矩阵补零升阶变为单元贡献矩阵,其在总体刚度矩阵中的位置为:节点号→1 2 3 41 [])2(11K[])2(12K[])2(13K2 [])2(21K[])2(22K[])2(23K3 [])2(31K [])2(32K [])2(33K 4三、计算总体刚度矩阵总体刚度矩阵是由各单元的贡献矩阵迭加而成)2()1(][][][][K K K K e +==∑四、进行节点约束处理根据节点约束情况,在总刚矩阵中可采用划行划列处理约束的方法,由题目易知,节点3和4的已知水平位移和垂直位移都为零,划去其相对应的行和列,则总刚矩阵由8阶变为4阶,矩阵如下:⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧--=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------2/02/03/13043/203/73/23/443/23/133/43/23/43/43/73292211p p v u v u Et329][Et K =1 2 3 413/133/43/43/743/23/23/4----3/13/23/21----000243/23/23/4----3/13003/73/43/403/13/23/21----33/13/23/21----3/43/403/13003/743/23/23/4----40003/13/23/21----43/23/23/4----3/133/43/43/7化简⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧--=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------Et p Et p v u v u 3/1603/160130122072412213424472211 五、求解线性方程组方法:采用LU 分解法 1.求解矩阵[]U 各元素⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-----⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------75/10775/640075/6475/353007/767/27/7502447~7/877/87/7607/87/337/207/767/27/7502447~13012207241221342447⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----353/44900075/6475/353007/767/27/7502447~ 得到的[]U 矩阵如下:[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=353/44900075/6475/353007/767/27/7502447U 2.求解矩阵[]L 各元素⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----13012207241221342447353/44900075/6475/353007/767/27/75024471353/6475/767/20175/27/40017/40001 得到的[]L 矩阵如下:[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------=13012207241221342447L3.进行求解⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧---=⇒⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧--=Et p Et p Et p y Et p Et p Ly 79425/850800225/323/1603/1603/160⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧---=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----⇒=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡Et p Et p Et p v u v u y v u v u U 79425/850800225/323/160353/44900075/6475/353007/7675/27/750244722112211 解得Et p v /422.82-= Et p u /497.12-= Et p v /028.91-= Et p u /897.11=于是求得各节点的位移为:⎩⎨⎧-==Etp v Etp u /028.9/897.111 ⎩⎨⎧-=-=Etp v Etp u /422.8/497.122 ⎩⎨⎧==033v u ⎩⎨⎧==044v u 六、求解相应的支反力(运用静力学的平衡方程进行求解)3号节点和4号节点的支反力如下图所示:。

有限元大作业报告

有限元大作业报告

有限元计算分析报告***********院&……*设计专业2008年六月目录试题一 (1)试题二 (4)试题三 (6)试题四 (8)试题一问题描述:图示为一带圆孔的单位厚度的正方形平板,在x 方向作用均布压力0.25Mpa ,试用三节点常应变单元和六节点三角形单元对平板进行有限元分析。

图1.(1)分析:(1)该平板的受力属于平面应力问题,在对称外载荷的作用下,模型的受力也是对称的。

故而对该问题也可以简化为4/1平板倒圆角的模型。

在对称均布压力作用下,平板的内部应力图也应该是对称的,平板受到沿x 方向的正应力x σ和y 方向的应力y σ,板面上没有力的作用,即:0=zy τ,0=zx τ。

垂直于板面没有正应力作用,z σ=0。

即该数学模型与z 无关,仅仅是x 、y 的函数。

(2)简化为41平板计算时,可以将倒圆左边界视作x 约束而y 方向无约束,下边界视作y 方向有约束而x 方向没有约束。

分别用三节点和六节点画出单元网格。

(3) 由于平板的中间孔存在集中应力,所以在孔的附近的有限元网格需要细化,而远离孔的网格就可以不画那么细了。

有限元网格划分结果如下图1.(2)所示;数学建模:按照1/4计算,取点(0,0,0),由矢量(0.024,0.024,0)通过surface 创建平面,再建立2DArcangle 画出90度圆弧,将平面打断删去多余部分便得到了几何模型。

六节点应力图六节点应变图三节点应力图三节点应变图三节点不同数目网格应力图三节点不同数目网格数目应变图试题二问题描述:图示 2.(1)为带方孔(边长为80mm)的悬臂梁,其上受部分均布载荷(p=10Kn/m)作用,试采用一种平面单元,对图示两种结构进行有限元分析,并就方孔的布置进行分析比较,如将方孔设计为圆孔,结果有何变化?(板厚为1mm,材料为钢)分析:(1)该题同第一个问题一样是属于平面问题,是平面的应变问题。

(2)平面及面内无z方向应力分量,且限制了z向位移。

(完整word版)有限元分析大作业报告要点

(完整word版)有限元分析大作业报告要点

有限元分析大作业报告试题1:一、问题描述及数学建模图示无限长刚性地基上的三角形大坝,受齐顶的水压力作用,试用三节点常应变单元和六节点三角形单元对坝体进行有限元分析,并对以下几种计算方案进行比较:(1)分别采用相同单元数目的三节点常应变单元和六节点三角形单元计算;(2)分别采用不同数量的三节点常应变单元计算;(3)当选常应变三角单元时,分别采用不同划分方案计算。

该问题属于平面应变问题,大坝所受的载荷为面载荷,分布情况及方向如图所示。

二、采用相同单元数目的三节点常应变单元和六节点三角形单元计算1、有限元建模(1)设置计算类型:两者因几何条件和载荷条件均满足平面应变问题,故均取Preferences 为Structural(2)选择单元类型:三节点常应变单元选择的类型是Solid Quad 4 node182;六节点三角形单元选择的类型是Solid Quad 8 node183。

因研究的问题为平面应变问题,故对Element behavior(K3)设置为plane strain。

(3)定义材料参数:弹性模量E=2.1e11,泊松比σ=0.3(4)建几何模型:生成特征点;生成坝体截面(5)网格化分:划分网格时,拾取lineAB和lineBC,设定input NDIV 为15;拾取lineAC,设定input NDIV 为20,选择网格划分方式为Tri+Mapped,最后得到600个单元。

(6)模型施加约束:约束采用的是对底面BC 全约束。

大坝所受载荷形式为Pressure ,作用在AB 面上,分析时施加在L AB 上,方向水平向右,载荷大小沿L AB 由小到大均匀分布。

以B 为坐标原点,BA 方向为纵轴y ,则沿着y 方向的受力大小可表示为:}{*980098000)10(Y y g gh P -=-==ρρ2、 计算结果及结果分析 (1) 三节点常应变单元三节点常应变单元的位移分布图三节点常应变单元的应力分布图(2)六节点三角形单元六节点三角形单元的变形分布图六节点三角形单元的应力分布图①最大位移都发生在A点,即大坝顶端,最大应力发生在B点附近,即坝底和水的交界处,且整体应力和位移变化分布趋势相似,符合实际情况;②结果显示三节点和六节点单元分析出来的最大应力值相差较大,原因可能是B点产生了虚假应力,造成了最大应力值的不准确性。

有限元分析大作业

有限元分析大作业

《有限元分析及应用》大作业——齿根弯曲应力计算报告班级:无可奉告姓名:无可奉告学号:无可奉告指导老师:无可奉告目录目录 (2)1.概述 (3)1.1工程问题描述 (3)1.2问题分析 (3)2.建模过程 (4)2.1几何建模 (4)2.2CAE网格划分与计算 (5)2.3后处理 (8)3.多方案比较与结果分析 (9)3.1多方案比较 (9)3.2结果分析 (11)1.概述1.1工程问题描述我在本次作业中的选题为齿根弯曲应力的计算与校核。

通过对机械设计的学习,我们可以知道,齿轮的失效形式主要是齿面接触疲劳和齿根弯曲断裂,而闭式传动硬齿面齿轮的失效形式以齿根弯曲断裂,这个时候进行齿根弯曲应力的校核才比较有意义,在设计问题的时候应当选取这种类型的算例。

设计计算的另一个主要思路是将有限元计算的结果与传统机械设计的结算结果进行对比,以从多方面验证计算结果的准确性。

综上,我们最终选取了《机械原理》(第三版)P50例3-1中的问题进行校核计算。

已知起重机械用的一对闭式直齿圆柱齿轮,传动,输入转速n1=730r/min,输入功率P1=35kW,每天工作16小时,使用寿命5年,齿轮为非对称布置,轴的刚性较大,原动机为电动机,工作机载荷为中等冲击。

z1=29,z2=129,m=2.5mm,b1=48mm,b2=42mm,大、小齿轮均为20CrMnTi,渗碳淬火,齿面硬度为58~62HRC,齿轮精度为7级,试验算齿轮强度。

齿面为硬齿面,传动方式为闭式传动。

根据设计手册查出的许用接触应力为1363.6Mpa,计算结果为1260Mpa,强度合格。

根据设计手册查出的许用弯曲应力为613.3MPa,计算结果为619Mpa,强度略显不够。

1.2问题分析大小齿轮啮合,小齿轮受载荷情况较为严峻,故分析对象应当为小齿轮。

可以看出,由于齿轮单侧受载荷,传动过程中每个齿上载荷的变化过程是相同的,故问题可被简化为反对称问题,仅需研究单个齿。

有限元法理论及应用参考答案

有限元法理论及应用参考答案

有限元法理论及应用大作业1、试简要阐述有限元理论分析的基本步骤主要有哪些?答:有限元分析的主要步骤主要有:(1)结构的离散化,即单元的划分;(2)单元分析,包括选择位移模式、根据几何方程建立应变与位移的关系、根据虚功原理建立节点力与节点位移的关系,最后得到单元刚度方程;(3)等效节点载荷计算;(4)整体分析,建立整体刚度方程;(5)引入约束,求解整体平衡方程。

2、有限元网格划分的基本原则是什么?指出图示网格划分中不合理的地方。

题2图答:一般选用三角形或四边形单元,在满足一定精度情况,尽可能少一些单元。

有限元划分网格的基本原则:1.拓扑正确性原则。

即单元间是靠单元顶点、或单元边、或单元面连接2.几何保持原则。

即网络划分后,单元的集合为原结构近似3.特性一致原则。

即材料相同,厚度相同4.单元形状优良原则。

单元边、角相差尽可能小5.密度可控原则。

即在保证一定精度的前提下,网格尽可能的稀疏一些。

(a)(b)中节点没有有效的连接,且(b)中单元边差相差很大。

(c)中没有考虑对称性,单元边差很大。

3、分别指出图示平面结构划分为什么单元?有多少个节点?多少个自由度?题3图答:(a )划分为杆单元, 8个节点,12个自由度。

(b )划分为平面梁单元,8个节点,15个自由度。

(c )平面四节点四边形单元,8个节点,13个自由度。

(d )平面三角形单元,29个节点,38个自由度。

4、什么是等参数单元?。

答:如果坐标变换和位移插值采用相同的节点,并且单元的形状变换函数与位移插值的形函数一样,则称这种变换为等参变换,这样的单元称为等参单元。

5、在平面三节点三角形单元中,能否选取如下的位移模式,为什么?(1).⎪⎩⎪⎨⎧++=++=26543221),(),(y x y x v yx y x u αααααα (2). ⎪⎩⎪⎨⎧++=++=2652423221),(),(yxy x y x v yxy x y x u αααααα 答:(1)不能,因为位移函数要满足几何各向同性,即单元的位移分布不应与人为选取的 坐标方位有关,即位移函数中的坐标x,y 应该是能够互换的。

有限元法大作业

有限元法大作业

有限元法大作业一平面刚架的程序用Visual C++编制的平面刚架的源程序如下:///////////////////////////////////////////////////////程序开始//////////////////////////////////////////////////////////////////#include"iostream.h"#include"math.h"#include"stdlib.h"#include"conio.h"//*****************//声明必要变量//*****************#define PI 3.141592654int NE; //单元数int NJ; //节点数int NZ; //支承数int NPJ; //有节点载荷作用的节点数int NPF; //非节点载荷数int HZ; //载荷码int E; //单元码int fangchengshu;double F[303]; //各节点等效总载荷数值int dym_jdm[100][2]; //单元码对应的节点码:dym_jdm[][0], dym_jdm[][1]分为前后节点总码int zhichengweizhi[300]; //记录支持节点作用点的数组int fjzhzuoyongdanyuan[100]; //非节点载荷作用单元int fjzhleixing[100]; //非节点载荷类型:1-均布,2-垂直集中,3-平行集中,4-力偶,5-角度集中double fjzhzhi[100]; //非节点载荷的值double fjzhzuoyongdian[100]; //非节点载荷在各竿的作用点double fjzhjiaodu[100]; //非节点载荷作用角度int jdzhzuoyongdian[100]; //节点载荷作用的节点数组double jiedianzaihe[101][3];//节点载荷值,其jiedianzaihe[][0]-- jiedianzaihe[][2]分别为U, V, Mdouble zhengtigangdu[303][303]; //整体刚度数组double changdu[100]; //各单元竿长数组double jiaodu[100]; //各单元角度数组double tanxingmoliang[100]; //各单元弹性模量数组double J_moliang[100]; //各单元J模量数组double mianji[100]; //各单元面积数组double weiyi[303]; //记录各个节点位移的数组double dy_weiyi[100][6]; //各个单元在局部坐标系中的位移数组dy_weiyi[i][0]-dyweiyi[i][6]分别为第i+1单元的u1,v1,@1,u2,v2,@2double dy_neili[100][6];//各个单元在局部坐标系中的固端内力dy_weiyi[i][0]-dyweiyi[i][6]分别为第i+1单元的U1,V1,M1,U2,V2,M2double gan_neili[100][6];//各个单元的竿端内力数组,gan_neili[i][6]表示第i+1单元的6内力.//*******************//*******************void input(); //数据的输入void zonggang(); //计算总刚度,存放于zhengtigangdu[][]数组中void zongzaihe(); //计算等效总节点载荷void zhichengyinru(); //引入支承条件void jsweiyi(); //求各个节点位移void js_dy_weiyi(); //求局部坐标系中的位移void ganduanneili(); //求竿端内力void dy_gangdu(int i,double dg[6][6]); //求单元在局部坐标系中的单刚void js_T_T1(int i,double T[6][6],double T1[6][6]); //求单元的转换矩阵及其逆阵//************//主函数//************void main(){input();cout<<" 输出结果"<<endl;cout<<"==============================”<<endl;zonggang();zongzaihe();zhichengyinru();jsweiyi();cout<<"输出位移:"<<endl;for( int i=0;i<NJ;i++){cout<<"u["<<i+1<<"]:"<<weiyi[i*3]<<endl;cout<<"v["<<i+1<<"]:"<<weiyi[i*3+1]<<endl;cout<<"@["<<i+1<<"]:"<<weiyi[i*3+2]<<endl;}js_dy_weiyi();ganduanneili();//输出竿端内力cout<<"==============================”<<endl;cout<<"输出竿内力:"<<endl;for( i=0;i<NE;i++){cout<<"第"<<i+1<<"单元:"<<endl;cout<<" U1:"<<gan_neili[i][0]<<" V1:"<<gan_neili[i][01]<<" M1:"<<gan_neili[i][2]<<endl;cout<<" U2:"<<gan_neili[i][3]<<" V2:"<<gan_neili[i][4]<<" M2:"<<gan_neili[i][5]<<endl;}cout<<" ==========================================="<<endl;cout<<"输出完毕,按任意键退出!"<<endl;getchi();}//***************************//***************************void dy_gangdu(int i,double dg[6][6]){dg[0][0]=tanxingmoliang[i]*mianji[i]/changdu[i]; dg[0][1]=0; dg[0][2]=0;dg[0][3]=-tanxingmoliang[i]*mianji[i]/changdu[i]; dg[0][4]=0; dg[0][5]=0; dg[1][0]=0;dg[1][1]=12*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]*changdu[i]);dg[1][2]=-6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]); dg[1][3]=0;dg[1][4]=-12*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]*changdu[i]);dg[1][5]=-6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]); dg[2][0]=0;dg[2][1]=-6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]);dg[2][2]=4*tanxingmoliang[i]*J_moliang[i]/changdu[i]; dg[2][3]=0;dg[2][4]=6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]);dg[2][5]=2*tanxingmoliang[i]*J_moliang[i]/changdu[i];dg[3][0]=-tanxingmoliang[i]*mianji[i]/changdu[i]; dg[3][1]=0; dg[3][2]=0;dg[3][3]=tanxingmoliang[i]*mianji[i]/changdu[i]; dg[3][4]=0; dg[3][5]=0; dg[4][0]=0;dg[4][1]=-12*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]*changdu[i]);dg[4][2]=6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]); dg[4][3]=0;dg[4][4]=12*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]*changdu[i]);dg[4][5]=6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]); dg[5][0]=0;dg[5][1]=-6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]);dg[5][2]=2*tanxingmoliang[i]*J_moliang[i]/changdu[i]; dg[5][3]=0;dg[5][4]=6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]);dg[5][5]=4*tanxingmoliang[i]*J_moliang[i]/changdu[i];}//*************************//求单元的转换矩阵及其逆阵//*************************void js_T_T1(int i,double T[6][6],double T1[6][6]){int Tti,Tti2;for( Tti=0;Tti<6;Tti++){//计算单元的转换矩阵Tfor( Tti2=0;Tti2<6;Tti2++){if(Tti==0&&Tti2==0)T[Tti][Tti2]=cos(jiaodu[i]);else if(Tti==0&&Tti2==1)T[Tti][Tti2]=sin(jiaodu[i]);else if(Tti==1&&Tti2==0)T[Tti][Tti2]=-sin(jiaodu[i]);else if(Tti==1&&Tti2==1)T[Tti][Tti2]=cos(jiaodu[i]);else if(Tti==2&&Tti2==2)T[Tti][Tti2]=1;T[Tti][Tti2]=cos(jiaodu[i]);else if(Tti==3&&Tti2==4)T[Tti][Tti2]=sin(jiaodu[i]);else if(Tti==4&&Tti2==3)T[Tti][Tti2]=-sin(jiaodu[i]);else if(Tti==4&&Tti2==4)T[Tti][Tti2]=cos(jiaodu[i]);else if(Tti==5&&Tti2==5)T[Tti][Tti2]=1;elseT[Tti][Tti2]=0;if(fabs(T[Tti][Tti2])<0.0000001)T[Tti][Tti2]=0;if(T[Tti][Tti2]>0.99999999)T[Tti][Tti2]=1;if(T[Tti][Tti2]<-0.99999999)T[Tti][Tti2]=-1;}}//计算T的转置矩阵T1for(Tti=0;Tti<6;Tti++){for(Tti2=0;Tti2<6;Tti2++)T1[Tti2][Tti]=T[Tti][Tti2];}}//*******************//计算竿端内力的函数//*******************void ganduanneili(){int Ti,Ti2;double chengji=0;double danyuan_gangdu[6][6]; //单元刚度数组for(int i=0;i<NE;i++){dy_gangdu(i,danyuan_gangdu);for(Ti=0;Ti<6;Ti++){for(Ti2=0;Ti2<6;Ti2++)chengji=chengji+danyuan_gangdu[Ti][Ti2]*dy_weiyi[i][Ti2];gan_neili[i][Ti]=chengji+dy_neili[i][Ti];chengji=0;}}//************************************************************//计算各个单元在局部坐标系中的位移的函数,存放在dy_weiyi[][]中//************************************************************void js_dy_weiyi(){int Ti,Ti3;double T[6][6],T1[6][6];double chengji=0;for(int i=0;i<NE;i++){js_T_T1(i,T,T1);int qianjiedian,houjiedian;qianjiedian=dym_jdm[i][0];houjiedian=dym_jdm[i][1];for(Ti=0;Ti<6;Ti++){dy_weiyi[i][Ti]=T[Ti][0]*weiyi[(qianjiedian-1)*3]+T[Ti][1]*weiyi[(qianjiedian-1)*3+1]+T[Ti][2]* weiyi[(qianjiedian-1)*3+2]+T[Ti][3]*weiyi[(houjiedian-1)*3]+T[Ti][4]*weiyi[(houjiedian-1)*3+1]+T[Ti][5]* weiyi[(houjiedian-1)*3+2];}}}//***********************//计算得位移的函数//***********************void jsweiyi(){int Ti,Ti2;int i,j,k;int i_max;double max,zhongjian;double gd_kz[303][303];for(Ti=0;Ti<NJ*3;Ti++){for(Ti2=0;Ti2<NJ*3;Ti2++)gd_kz[Ti][Ti2]=zhengtigangdu[Ti][Ti2];gd_kz[Ti][NJ*3]=F[Ti];}//用高斯主列消去法解方程,先计算增广矩阵的变形矩阵for(Ti=0;Ti<NJ*3-1;Ti++){max=fabs(gd_kz[Ti][Ti]);i_max=Ti;for(Ti2=Ti;Ti2<NJ*3;Ti2++){{i_max=Ti2-1;max=gd_kz[Ti2][Ti];}}if(i_max!=Ti){for(i=Ti;i<NJ*3+1;i++){zhongjian=gd_kz[Ti][i];gd_kz[Ti][i]=gd_kz[i_max][i];gd_kz[i_max][i]=zhongjian;}}for(j=Ti+1;j<NJ*3;j++){double shang=gd_kz[j][Ti]/gd_kz[Ti][Ti];for(k=0;k<NJ*3+1;k++)gd_kz[j][k]=-shang*gd_kz[Ti][k]+gd_kz[j][k];}}//反代求值double chengji=0;weiyi[NJ*3-1]=gd_kz[NJ*3-1][NJ*3]/gd_kz[NJ*3-1][NJ*3-1];for(i=NJ*3-2;i>=0;i--){for(j=i+1;j<NJ*3;j++)chengji=chengji+gd_kz[i][j]*weiyi[j];weiyi[i]=(gd_kz[i][NJ*3]-chengji)/gd_kz[i][i];chengji=0;}}//********************//引入支承条件的函数//********************void zhichengyinru(){int i,Ti;int weizhi;fangchengshu=NJ*3-NZ;for(i=0;i<NZ;i++){weizhi=zhichengweizhi[i]-1;for(Ti=0;Ti<NJ*3;Ti++)zhengtigangdu[weizhi][Ti]=0;zhengtigangdu[Ti][weizhi]=0;zhengtigangdu[weizhi][weizhi]=1;F[weizhi]=0;}}//*************************//计算等效总节点载荷的函数//*************************void zongzaihe(){int i,Ti;double U1=0,V1=0,M1=0,U2=0,V2=0,M2=0;double c,G,l,d;double Ff[6];double T[6][6],T1[6][6]; //T为坐标转换矩阵,T1为T的转置矩阵for(i=0;i<NJ*3;i++)F[i]=0;for(i=0;i<6;i++)Ff[i]=0;for( i=0;i<NE;i++){for(int jl=0;jl<6;jl++)dy_weiyi[i][jl]=0;}for(i=0;i<NPF;i++){int dym=fjzhzuoyongdanyuan[i];U1=V1=M1=U2=V2=M2=0;for(int t=0;t<6;t++)Ff[t]=0;c=fjzhzuoyongdian[i];l=changdu[(fjzhzuoyongdanyuan[i]-1)];d=l-c;G=fjzhzhi[i];if(fjzhleixing[i]==1){V1=-0.5*G*c*(2-2*c*c/(l*l)+c*c*c/(l*l*l));V2=-G*c-V1;M1=G*c*c*(6-8*c/l+3*c*c/(l*l))/12;M2=-G*c*c*c*(4-3*c/l)/(12*l);U1=U2=0;}else if(fjzhleixing[i]==2){V1=-G*(l+2*c)*d*d/(l*l*l);M1=G*c*d*d/(l*l);M2=-G*c*c*d/(l*l);U1=U2=0;}else if(fjzhleixing[i]==3){U1=-G*d/l;U2=-G*c/l;M1=M2=V1=V2=0;}else if(fjzhleixing[i]==4){U1=U2=0;V1=-6*G*c*d/(l*l*l);V2=-V1;M1=G*d*(2-3*d/l)/l;M2=G*c*(2-3*c/l)/l;}else if(fjzhleixing[i]==5){V1=-G*sin(fjzhjiaodu[i])*(l+2*c)*d*d/(l*l*l);V2=-G*sin(fjzhjiaodu[i])*(l+2*d)*c*c/(l*l*l);M1=G*sin(fjzhjiaodu[i])*c*d*d/(l*l);M2=-G*sin(fjzhjiaodu[i])*c*c*d/(l*l);U1=U2=0;U1=-G*cos(fjzhjiaodu[i])*d/l+U1;U2=-G*cos(fjzhjiaodu[i])*c/l+U2;}//记录竿的固端内力dy_neili[dym-1][0]=U1;dy_neili[dym-1][1]=V1;dy_neili[dym-1][2]=M1;dy_neili[dym-1][3]=U2;dy_neili[dym-1][4]=V2;dy_neili[dym-1][5]=M2;js_T_T1( fjzhzuoyongdanyuan[i]-1, T,T1);for(Ti=0;Ti<6;Ti++)Ff[Ti]=-(T1[Ti][0]*U1+T1[Ti][1]*V1+T1[Ti][2]*M1+T1[Ti][3]*U2+T1[Ti][4]*V2+T1[Ti][5]*M2);//将载荷转换到整体坐标系中F[(dym_jdm[dym-1][0]-1)*3]=F[(dym_jdm[dym-1][0]-1)*3]+Ff[0];F[(dym_jdm[dym-1][0]-1)*3+1]=F[(dym_jdm[dym-1][0]-1)*3+1]+Ff[1];F[(dym_jdm[dym-1][0]-1)*3+2]=F[(dym_jdm[dym-1][0]-1)*3+2]+Ff[2];F[(dym_jdm[dym-1][1]-1)*3]=F[(dym_jdm[dym-1][1]-1)*3]+Ff[3];F[(dym_jdm[dym-1][1]-1)*3+1]=F[(dym_jdm[dym-1][1]-1)*3+1]+Ff[4];F[(dym_jdm[dym-1][1]-1)*3+2]=F[(dym_jdm[dym-1][1]-1)*3+2]+Ff[5];for(i=0;i<NPJ;i++){F[(jdzhzuoyongdian[i]-1)*3]=F[(jdzhzuoyongdian[i]-1)*3]+jiedianzaihe[i][0];F[(jdzhzuoyongdian[i]-1)*3+1]=F[(jdzhzuoyongdian[i]-1)*3+1]+jiedianzaihe[i][1];F[(jdzhzuoyongdian[i]-1)*3+2]=F[(jdzhzuoyongdian[i]-1)*3+2]+jiedianzaihe[i][2];}}//***********************************************//计算总刚度的函数,存放于zhengtigangdu[][]数组中//***********************************************void zonggang(){double gd[6][6],gd2[6][6]; // gd1为中间变量.double T[6][6],T1[6][6]; //T为坐标转换矩阵,T1为T的转置矩阵double danyuan_gangdu[6][6];int Ti,Ti2,Ti3;double chengji=0; //中间变量for(Ti=0;Ti<NJ*3;Ti++) //先将整体刚度付0{for(Ti2=0;Ti2<NJ*3;Ti2++)zhengtigangdu[Ti][Ti2]=0;}for(int i=0;i<NE;i++){dy_gangdu(i,danyuan_gangdu); //局部坐标系中的单元刚度js_T_T1( i, T,T1); //计算单元的转换矩阵T//以下为计算单元在整体坐标系中的单刚存在gd[]数组中.for(Ti=0;Ti<6;Ti++){for(Ti2=0;Ti2<6;Ti2++){gd[Ti][Ti2]=0;gd2[Ti][Ti2]=0;}}chengji=0;for(Ti=0;Ti<6;Ti++){for(Ti2=0;Ti2<6;Ti2++){for(Ti3=0;Ti3<6;Ti3++)chengji=chengji+T1[Ti][Ti3]*danyuan_gangdu[Ti3][Ti2]; //改过gd2[Ti][Ti2]=chengji;chengji=0;}chengji=0;for(Ti=0;Ti<6;Ti++){for(Ti2=0;Ti2<6;Ti2++){for(Ti3=0;Ti3<6;Ti3++)chengji=chengji+gd2[Ti][Ti3]*T[Ti3][Ti2]; //改过gd[Ti][Ti2]=chengji;chengji=0;}}chengji=0;//以下为将每个单元刚度集成到整体刚度中计算整体刚度int qianjiedian,houjiedian;qianjiedian=dym_jdm[i][0];houjiedian=dym_jdm[i][1];for(Ti=0;Ti<3;Ti++){for(Ti2=0;Ti2<3;Ti2++)zhengtigangdu[(qianjiedian-1)*3+Ti][(qianjiedian-1)*3+Ti2]=zhengtigangdu[(qianjiedian-1)*3+Ti][(qian jiedian-1)*3+Ti2]+gd[Ti][Ti2];}for(Ti=0;Ti<3;Ti++){for(Ti2=0;Ti2<3;Ti2++)zhengtigangdu[(houjiedian-1)*3+Ti][(houjiedian-1)*3+Ti2]=zhengtigangdu[(houjiedian-1)*3+Ti][(houjie dian-1)*3+Ti2]+gd[Ti+3][Ti2+3];}for(Ti=0;Ti<3;Ti++){for(Ti2=0;Ti2<3;Ti2++)zhengtigangdu[(qianjiedian-1)*3+Ti][(houjiedian-1)*3+Ti2]=zhengtigangdu[(qianjiedian-1)*3+Ti][(houji edian-1)*3+Ti2]+gd[Ti][Ti2+3];}for(Ti=0;Ti<3;Ti++){for(Ti2=0;Ti2<3;Ti2++)zhengtigangdu[(houjiedian-1)*3+Ti][(qianjiedian-1)*3+Ti2]=zhengtigangdu[(houjiedian-1)*3+Ti][(qianji edian-1)*3+Ti2]+gd[Ti+3][Ti2];}}}//**************************//输入数据的函数//**************************{char queding;//基本参数输入cout<<"请输入基本参数:单元数,节点数,支持数,有节点载荷的节点数,非节点载荷数;"<<endl<<endl;cin>>NE>>NJ>>NZ>>NPJ>>NPF;cout<<endl<<endl;//输入单元几何参数cout<<"请输入单元的几何参数:"<<endl;for(int i=0;i<NE;i++){cout<<"第"<<i+1<<"单元"<<endl;cout<<"长度:"; cin>>changdu[i];cout<<"角度:"; cin>>jiaodu[i]; jiaodu[i]=jiaodu[i]*PI/180;cout<<"弹性模量E:"; cin>>tanxingmoliang[i];cout<<"J模量:"; cin>>J_moliang[i];cout<<"面积:"; cin>>mianji[i];cout<<endl;//cout<<"第"<<i+1<<"单元长度:"<<changdu[i]<<" 角度:"<<jiaodu[i]<<" 弹性模量E:"<<tanxingmoliang[i]<<" J量:"<<J_moliang[i]<<" 面积:"<<mianji[i]<<endl<<endl;}//输入单元对应的节点总码cout<<"请输入个单元对应的节点总码"<<endl;for(i=0;i<NE;i++){cout<<"第"<<i+1<<"单元前面节点的总码:"; cin>>dym_jdm[i][0];cout<<"第"<<i+1<<"单元后面节点的总码:"; cin>>dym_jdm[i][1];}cout<<endl;//输入支承点位置cout<<"请输入支承点位置:"<<endl;for(i=0;i<NZ;i++){cout<<"第"<<i+1<<"支承位:"; cin>>zhichengweizhi[i];}cout<<endl;//输入节点载荷值cout<<"请输入节点载荷作用节点及其载荷值:"<<endl;for(i=0;i<NPJ;i++){cout<<"载荷所在节点位置:"; cin>>jdzhzuoyongdian[i];}for(i=0;i<NPJ;i++){cout<<"第"<<jdzhzuoyongdian[i]<<"节点的载荷值:"<<endl;cout<<"U:"; cin>>jiedianzaihe[i][0];cout<<"V:"; cin>>jiedianzaihe[i][1];cout<<"M:"; cin>>jiedianzaihe[i][2];}cout<<endl;cout<<"请输入非节点载荷作单元,载荷类型,载荷值,作用点(指离第一个节点的距离),作用角度:"<<endl;cout<<"其中1-均布,2-垂直集中,3-平行集中,4-力偶,5-角度集中"<<endl;for(i=0;i<NPF;i++){cout<<"第"<<i+1<<"个非节点载荷作用单元:"; cin>>fjzhzuoyongdanyuan[i];cout<<"载荷类型:"; cin>>fjzhleixing[i];cout<<"载荷值(如果是类型5的力,只输入力的绝对值) :"; cin>>fjzhzhi[i];cout<<"离第一个节点的距离:"; cin>>fjzhzuoyongdian[i];cout<<"力在该单元坐标中的作用角度:"; cin>>fjzhjiaodu[i];fjzhjiaodu[i]=fjzhjiaodu[i]*PI/180; cout<<endl;}cout<< endl<<endl;cout<<"按任意键显示输入情况:"<<endl;getch();//显示输入的情况cout<<"输入的数据如下:"<<endl;cout<<"===========================================:"<<endl;cout<<"单元数:"<<NE<<" 节点数:"<<NJ<<" 支持数:"<<NZ<<"有节点载荷的节点数:"<<NPJ<<" 非节点载荷数:"<<NPF<<endl;cout<<"============================================"<<endl;cout<<"各单元几何参数如下:"<<endl;cout<<"单元"<<"\t"<<"单元长度:"<<"\t"<<"角度: "<<"\t"<<"弹性模量E:"<<"\t"<<"J量: "<<"\t"<<"面积:"<<endl;for(i=0;i<NE;i++)cout<<i+1<<"\t"<<changdu[i]<<"\t\t"<<jiaodu[i]<<"\t\t"<<tanxingmoliang[i]<<"\t\t"<<J_moliang[i]<<"\ t\t"<<mianji[i]<<endl;cout<<"==========================================="<<endl;cout<<"各单元对应节点总码如下:"<<endl;cout<<"单元号"<<"\t"<<"前节点总码"<<"\t"<<"后节点总码"<<endl;for(i=0;i<NE;i++)cout<<i+1<<"\t\t"<<dym_jdm[i][0]<<"\t\t"<<dym_jdm[i][1]<<endl;cout<<" ==========================================="<<endl;cout<<"支承位置如下:"<<endl;cout<<"支承位置"<<"\t";for(i=0;i<NZ;i++)cout<<zhichengweizhi[i]<<" ";cout<<endl;cout<<" ==========================================="<<endl;cout<<"节点载荷输入如下:"<<endl;cout<<"节点位置"<<"\t"<<"U "<<"\t\t"<<"V "<<"\t\t"<<"M "<<endl;for(i=0;i<NPJ;i++){cout<<jdzhzuoyongdian[i]<<"\t\t"<<jiedianzaihe[i][0]<<"\t\t"<<jiedianzaihe[i][1]<<"\t\t"<<jiedianzaihe[ i][2];}cout<<" ==========================================="<<endl;cout<<"非节点载荷的作用情况如下:"<<endl;cout<<"作用单元"<<"\t"<<"载荷类型"<<"\t"<<"载荷值"<<"\t\t"<<"作用点"<<"\t\t"<<"作用角度"<<endl;for(i=0;i<NPF;i++)cout<<fjzhzuoyongdanyuan[i]<<"\t\t"<<fjzhleixing[i]<<"\t\t"<<fjzhzhi[i]<<"\t\t"<<fjzhzuoyongdian[i]<< "\t\t"<<fjzhjiaodu[i]<<endl;cout<<" ==========================================="<<endl;cout<<"按任意键继续!"<<endl;getch();cout<<endl;}///////////////////////////////////////////////////////程序结束////////////////////////////////////////////////////////////////// 用以上程序计算教材例题3-7以及3-8结果都和例题的结果一样。

ABAQUS 有限元分析大作业 凹槽成型实例

ABAQUS 有限元分析大作业 凹槽成型实例

有限元分析大作业-凹槽成型一:前处理-利用ABAQUS/CAE创建模型。

1、定义并创建四个部件如下图:图1:可变形毛坯Blank图2:刚性冲头Punch图3:刚性夹具Holder2、定义材料及截面特性图5:分别根据提供的数据定义Elastic和Plastic两个材料特性图6:创建均匀实体截面提交给材料Steel并赋予Blank此截面属性图7:建立一个局部数据坐标系(在随着毛坯运动的共旋坐标系下显示应力和应变的值): 2、装配部件图8:装配图(根据相互关系进行装配)3、创建几何集合创建6个几何集合:每个刚性体参考点各一个,毛坯对称面一个,毛坯中面的每段各一个。

图9:创建六个几何集合4、定义分析步和输出要求图10:创建分析步1 Establish contact 1图11:创建分析步2,Remove right constraint图12:创建分析步3 Holder force图13:创建分析步4 Establish contact 2图14:创建分析步5 Move punch图15:编辑场输出图16:编辑历史输出5、监控自由度的值图17:定义RefPunch的监控自由度6、定义接触相互作用图18:首先在Interaction模块中定义以上5个表面图19:定义一个无摩擦接触相互作用属性图20:再定义一个有摩擦相互作用属性,摩擦系数取为0.1 最后定义三个表面间的相互作用:图21:定义三个表面之间的相互作用(具体见下面三个图)图22:定义Die-Blank 相互作用 图23:定义Holder-Blank 相互作用图24:定义Punch-Blank 相互作用7、各分析步的边界条件首先在STEP1总添加每个初始边界条件:图25:Center、MidLeft、MidRight边界条件图26:RefDie、RefHolder、RefPunch边界条件图27:分析步1边界条件图28:Step2边界条件图29:Step3边界条件图30:Step4边界条件图31:Step5边界条件图:32集中力与负压力施加8、划分网格和定义作业图33:首先在毛坯上下与左右表面分别撒种图34:选择单元类型图35:划分网格二:求解1、创建作业图36:创建作业2、提交求解并监视求解过程图:37:求解完成图:38监视点U2向位移图三:后处理1、成型过程图39:开始图40:开始变形图41:变形扩大图42:最后成型结果2、绘制塑性应变等值线图3、绘制冲头上的反作用力4、绘制接触压力等值线图图45:接触压力等值线图。

有限元分析大作业

有限元分析大作业

一、有限元方法的手工计算结果与ansys分析结果的对比1分析的问题描述如图1所示,桁架的杆截面面积为8,由钢制成(E=200GPa)。

用有限元法计算出每个节点的位移以及反作用力。

(1)(2)(3)图1对于上述问题,本文将用手工计算和ansys软件分别计算出结果,对计算出来的结果进行对比。

2手工计算2.1桁架结构的有限元计算方法对于桁架结构,每个单元的刚度矩阵为,(2-1)YX图2其中,为桁架单元在整体坐标系中与X轴的夹角;,A为桁架的截面积,E 为弹性模量,L为桁架长度。

在固体力学问题中,有限元公式通常由如下的一般形式,Ku=F(2-2)其中,K为刚度矩阵,u为位移矩阵,F为载荷矩阵。

运用公式(2-3),就能求出反作用力,R=Ku-F(2-3)其中,R为反作用力矩阵。

2.2计算过程计算每个桁架单元的刚度,用公式(2-1)计算每个每个桁架单元的刚度矩阵,将每个单元放入总刚度矩阵,他们的位置分别为:10-100000 00000000 -10100000 00000000 00000000 00000000 00000000 0000000000000000 00000000 0010-1000 0000000000-101000 00000000 00000000 000000003.9-4.90000-3.9 4.9-4.9 6.10000 4.9-6.1 00000000 00000000 00000000 00000000-3.9 4.90000 3.9-4.9 4.9-6.10000-4.9 6.100000000 00000000 00000000 000 1.28000-1.28 00000000 00000000 00000000 000-1.28000 1.2800000000 00000000 00000000 00000000 0000 3.9 4.9-3.9-4.9 0000 4.9 6.1-4.9-6.1 0000-3.9-4.9 3.9 4.9 0000-4.9-6.1 4.9 6.1将个刚度矩阵相加得到总刚度矩阵为,19.9-4.90-16000-3.9 4.9 -4.9 6.100000 4.9-6.1 -160320-16000 00012.8000-12.8 00-16019.9 4.9-3.9-4.9 0000 4.9 6.1-4.9-6.1 -3.9 4.900-3.9-4.97.80 4.9-6.10-12.8-4.9-6.1025应用边界条件施加载荷,将总刚度矩阵带入式(2-2)得:19.9-4.90-16000-3.9 4.9Ux1 -4.9 6.100000 4.9-6.1Uy1 -160320-16000Ux2 00012.8000-12.8Uy200-16019.9 4.9-3.9-4.9Ux3 0000 4.9 6.1-4.9-6.1Uy3 -3.9 4.900-3.9-4.97.80Ux4 4.9-6.10-12.8-4.9-6.1025Uy4带入边界条件解得:将结果带入(2-3)得:=Fx1Fy1Fx2Fy2Fx3Fy3Fx4Fy43用ansys软件求解(单位统一N,mm,Mpa)(1)选择单元(图3)图3(2)附材料属性(图4)图4(3)创建模型(图5)图5(4)施加载荷(图6)图6(5)求解每个节点的位移(图7)图7节点的反力(图8)图8(6)模型变形图(7)位移等值线分布图4结果对比及分析手算结果ansys 计算结果位移(mm)Ux100Uy100Ux2-0.0016-0.0016Uy2-0.0468-0.0468Ux300Uy300Ux4-0.0066-0.0066Uy4-0.0317-0.0317表1手算结果ansys计算结果节点反力(N)Fx1-1027.8-1027.8 Fy11608.31608.3 Fx2 5.60 Fy2-100 Fx32066.72063.1 Fy32257.12255.4 Fx48.80 Fy4-2.70表2由表1和表2可以看出,手工计算的结果与ansys计算的结果基本一致。

国开最新《离散数学(本)》形考任务:大作业word版

国开最新《离散数学(本)》形考任务:大作业word版

离散数学大作业大作业时间为第1周到第17周,满分100分,由两部分组成。

提交作业方式有以下三种,请务必与辅导教师沟通后选择:1. 将此次作业用A4纸打印出来,手工书写答题,字迹工整,解答题要有解答过程,完成作业后交给辅导教师批阅。

注意选择此种提交方式时仍然需要在网络课提交作业入口处上传说明文档,文档内注明“作业已由线下提交给辅导老师”。

2. 在线提交word文档.3. 自备答题纸张,将答题过程手工书写,并拍照上传.第一部分一、公式翻译题(每小题2分,共10分)1.将语句“我会英语,并且会德语.”翻译成命题公式.设p.我学英语Q:我学法语则命题公式为:pΛQ2.将语句“如果今天是周三,则昨天是周二.”翻译成命题公式.设P:今天是周三Q:昨天是周二则命题公式为:P→Q3.将语句“小王是个学生,小李是个职员.”翻译成命题公式.设P:小王是个学生Q:小李是个职员则命题公式为:P∧Q4.将语句“如果明天下雨,我们就去图书馆.”翻译成命题公式.设 P 表示“明天下雨”Q 表示“我们就去图书馆”命题公式:P → Q5.将语句“当大家都进入教室后,讨论会开始进行.”翻译成命题公式.设 P :大家都进入教室后Q :讨论会开始进行命题公式:P → Q二、计算题(每小题10分,共50分)1.设集合A={1, 2, 3},B={2, 3, 4},C={2, {3}},试计算(1)A-C;(2)A∩B;(3)(A∩B)×C.答:(1)A-C ={1,3};(2)A∩B={2,3};(3)(A∩B)×C={<2,2>,<2,{3}>,<3,2>,<3,{3}>}。

2. 设G =<V ,E >,V ={v 1, v 2, v 3, v 4, v 5},E ={(v 1,v 3) , (v 1,v 5) , (v 2,v 3) , (v 3,v 4) , (v 4,v 5) },试(1)给出G 的图形表示;(2)求出每个结点的度数;(3)画出其补图的图形.答:(1)G 的图形表示如图所示(2)v1, v2, v3, v4, v5结点的度数依次为2,1,3,2,2(3)补图的图形3.试画一棵带权为1, 2, 3, 3, 4的最优二叉树,并计算该最优二叉树的权.最优二叉树的权为1×3+2×3+3×2+3×2+4×2=294.求出如下所示赋权图中的最小生成树(要求写出求解步骤),并求此最小生成树的权.W(v2,v6)=1,选(v2,v6)W(v4,v5)=1,选(v4,v5)W(v1,v6)=2,选(v1,v6)W(v3,v5)=2,选(v3,v5)W(v2,v3)=4,选(v2,v3)最小生成树,如图生成树的权W(T)=1+1+2+2+4=10 ο ο ο ο ο v 6 v 1 v 2 v 5 v 3 ο v 4 1 6 2 4 5 7 9 3 1 5 2 ο ο ο ο ο v 6 v 1 v 2 v 5 v 3ο v 4 1 6 2 4 57 9 3 1 5 25.求P→(Q∧R) 的析取范式与合取范式.P→(Q∧R)=┐P∨(Q∧R)=(┐P∨Q)∧(┐P∨R)合取范式=(┐P∨Q)∨(R∧┐R)∧(┐P∨R)=(┐P∨Q)∨(R∧┐R)∧(┐P∨R)∨(Q∧┐Q)=(┐P∨Q∨R)∧(┐P∨Q∨┐R)∧(┐P∨┐Q∨R)主合取范式=(┐P∧┐Q∧┐R)∨(┐P∧┐Q∧R)∨(┐P∧┐Q∧┐R)(┐P∧Q∧R)∨(P∧┐Q∧R)∨(P∧Q∧┐R)∨(P∧Q∧R)主析取范式第二部分从下列选题中选择一个感兴趣的主题,自主查阅文献资料进行深入的研究和学习,并形成一份至少一千字的总结报告。

《结构分析中的有限元法》2015-有限元习题-参考答案

《结构分析中的有限元法》2015-有限元习题-参考答案
5、简述有限元法的发展和现状。
近几十年,伴随着计算机科学和技术的快速发展,有限元法作为工程分析的 有效方法在理论、方法的研究、计算机程序的开发以及应用领域的开拓者方面均 取得了根本性的发展。
(1)单元的类型和形式 为了扩大有限元法的应用领域,新的单元类型和形式不断涌现(等参元,梁板 壳,复合材料) (2)有限元法的理论基础和离散格式 将 Hellinger-Reissner、Hu—Washizu(多场变量变分原理)应用于有限元分析, 发展了混合模型、杂交型的有限元表达格式,应研究了各自的收敛条件;将加权 余量法用于建立有限元的表达格式;进一步研究发展有限元解的后验误差估计和 应力磨平方法。 (3)有限元方程的解法(大型复杂工程结构问题——静态, 特征值, 瞬态等) (4)有限元法的计算机软件(专用软件, 通用软件)
4、说明用有限单元法解题的主要步骤。 答:研究问题的力学建模;结构离散;单元分析;整体分析与求解;结果分析及 后处理。
5、推导基于变分原理的总势能泛函极值条件。 解:有积分形式确立的标量泛函有
Π
F
u,
u x
,

E F 和 E 是特定的算子, 是求解域, 是 的边界。 Π 称 为未知函数 u 的泛函,随函数 u 的变化而变化。连续介质问题的解 u 使泛函 Π 对 于微小的变化u 取驻值,即泛函的“变分”等于零 Π 0 ,此为变分法。
物理意义:应力分量与体力分量之间的关系。 (2)几何方程:
x
u x
, y
v y
,z
w z
xy
u y
v x
,
yz
v z
w y
,
zx
w x
u z
物理意义:应变分量与位移分量之间的关系。 (3)物理方程:

有限元分析大作业程序部分

有限元分析大作业程序部分
>> L2=0.96;
>> k1=Beam2D2Node_Stiffness(E,I,A,L1);
>> k2=Beam2D2Node_Stiffness(E,I,A,L2);
7 0 0 -1.4167 0 0
0 0.0078 0.0056 0 -0.0078 0.0056
3,建立整体刚度方程:
>> T=[0,1,0,0,0,0;-1,0,0,0,0,0;0,0,1,0,0,0;0,0,0,0,1,0;0,0,0,-1,0,0;0,0,0,0,0,1];
>> k3=T'*k2*T;
>> KK=zeros(12,12);
>> KK=Beam2D2Node_Assemble(KK,k1,1,2);
functionkbeam2d2nodestiffnesseial2矩阵组装functionbeam2d2nodeassemblekkkij以上函数进行单元刚度矩阵的组装输入单元刚度矩阵k单元的节点编号ij输出整体刚度矩阵kkdof13i2
程序部分
1,计算单元刚度矩阵(输入E,A,I,L)然后输出k(6*6)刚度矩阵:
functionk=Beam2D2Node_Stiffness(E,I,A,L)
k=[E*A/L,0,0,-E*A/L,0,0;
0,12*E*I/(L^3),6*E*I/(L^2),0,-12*E*I/(L^3),6*E*I/(L^2);
0,6*E*I/(L^2),4*E*I/L,0,-6*E*I/(L^2),2*E*I/L;
DOF(6)=3*j;
for n1=1:6
for n2=1:6
KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);

ANSYS大作业扳手有限元分析

ANSYS大作业扳手有限元分析

弹性力学与有限元法扳手有限元分析[摘要]运用 ANSYS 软件中参数化语言设计与蒙特卡罗法相结合的随机有限元法,以扳手为研究对象,选择扳手的多种设计参数和作用载荷为随机变量,在假定的统计分布情况下进行扳手的强度可靠性分析。

运用增量求解的非线性有限元分析方法,采用扩展的拉格朗日算法和库仑摩擦模型计算扳手中的摩擦接触问题,得到了扳手应力应变精确的数值解,通过分析发现该产品应力应变不均布,存在应力集中区域,而其余不少区域应力和应变较小,材料富余。

应用UG NX对扳手进行三维建模,基于ANSYS有限元分析平台下显示静力学模块(Static Structural)建立扳手静力学模型,对扳手进行应力分析和变形量分析。

经分析计算,找出设计结构齿根断裂破坏原因,从而对原结构进行改进。

[1]关键词:扳手;有限元法;ANSYS;引言扳手是一种常用的安装与拆卸工具,是利用杠杆原理拧转螺栓、螺钉、螺母和其他螺纹紧持螺栓或螺母的开口或套孔固件的手工工具。

扳手通常用碳素或合金材料的结构钢制造。

蒙特卡罗法亦称统计模拟方法,是 Metropolis 在第二次世界大战期间提出的用于研究与原子弹有关的中子输运过程的一种方法。

该方法目前已经广泛应用于多种学科对随机过程的统计模拟,如计量过程中不确定度的评定和分析、高等数学中数值积分和规划问题的求解、实验探测器的模拟以及数值分析、工程结构的可靠性及灵敏度分析等。

[2]采用蒙特卡罗法对大量的随机过程进行建模以及仿真模拟分析,可以使复杂的随机问题得到解决。

特别是计算机的出现,通过计算机软件使大量繁琐运算通过计算机得以完成,大大提高了蒙特卡罗方法的应用范围和工作效率。

[3]1 扳手模型建立在UG NX中建立扳手模型,UG NX版本太高,ANSYS不能直接识别其文件,需要另存为Step格式。

图1 在UG NX 中绘制的扳手模型2 在ANSYS Workbench 中建立工程文件2.1 打开Static Structural在Geometry 中导入Step 文件,打开Model 选项。

有限元分析与应用大作业

有限元分析与应用大作业

有限元分析及应用大作业课程名称: 有限元分析及应用班级:姓名:试题2:图示薄板左边固定,右边受均布压力P=100Kn/m作用,板厚度为0.3cm;试采用如下方案,对其进行有限元分析,并对结果进行比较。

1)三节点常应变单元;(2个和200个单元)2)四节点矩形单元;(1个和50个单元)3)八节点等参单元。

(1个和20个单元)图2-1 薄板结构及受力图一、建模由图2-1可知,此薄板长和宽分别为2m和1.5m,厚度仅为0.3cm,本题所研究问题为平面应力问题。

经计算,平板右边受均匀载荷P=33.33MPa,而左边被固定,所以要完全约束个方向的自由度,如图2-2所示。

取弹性模量E=2.1×11Pa,泊松比μ=0.3。

P=33.33MPa图2-2 数学模型二、第一问三节点常应变单元(2个和200个单元)三节点单元类型为PLANE42,设置好单元类型后,实常数设置板厚为0.3M。

采用2个单元的网格划分后的结果如图2-3,200个单元的网格划分图如图2-6所示。

约束的施加方式和载荷分布如图2-2中所示。

约束右边线上节点全部自由度。

计算得到的位移云图分别如图2-4、7所示,应力云图如图2-5、8所示。

图2-3 2个三角形单元的网格划分图图2-4 2个三角形单元的位移云图图2-5 2个三角形单元的应力云图图2-6 200个三角形单元的网格划分图图2-7 200个三角形单元的位移云图图2-8 200个三角形单元的应力云图三、第二问四节点矩形单元的计算四节点单元类型为PLANE42,设置好单元类型后,实常数设置板厚为0.3M。

采用1个单元的网格划分后的结果如图2-9,50个单元的网格划分图如图2-12所示。

约束的施加方式和载荷分布如图2-2中所示。

约束右边线上节点全部自由度。

计算得到的位移云图分别如图2-10、11所示,应力云图如图2-13、14所示。

图2-9 1个四边形单元的网格划分图图2-10 1个四边形单元的位移云图图2-11 1个四边形单元的应力云图图2-12 50个四边形单元的网格划分图图2-13 50个四边形单元的位移云图图2-14 50个四边形单元的应力云图四、第三问八节点等参单元的计算四节点单元类型为PLANE82,设置好单元类型后,实常数设置板厚为0.3M。

ANSYS大作业_轴承座有限元分析

ANSYS大作业_轴承座有限元分析

轴承座轴瓦 轴四个安装孔径向约束 (对称) 轴承座底部约束 (UY=0)沉孔上的推力 (3000 psi.) 向下作用力 (15000 psi.) 基于ANSYS 的轴承座有限元分析一、 问题描述在我们机械设计课程中曾经学习过轴系,主要是学习了轴的设计、受力分析以及轴承的设计等等。

但没有对轴承座的承受能力进行分析,所以我在这里主要是对一种简单的轴承座进行了有限元分析。

在查阅了相关资料之后,可将分析的轴承座示意如下图。

在实际当中,考虑到工艺的要求,图中相应的边缘处须设置有圆角、倒边等等。

但在有限元模型中忽略了这些要素。

二、 力学模型的分析与建立如下图所示在查阅了相关资料后可将上面描述的问题简化成上述模型,其中的载荷参考了网上的相关资料,在沉孔面上垂直于沉孔面上作用有3000psi.的推力载荷,在轴承孔的下半部分施加15000psi.的径向压力载荷,这个载荷是由于受重载的轴承受到支撑作用而产生的。

由于轴承座一般固定于机身上,所以可以在其底部施加法向位移约束,并且四个安装孔要受到螺栓的约束,所以可以在四个螺栓孔中施加径向对称约束(在ansys中体现为Symmetry B.C.)三、力学模型的有限元分析1.建立模型1)创建基座模型生成长方体Main Menu:Preprocessor->Modeling->Create->Volumes->Block->By Dimensions输入x1=0,x2=3,y1=0,y2=1,z1=0,z2=3平移并旋转工作平面Utility Menu>WorkPlane->Offset WP by IncrementsX,Y,Z Offsets 输入2.25,1.25,.75 点击ApplyXY,YZ,ZX Angles输入0,-90点击OK。

创建圆柱体Main Menu:Preprocessor->Modeling->Create->Volumes->Cylinder> Solid CylinderRadius输入0.75/2, Depth输入-1.5,点击OK。

最新《有限单元法》复习参考题

最新《有限单元法》复习参考题

精品资料《有限单元法》复习参考题........................................《有限单元法》复习参考题一、简答题:1、简述应用有限单元法解决具体问题的要点。

(1) 将一个表示结构或者连续体的求解域离散为若干个子域(单元),并通过他们边界上的结点相互结合为组合体。

(2) 用每个单元内所假设的近似函数来分片地表示全求解域内待求的未知场变量。

而每个单元内的近似函数由未知场函数(或及其导数,为了叙述方便,后面略去此加注)在单元各个节点上的数值与其对应的插值函数来表达。

(3) 通过和原问题数学模型(基本方程、边界条件)等效的变分原理或者加权余量法,建立求解基本未知量(场函数的结点值)的代数方程或者常微分方程组。

2、等效积分形式和等效积分“弱”形式的区别何在?为什么等效积分“弱”形式在数值分析中得到更多的应用?在很多情况下对微分方程的等效积分形式进行分部积分可以得到等效积分的弱形式,如下式T T C D E ()F()d 0ΩΓυΩ+υυΓ=⎰⎰()(u)d ,其中C 、D 、E 、F 是微分算子。

像这种通过适当提高对任意函数和υ 的连续性要求,以降低对微分方程场函数u 的连续性要求所建立的等效积分形式称为微分方程的等效积分“弱”形式。

值得指出的是,从形式上看“弱”形式对函数u 的连续性要求降低了,但对于实际的物理问题却常常较原始的微分方程更逼近真正的解,因为原始微分方程往往对解提出了过分的要求。

所以等效积分“弱”形式在数值分析中得到更多的应用。

3、什么是Ritz (里兹)方法?其优缺点是什么?收敛的条件是什么?基于变分原理的近似解法称为Ritz (里兹),解法如下:优缺点:一般来说,使用里兹方法求解,当试探函数族的范围扩大以及待定参数的数目增多时,近似解的精度将会提高。

局限性:(1) 在求解域比较复杂的情况下,选取满足边界条件的试探函数,往往会产生难以克服的困难。

(2) 为了提高近似解的精度,需要增加待定参数,即增加试探函数的项数,这就增加了求解的复杂性,而且由于试探函数定义于全域,因此不可能根据问题的要求在求解域的不同部位对试探函数提出不同精度的要求,往往由于局部精度的要求使整个问题求解增加许多困难。

有限元考试试题

有限元考试试题

有限元考试试题一、选择题(每题5分,共30分)1、在有限元分析中,我们通常使用什么方法来求解偏微分方程?A.积分法B.差分法C.有限差分法D.有限元法2、下列哪个不是有限元法的优点?A.可以处理复杂几何形状B.可以处理非线性问题C.可以处理大规模问题D.可以处理不稳定问题3、在有限元分析中,我们通常将连续的物理场离散化为一系列的什么?A.有限个点B.无限个小段C.有限个小段D.无限个点4、下列哪个不是有限元分析的基本步骤?A.划分网格B.建立模型C.执行计算D.编写代码5、在有限元分析中,我们通常使用什么来描述物理场的性质?A.偏微分方程B.泛函方程C.常微分方程D.边界条件6、下列哪个不是有限元分析的应用领域?A.结构分析B.流体动力学C.电磁学D.社会科学二、填空题(每题10分,共40分)7、______是一种将连续的物理场离散化为一系列有限个点的方法,是有限元分析的基础。

8、在有限元分析中,我们通常使用______来对物理场进行离散化处理。

9、______是一种求解偏微分方程的数值方法,广泛应用于有限元分析。

10、在有限元分析中,我们通常使用______来描述物理场的性质。

三、解答题(每题20分,共60分)11、请简述有限元分析的基本步骤,并解释其在结构分析中的应用。

12、请说明在有限元分析中,如何处理边界条件,并举例说明。

13、请简述有限元分析的优点和局限性。

有限空间培训考试试题及答案一、选择题1、在有限空间内,以下哪个行为是危险的?A.带压操作B.穿著宽松衣服C.使用电动工具D.所有上述答案:D.所有上述。

在有限空间内,带压操作、穿著宽松衣服和使用电动工具都是危险的。

2、当进入有限空间前,应该进行哪项操作?A.排放内部气体B.测试内部气体C.对内部进行冲洗D.所有上述答案:D.所有上述。

在进入有限空间前,应该进行排放内部气体、测试内部气体并对内部进行冲洗。

3、有限空间内的危险因素不包括以下哪个?A.缺氧B.有毒气体C.电击D.所有上述答案:C.电击。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

有限元法大作业一平面刚架的程序用Visual C++编制的平面刚架的源程序如下:///////////////////////////////////////////////////////程序开始//////////////////////////////////////////////////////////////////#include"iostream.h"#include"math.h"#include"stdlib.h"#include"conio.h"//*****************//声明必要变量//*****************#define PI 3.141592654int NE; //单元数int NJ; //节点数int NZ; //支承数int NPJ; //有节点载荷作用的节点数int NPF; //非节点载荷数int HZ; //载荷码int E; //单元码int fangchengshu;double F[303]; //各节点等效总载荷数值int dym_jdm[100][2]; //单元码对应的节点码:dym_jdm[][0], dym_jdm[][1]分为前后节点总码int zhichengweizhi[300]; //记录支持节点作用点的数组int fjzhzuoyongdanyuan[100]; //非节点载荷作用单元int fjzhleixing[100]; //非节点载荷类型:1-均布,2-垂直集中,3-平行集中,4-力偶,5-角度集中double fjzhzhi[100]; //非节点载荷的值double fjzhzuoyongdian[100]; //非节点载荷在各竿的作用点double fjzhjiaodu[100]; //非节点载荷作用角度int jdzhzuoyongdian[100]; //节点载荷作用的节点数组double jiedianzaihe[101][3];//节点载荷值,其jiedianzaihe[][0]-- jiedianzaihe[][2]分别为U, V, Mdouble zhengtigangdu[303][303]; //整体刚度数组double changdu[100]; //各单元竿长数组double jiaodu[100]; //各单元角度数组double tanxingmoliang[100]; //各单元弹性模量数组double J_moliang[100]; //各单元J模量数组double mianji[100]; //各单元面积数组double weiyi[303]; //记录各个节点位移的数组double dy_weiyi[100][6]; //各个单元在局部坐标系中的位移数组dy_weiyi[i][0]-dyweiyi[i][6]分别为第i+1单元的u1,v1,@1,u2,v2,@2double dy_neili[100][6];//各个单元在局部坐标系中的固端内力dy_weiyi[i][0]-dyweiyi[i][6]分别为第i+1单元的U1,V1,M1,U2,V2,M2double gan_neili[100][6];//各个单元的竿端内力数组,gan_neili[i][6]表示第i+1单元的6内力.//*******************//*******************void input(); //数据的输入void zonggang(); //计算总刚度,存放于zhengtigangdu[][]数组中void zongzaihe(); //计算等效总节点载荷void zhichengyinru(); //引入支承条件void jsweiyi(); //求各个节点位移void js_dy_weiyi(); //求局部坐标系中的位移void ganduanneili(); //求竿端内力void dy_gangdu(int i,double dg[6][6]); //求单元在局部坐标系中的单刚void js_T_T1(int i,double T[6][6],double T1[6][6]); //求单元的转换矩阵及其逆阵//************//主函数//************void main(){input();cout<<" 输出结果"<<endl;cout<<"==============================”<<endl;zonggang();zongzaihe();zhichengyinru();jsweiyi();cout<<"输出位移:"<<endl;for( int i=0;i<NJ;i++){cout<<"u["<<i+1<<"]:"<<weiyi[i*3]<<endl;cout<<"v["<<i+1<<"]:"<<weiyi[i*3+1]<<endl;cout<<"@["<<i+1<<"]:"<<weiyi[i*3+2]<<endl;}js_dy_weiyi();ganduanneili();//输出竿端内力cout<<"==============================”<<endl;cout<<"输出竿内力:"<<endl;for( i=0;i<NE;i++){cout<<"第"<<i+1<<"单元:"<<endl;cout<<" U1:"<<gan_neili[i][0]<<" V1:"<<gan_neili[i][01]<<" M1:"<<gan_neili[i][2]<<endl;cout<<" U2:"<<gan_neili[i][3]<<" V2:"<<gan_neili[i][4]<<" M2:"<<gan_neili[i][5]<<endl;}cout<<" ==========================================="<<endl;cout<<"输出完毕,按任意键退出!"<<endl;getchi();}//***************************//***************************void dy_gangdu(int i,double dg[6][6]){dg[0][0]=tanxingmoliang[i]*mianji[i]/changdu[i]; dg[0][1]=0; dg[0][2]=0;dg[0][3]=-tanxingmoliang[i]*mianji[i]/changdu[i]; dg[0][4]=0; dg[0][5]=0; dg[1][0]=0;dg[1][1]=12*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]*changdu[i]);dg[1][2]=-6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]); dg[1][3]=0;dg[1][4]=-12*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]*changdu[i]);dg[1][5]=-6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]); dg[2][0]=0;dg[2][1]=-6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]);dg[2][2]=4*tanxingmoliang[i]*J_moliang[i]/changdu[i]; dg[2][3]=0;dg[2][4]=6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]);dg[2][5]=2*tanxingmoliang[i]*J_moliang[i]/changdu[i];dg[3][0]=-tanxingmoliang[i]*mianji[i]/changdu[i]; dg[3][1]=0; dg[3][2]=0;dg[3][3]=tanxingmoliang[i]*mianji[i]/changdu[i]; dg[3][4]=0; dg[3][5]=0; dg[4][0]=0;dg[4][1]=-12*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]*changdu[i]);dg[4][2]=6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]); dg[4][3]=0;dg[4][4]=12*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]*changdu[i]);dg[4][5]=6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]); dg[5][0]=0;dg[5][1]=-6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]);dg[5][2]=2*tanxingmoliang[i]*J_moliang[i]/changdu[i]; dg[5][3]=0;dg[5][4]=6*tanxingmoliang[i]*J_moliang[i]/(changdu[i]*changdu[i]);dg[5][5]=4*tanxingmoliang[i]*J_moliang[i]/changdu[i];}//*************************//求单元的转换矩阵及其逆阵//*************************void js_T_T1(int i,double T[6][6],double T1[6][6]){int Tti,Tti2;for( Tti=0;Tti<6;Tti++){//计算单元的转换矩阵Tfor( Tti2=0;Tti2<6;Tti2++){if(Tti==0&&Tti2==0)T[Tti][Tti2]=cos(jiaodu[i]);else if(Tti==0&&Tti2==1)T[Tti][Tti2]=sin(jiaodu[i]);else if(Tti==1&&Tti2==0)T[Tti][Tti2]=-sin(jiaodu[i]);else if(Tti==1&&Tti2==1)T[Tti][Tti2]=cos(jiaodu[i]);else if(Tti==2&&Tti2==2)T[Tti][Tti2]=1;T[Tti][Tti2]=cos(jiaodu[i]);else if(Tti==3&&Tti2==4)T[Tti][Tti2]=sin(jiaodu[i]);else if(Tti==4&&Tti2==3)T[Tti][Tti2]=-sin(jiaodu[i]);else if(Tti==4&&Tti2==4)T[Tti][Tti2]=cos(jiaodu[i]);else if(Tti==5&&Tti2==5)T[Tti][Tti2]=1;elseT[Tti][Tti2]=0;if(fabs(T[Tti][Tti2])<0.0000001)T[Tti][Tti2]=0;if(T[Tti][Tti2]>0.99999999)T[Tti][Tti2]=1;if(T[Tti][Tti2]<-0.99999999)T[Tti][Tti2]=-1;}}//计算T的转置矩阵T1for(Tti=0;Tti<6;Tti++){for(Tti2=0;Tti2<6;Tti2++)T1[Tti2][Tti]=T[Tti][Tti2];}}//*******************//计算竿端内力的函数//*******************void ganduanneili(){int Ti,Ti2;double chengji=0;double danyuan_gangdu[6][6]; //单元刚度数组for(int i=0;i<NE;i++){dy_gangdu(i,danyuan_gangdu);for(Ti=0;Ti<6;Ti++){for(Ti2=0;Ti2<6;Ti2++)chengji=chengji+danyuan_gangdu[Ti][Ti2]*dy_weiyi[i][Ti2];gan_neili[i][Ti]=chengji+dy_neili[i][Ti];chengji=0;}}//************************************************************//计算各个单元在局部坐标系中的位移的函数,存放在dy_weiyi[][]中//************************************************************void js_dy_weiyi(){int Ti,Ti3;double T[6][6],T1[6][6];double chengji=0;for(int i=0;i<NE;i++){js_T_T1(i,T,T1);int qianjiedian,houjiedian;qianjiedian=dym_jdm[i][0];houjiedian=dym_jdm[i][1];for(Ti=0;Ti<6;Ti++){dy_weiyi[i][Ti]=T[Ti][0]*weiyi[(qianjiedian-1)*3]+T[Ti][1]*weiyi[(qianjiedian-1)*3+1]+T[Ti][2]* weiyi[(qianjiedian-1)*3+2]+T[Ti][3]*weiyi[(houjiedian-1)*3]+T[Ti][4]*weiyi[(houjiedian-1)*3+1]+T[Ti][5]* weiyi[(houjiedian-1)*3+2];}}}//***********************//计算得位移的函数//***********************void jsweiyi(){int Ti,Ti2;int i,j,k;int i_max;double max,zhongjian;double gd_kz[303][303];for(Ti=0;Ti<NJ*3;Ti++){for(Ti2=0;Ti2<NJ*3;Ti2++)gd_kz[Ti][Ti2]=zhengtigangdu[Ti][Ti2];gd_kz[Ti][NJ*3]=F[Ti];}//用高斯主列消去法解方程,先计算增广矩阵的变形矩阵for(Ti=0;Ti<NJ*3-1;Ti++){max=fabs(gd_kz[Ti][Ti]);i_max=Ti;for(Ti2=Ti;Ti2<NJ*3;Ti2++){{i_max=Ti2-1;max=gd_kz[Ti2][Ti];}}if(i_max!=Ti){for(i=Ti;i<NJ*3+1;i++){zhongjian=gd_kz[Ti][i];gd_kz[Ti][i]=gd_kz[i_max][i];gd_kz[i_max][i]=zhongjian;}}for(j=Ti+1;j<NJ*3;j++){double shang=gd_kz[j][Ti]/gd_kz[Ti][Ti];for(k=0;k<NJ*3+1;k++)gd_kz[j][k]=-shang*gd_kz[Ti][k]+gd_kz[j][k];}}//反代求值double chengji=0;weiyi[NJ*3-1]=gd_kz[NJ*3-1][NJ*3]/gd_kz[NJ*3-1][NJ*3-1];for(i=NJ*3-2;i>=0;i--){for(j=i+1;j<NJ*3;j++)chengji=chengji+gd_kz[i][j]*weiyi[j];weiyi[i]=(gd_kz[i][NJ*3]-chengji)/gd_kz[i][i];chengji=0;}}//********************//引入支承条件的函数//********************void zhichengyinru(){int i,Ti;int weizhi;fangchengshu=NJ*3-NZ;for(i=0;i<NZ;i++){weizhi=zhichengweizhi[i]-1;for(Ti=0;Ti<NJ*3;Ti++)zhengtigangdu[weizhi][Ti]=0;zhengtigangdu[Ti][weizhi]=0;zhengtigangdu[weizhi][weizhi]=1;F[weizhi]=0;}}//*************************//计算等效总节点载荷的函数//*************************void zongzaihe(){int i,Ti;double U1=0,V1=0,M1=0,U2=0,V2=0,M2=0;double c,G,l,d;double Ff[6];double T[6][6],T1[6][6]; //T为坐标转换矩阵,T1为T的转置矩阵for(i=0;i<NJ*3;i++)F[i]=0;for(i=0;i<6;i++)Ff[i]=0;for( i=0;i<NE;i++){for(int jl=0;jl<6;jl++)dy_weiyi[i][jl]=0;}for(i=0;i<NPF;i++){int dym=fjzhzuoyongdanyuan[i];U1=V1=M1=U2=V2=M2=0;for(int t=0;t<6;t++)Ff[t]=0;c=fjzhzuoyongdian[i];l=changdu[(fjzhzuoyongdanyuan[i]-1)];d=l-c;G=fjzhzhi[i];if(fjzhleixing[i]==1){V1=-0.5*G*c*(2-2*c*c/(l*l)+c*c*c/(l*l*l));V2=-G*c-V1;M1=G*c*c*(6-8*c/l+3*c*c/(l*l))/12;M2=-G*c*c*c*(4-3*c/l)/(12*l);U1=U2=0;}else if(fjzhleixing[i]==2){V1=-G*(l+2*c)*d*d/(l*l*l);M1=G*c*d*d/(l*l);M2=-G*c*c*d/(l*l);U1=U2=0;}else if(fjzhleixing[i]==3){U1=-G*d/l;U2=-G*c/l;M1=M2=V1=V2=0;}else if(fjzhleixing[i]==4){U1=U2=0;V1=-6*G*c*d/(l*l*l);V2=-V1;M1=G*d*(2-3*d/l)/l;M2=G*c*(2-3*c/l)/l;}else if(fjzhleixing[i]==5){V1=-G*sin(fjzhjiaodu[i])*(l+2*c)*d*d/(l*l*l);V2=-G*sin(fjzhjiaodu[i])*(l+2*d)*c*c/(l*l*l);M1=G*sin(fjzhjiaodu[i])*c*d*d/(l*l);M2=-G*sin(fjzhjiaodu[i])*c*c*d/(l*l);U1=U2=0;U1=-G*cos(fjzhjiaodu[i])*d/l+U1;U2=-G*cos(fjzhjiaodu[i])*c/l+U2;}//记录竿的固端内力dy_neili[dym-1][0]=U1;dy_neili[dym-1][1]=V1;dy_neili[dym-1][2]=M1;dy_neili[dym-1][3]=U2;dy_neili[dym-1][4]=V2;dy_neili[dym-1][5]=M2;js_T_T1( fjzhzuoyongdanyuan[i]-1, T,T1);for(Ti=0;Ti<6;Ti++)Ff[Ti]=-(T1[Ti][0]*U1+T1[Ti][1]*V1+T1[Ti][2]*M1+T1[Ti][3]*U2+T1[Ti][4]*V2+T1[Ti][5]*M2);//将载荷转换到整体坐标系中F[(dym_jdm[dym-1][0]-1)*3]=F[(dym_jdm[dym-1][0]-1)*3]+Ff[0];F[(dym_jdm[dym-1][0]-1)*3+1]=F[(dym_jdm[dym-1][0]-1)*3+1]+Ff[1];F[(dym_jdm[dym-1][0]-1)*3+2]=F[(dym_jdm[dym-1][0]-1)*3+2]+Ff[2];F[(dym_jdm[dym-1][1]-1)*3]=F[(dym_jdm[dym-1][1]-1)*3]+Ff[3];F[(dym_jdm[dym-1][1]-1)*3+1]=F[(dym_jdm[dym-1][1]-1)*3+1]+Ff[4];F[(dym_jdm[dym-1][1]-1)*3+2]=F[(dym_jdm[dym-1][1]-1)*3+2]+Ff[5];for(i=0;i<NPJ;i++){F[(jdzhzuoyongdian[i]-1)*3]=F[(jdzhzuoyongdian[i]-1)*3]+jiedianzaihe[i][0];F[(jdzhzuoyongdian[i]-1)*3+1]=F[(jdzhzuoyongdian[i]-1)*3+1]+jiedianzaihe[i][1];F[(jdzhzuoyongdian[i]-1)*3+2]=F[(jdzhzuoyongdian[i]-1)*3+2]+jiedianzaihe[i][2];}}//***********************************************//计算总刚度的函数,存放于zhengtigangdu[][]数组中//***********************************************void zonggang(){double gd[6][6],gd2[6][6]; // gd1为中间变量.double T[6][6],T1[6][6]; //T为坐标转换矩阵,T1为T的转置矩阵double danyuan_gangdu[6][6];int Ti,Ti2,Ti3;double chengji=0; //中间变量for(Ti=0;Ti<NJ*3;Ti++) //先将整体刚度付0{for(Ti2=0;Ti2<NJ*3;Ti2++)zhengtigangdu[Ti][Ti2]=0;}for(int i=0;i<NE;i++){dy_gangdu(i,danyuan_gangdu); //局部坐标系中的单元刚度js_T_T1( i, T,T1); //计算单元的转换矩阵T//以下为计算单元在整体坐标系中的单刚存在gd[]数组中.for(Ti=0;Ti<6;Ti++){for(Ti2=0;Ti2<6;Ti2++){gd[Ti][Ti2]=0;gd2[Ti][Ti2]=0;}}chengji=0;for(Ti=0;Ti<6;Ti++){for(Ti2=0;Ti2<6;Ti2++){for(Ti3=0;Ti3<6;Ti3++)chengji=chengji+T1[Ti][Ti3]*danyuan_gangdu[Ti3][Ti2]; //改过gd2[Ti][Ti2]=chengji;chengji=0;}chengji=0;for(Ti=0;Ti<6;Ti++){for(Ti2=0;Ti2<6;Ti2++){for(Ti3=0;Ti3<6;Ti3++)chengji=chengji+gd2[Ti][Ti3]*T[Ti3][Ti2]; //改过gd[Ti][Ti2]=chengji;chengji=0;}}chengji=0;//以下为将每个单元刚度集成到整体刚度中计算整体刚度int qianjiedian,houjiedian;qianjiedian=dym_jdm[i][0];houjiedian=dym_jdm[i][1];for(Ti=0;Ti<3;Ti++){for(Ti2=0;Ti2<3;Ti2++)zhengtigangdu[(qianjiedian-1)*3+Ti][(qianjiedian-1)*3+Ti2]=zhengtigangdu[(qianjiedian-1)*3+Ti][(qian jiedian-1)*3+Ti2]+gd[Ti][Ti2];}for(Ti=0;Ti<3;Ti++){for(Ti2=0;Ti2<3;Ti2++)zhengtigangdu[(houjiedian-1)*3+Ti][(houjiedian-1)*3+Ti2]=zhengtigangdu[(houjiedian-1)*3+Ti][(houjie dian-1)*3+Ti2]+gd[Ti+3][Ti2+3];}for(Ti=0;Ti<3;Ti++){for(Ti2=0;Ti2<3;Ti2++)zhengtigangdu[(qianjiedian-1)*3+Ti][(houjiedian-1)*3+Ti2]=zhengtigangdu[(qianjiedian-1)*3+Ti][(houji edian-1)*3+Ti2]+gd[Ti][Ti2+3];}for(Ti=0;Ti<3;Ti++){for(Ti2=0;Ti2<3;Ti2++)zhengtigangdu[(houjiedian-1)*3+Ti][(qianjiedian-1)*3+Ti2]=zhengtigangdu[(houjiedian-1)*3+Ti][(qianji edian-1)*3+Ti2]+gd[Ti+3][Ti2];}}}//**************************//输入数据的函数//**************************{char queding;//基本参数输入cout<<"请输入基本参数:单元数,节点数,支持数,有节点载荷的节点数,非节点载荷数;"<<endl<<endl;cin>>NE>>NJ>>NZ>>NPJ>>NPF;cout<<endl<<endl;//输入单元几何参数cout<<"请输入单元的几何参数:"<<endl;for(int i=0;i<NE;i++){cout<<"第"<<i+1<<"单元"<<endl;cout<<"长度:"; cin>>changdu[i];cout<<"角度:"; cin>>jiaodu[i]; jiaodu[i]=jiaodu[i]*PI/180;cout<<"弹性模量E:"; cin>>tanxingmoliang[i];cout<<"J模量:"; cin>>J_moliang[i];cout<<"面积:"; cin>>mianji[i];cout<<endl;//cout<<"第"<<i+1<<"单元长度:"<<changdu[i]<<" 角度:"<<jiaodu[i]<<" 弹性模量E:"<<tanxingmoliang[i]<<" J量:"<<J_moliang[i]<<" 面积:"<<mianji[i]<<endl<<endl;}//输入单元对应的节点总码cout<<"请输入个单元对应的节点总码"<<endl;for(i=0;i<NE;i++){cout<<"第"<<i+1<<"单元前面节点的总码:"; cin>>dym_jdm[i][0];cout<<"第"<<i+1<<"单元后面节点的总码:"; cin>>dym_jdm[i][1];}cout<<endl;//输入支承点位置cout<<"请输入支承点位置:"<<endl;for(i=0;i<NZ;i++){cout<<"第"<<i+1<<"支承位:"; cin>>zhichengweizhi[i];}cout<<endl;//输入节点载荷值cout<<"请输入节点载荷作用节点及其载荷值:"<<endl;for(i=0;i<NPJ;i++){cout<<"载荷所在节点位置:"; cin>>jdzhzuoyongdian[i];}for(i=0;i<NPJ;i++){cout<<"第"<<jdzhzuoyongdian[i]<<"节点的载荷值:"<<endl;cout<<"U:"; cin>>jiedianzaihe[i][0];cout<<"V:"; cin>>jiedianzaihe[i][1];cout<<"M:"; cin>>jiedianzaihe[i][2];}cout<<endl;cout<<"请输入非节点载荷作单元,载荷类型,载荷值,作用点(指离第一个节点的距离),作用角度:"<<endl;cout<<"其中1-均布,2-垂直集中,3-平行集中,4-力偶,5-角度集中"<<endl;for(i=0;i<NPF;i++){cout<<"第"<<i+1<<"个非节点载荷作用单元:"; cin>>fjzhzuoyongdanyuan[i];cout<<"载荷类型:"; cin>>fjzhleixing[i];cout<<"载荷值(如果是类型5的力,只输入力的绝对值) :"; cin>>fjzhzhi[i];cout<<"离第一个节点的距离:"; cin>>fjzhzuoyongdian[i];cout<<"力在该单元坐标中的作用角度:"; cin>>fjzhjiaodu[i];fjzhjiaodu[i]=fjzhjiaodu[i]*PI/180; cout<<endl;}cout<< endl<<endl;cout<<"按任意键显示输入情况:"<<endl;getch();//显示输入的情况cout<<"输入的数据如下:"<<endl;cout<<"===========================================:"<<endl;cout<<"单元数:"<<NE<<" 节点数:"<<NJ<<" 支持数:"<<NZ<<"有节点载荷的节点数:"<<NPJ<<" 非节点载荷数:"<<NPF<<endl;cout<<"============================================"<<endl;cout<<"各单元几何参数如下:"<<endl;cout<<"单元"<<"\t"<<"单元长度:"<<"\t"<<"角度: "<<"\t"<<"弹性模量E:"<<"\t"<<"J量: "<<"\t"<<"面积:"<<endl;for(i=0;i<NE;i++)cout<<i+1<<"\t"<<changdu[i]<<"\t\t"<<jiaodu[i]<<"\t\t"<<tanxingmoliang[i]<<"\t\t"<<J_moliang[i]<<"\ t\t"<<mianji[i]<<endl;cout<<"==========================================="<<endl;cout<<"各单元对应节点总码如下:"<<endl;cout<<"单元号"<<"\t"<<"前节点总码"<<"\t"<<"后节点总码"<<endl;for(i=0;i<NE;i++)cout<<i+1<<"\t\t"<<dym_jdm[i][0]<<"\t\t"<<dym_jdm[i][1]<<endl;cout<<" ==========================================="<<endl;cout<<"支承位置如下:"<<endl;cout<<"支承位置"<<"\t";for(i=0;i<NZ;i++)cout<<zhichengweizhi[i]<<" ";cout<<endl;cout<<" ==========================================="<<endl;cout<<"节点载荷输入如下:"<<endl;cout<<"节点位置"<<"\t"<<"U "<<"\t\t"<<"V "<<"\t\t"<<"M "<<endl;for(i=0;i<NPJ;i++){cout<<jdzhzuoyongdian[i]<<"\t\t"<<jiedianzaihe[i][0]<<"\t\t"<<jiedianzaihe[i][1]<<"\t\t"<<jiedianzaihe[ i][2];}cout<<" ==========================================="<<endl;cout<<"非节点载荷的作用情况如下:"<<endl;cout<<"作用单元"<<"\t"<<"载荷类型"<<"\t"<<"载荷值"<<"\t\t"<<"作用点"<<"\t\t"<<"作用角度"<<endl;for(i=0;i<NPF;i++)cout<<fjzhzuoyongdanyuan[i]<<"\t\t"<<fjzhleixing[i]<<"\t\t"<<fjzhzhi[i]<<"\t\t"<<fjzhzuoyongdian[i]<< "\t\t"<<fjzhjiaodu[i]<<endl;cout<<" ==========================================="<<endl;cout<<"按任意键继续!"<<endl;getch();cout<<endl;}///////////////////////////////////////////////////////程序结束////////////////////////////////////////////////////////////////// 用以上程序计算教材例题3-7以及3-8结果都和例题的结果一样。

相关文档
最新文档