平面四边形四节点等参单元Fortran源程序

C ************************************************

C * FINITE ELEMENT PROGRAM *

C * FOR Two DIMENSIONAL ELASticity PROBLEM *

C * WITH 4 NODE *

C ************************************************

PROGRAM ELASTICITY

character*32 dat,cch

DIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200),

& WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200)

COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN2/ N,MX,NH

COMMON /CMN3/ RF(8),SKE(8,8),NN(8)

WRITE(*,*)'PLEASE ENTER INPUT FILE NAME'

READ(*,'(A)')DAT

OPEN(4,FILE=dat,STATUS='OLD')

OPEN(7,FILE='OUT',STATUS='UNKNOWN')

READ(4,*)NP,NE,NM,NR

WRITE(7,'(A,I6)')'NUMBER OF NODE---------------------NP=',np

WRITE(7,'(A,I6)')'NUMBER OF ELEMENT------------------NE=',ne

WRITE(7,'(A,I6)')'NUMBER OF MATERIAL-----------------NM=',nm WRITE(7,'(A,I6)')'NUMBER OF surporting---------------NC=',Nr

CALL INPUT (JR,COOR,AE,MEL)

CALL CBAND (MA,JR,MEL)

DO I=1,NH

SK(I)=0.0

enddo

CALL SK0(SK,MEL,COOR,JR,MA,AE)

do I=1,N

R(I)=0.0

enddo

pause 'aaa'

stop

READ(4,*)NCP,NBE,iz

WRITE(*,'(5i8)')NCP,NBE,iz

WRITE(7,'(5i8)')NCP,NBE,iz

IF(NCP.GT.0)CALL CONCR(NCP,R,JR)

IF(NBE.GT.0) CALL BODYR(NBE,R,MEL,COOR,JR,AE) IF(iz.GT.0)then

do jj=1,iz

READ (4,*)Js,nse,(WG(I),I=1,4)

read(4,*)(iew(m),m=1,nse)

CALL FACER(iew,NSE,R,MEL,COOR,JR,WG)

enddo

endif

CALL DECOP (SK,MA)

CALL FOBA (SK,MA,R)

CALL OUTDISP(NP,R,JR)

CALL STRESS (COOR,MEL,JR,AE,R,STRE)

WRITE(7,'(A)')' PROGRAM SAFF HAS BEEN ENDED'

WRITE(*,'(A)')' PROGRAM SAFF HAS BEEN ENDED'

STOP

c RETURN

END

C *********************************************

SUBROUTINE INPUT (JR,COOR,AE,MEL)

DIMENSION JR(2,*),COOR(2,*),AE(4,*),MEL(5,*)

COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN2/ N,MX,NH

DO 70 I=1,NP

READ(4,*) IP,X,Y

COOR(1,IP)=X

COOR(2,IP)=Y

70 CONTINUE

DO 11 J=1,NE

READ(4,*)NEE,NME,(MEL(I,NEE),I=1,4)

MEL(5,NEE)=NME

11 CONTINUE

DO 10 I=1,NP

DO 10 J=1,2

10 JR(J,I)=1

DO 20 I=1,NR

READ(4,*) IP,IX,IY

JR(1,IP)=IX

JR(2,IP)=IY

20 CONTINUE

N=0

DO 30 I=1,NP

DO 30 J=1,2

IF (JR(J,I)) 30,30,25

25 N=N+1

JR(J,I)=N

30 CONTINUE

DO 55 J=1,NM

READ (4,*)JJ,(AE(I,JJ),I=1,4)

WRITE(*,910) JJ,(AE(I,JJ),I=1,4)

55 CONTINUE

910 FORMAT (/20X,'MATERIAL PROPERTIES'/(3X,I5,4(1x,E8.3))) RETURN

END

C **********************************************

SUBROUTINE CBAND (MA,JR,MEL)

DIMENSION MA(*),JR(2,*),MEL(5,*),NN(8)

COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN2/ N,MX,NH

DO 65 I=1,N

65 MA(I)=0

DO 90 IE=1,NE

DO 75 K=1,4

IEK=MEL(K,IE)

DO 95 M=1,2

JJ=2*(K-1)+M

NN(JJ)=JR(M,IEK)

95 CONTINUE

75 CONTINUE

L=N

DO 80 I=1,2*4

NNI=NN(I)

IF(NNI.EQ.0) GO TO 80

IF(NNI.LT.L) L=NNI

80 CONTINUE

DO 85 M=1,2*4

JP=NN(M)

IF(JP.EQ.0) GO TO 85

JPL=JP-L+1

IF(JPL.GT.MA(JP)) MA(JP)=JPL

85 CONTINUE

90 CONTINUE

MX=0

MA(1)=1

DO 10 I=2,N

IF(MA(I).GT.MX) MX=MA(I)

MA(I)=MA(I)+MA(I-1)

10 CONTINUE

NH=MA(N)

WRITE(7,'(A,I8)')'TOTAL DEGREES OF FREEDOM-----------N= ',N

WRITE(7,'(A,I8)')'MAX-SEMI-BANDWIDTH-----------------MX=',MX WRITE(7,'(A,I8)')'TOTAL-STORAGE----------------------NH=',NH 500 FORMAT (/5X,'FREEDOM N='

*,I5,3X,'SEMI-BANDWI. MX=',I5,3X,

* 'STORAGE NH=',I7)

RETURN

END

C **********************************************

SUBROUTINE SK0(SK,MEL,COOR,JR,MA,AE)

DIMENSION SK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*), * AE(4,*),XYZ(2,4),iven(4)

COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN2/ N,MX,NH

COMMON /CMN3/ RF(8),SKE(8,8),NN(8)

COMMON /CMN4/ NEE,NME

COMMON /GAUSS/ RSTG(3),H(3)

H(1)=0.5555555555555560

H(2)=0.8888888888888890

H(3)=H(1)

RSTG(1)=-0.7745966692414830

RSTG(2)=0.00

RSTG(3)=-RSTG(1)

DO 10 IE=1,NE

NEE=IE

NME=MEL(5,IE)

DO 75 K=1,4

IEK=MEL(K,IE)

iven(k)=IEK

DO 95 M=1,2

JJ=2*(K-1)+M

NN(JJ)=JR(M,IEK)

95 XYZ(M,K)=COOR(M,IEK)

75 CONTINUE

CALL STIF(XYZ,AE,iven)

DO 60 I=1,8

DO 60 J=1,8

II=NN(I)

JJ=NN(J)

IF ((JJ.EQ.0).OR.(II.LT.JJ)) GO TO 60

JN=MA(II)-(II-JJ)

SK(JN)=SK(JN)+SKE(I,J)

60 CONTINUE

70 CONTINUE

write(7,1111) ((ske(i,j),j=1,8),i=1,8)

1111 format(2x,8f12.2)

10 CONTINUE

RETURN

END

C *********************************************

SUBROUTINE STIF(XYZ,AE,iven)

DIMENSION AE(4,*),DNX(2,4),XYZ(2,*),iven(*),

* RJAC(2,2)

COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN2/ N,MX,NH

COMMON /CMN3/ RF(8),SKE(8,8),NN(8)

COMMON /CMN4/ NEE,NME

COMMON /GAUSS/ RSTG(3),H(3)

DO 40 I=1,8

RF(I)=0.00

DO 30 J=1,8

SKE(I,J)=0.00

30 CONTINUE

40 CONTINUE

E=AE(1,NME)

U=AE(2,NME)

GAMA=AE(3,NME)

D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U))

D2=E*U/((1.00+U)*(1.00-2.00*U))

D3=E*0.50/(1.00+U)

DO 120 I=1,4

II=2*(I-1)

I1=II+1

I2=II+2

DO 115 J=1,4

JJ=2*(J-1)

J1=JJ+1

J2=JJ+2

DXX=0

DXY=0

DYX=0

DYY=0

DO 99 IS=1,3

S=RSTG(IS)

SH=H(IS)

DO 98 IR=1,3

R=RSTG(IR)

RH=H(IR)

CALL FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DNIX=DNX(1,I)

DNIY=DNX(2,I)

DNJX=DNX(1,J)

DNJY=DNX(2,J)

DXX=DXX+DNIX*DNJX*DET*RH*SH

DXY=DXY+DNIX*DNJY*DET*RH*SH

DYX=DYX+DNIY*DNJX*DET*RH*SH

DYY=DYY+DNIY*DNJY*DET*RH*SH

98 CONTINUE

99 CONTINUE

SKE(I1,J1)=DXX*D1+DYY*D3

SKE(I2,J2)=DYY*D1+DXX*D3

SKE(I1,J2)=DXY*D2+DYX*D3

SKE(I2,J1)=DYX*D2+DXY*D3

115 CONTINUE

120 CONTINUE

RETURN

END

C *********************************************

SUBROUTINE CONCR(NCP,R,JR)

DIMENSION R(*),JR(2,*),XYZ(2)

DO 100 I=1,NCP

READ (4,*) IP,PX,PY

XYZ(1)=PX

XYZ(2)=PY

DO 95 J=1,2

L=JR(J,IP)

IF(L.EQ.0) GO TO 95

R(L)=R(L)+XYZ(J)

95 CONTINUE

100 CONTINUE

RETURN

END

C **********************************************

SUBROUTINE BODYR(NBE,R,MEL,COOR,JR,AE)

DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),

& AE(4,*),XYZ(2,4),iven(4)

COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN2/ N,MX,NH

COMMON /CMN3/ RF(8),SKE(8,8),NN(8)

COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)

COMMON /GAUSS/ RSTG(3),H(3)

H(1)=1.0

H(2)=1.0

RSTG(1)=-0.5773502691896260

RSTG(2)=-RSTG(1)

DO 10 IE=1,NBE

DO I=1,8

RF(I)=0.00

ENDDO

c READ(4,*)NEE

NEE=ie

NME=MEL(5,NEE)

GAMA=AE(3,NME)

DO 75 K=1,4

IEK=MEL(K,NEE)

iven(k)=iek

DO 95 M=1,2

JJ=2*(K-1)+M

NN(JJ)=JR(M,IEK)

95 XYZ(M,K)=COOR(M,IEK)

75 CONTINUE

DO 99 IS=1,2

S=RSTG(IS)

SH=H(IS)

DO 98 IR=1,2

RR=RSTG(IR)

RH=H(IR)

CALL FUN8 (XYZ,RR,S,DET)

DO 30 I=1,4

J=2*I

RF(J)=RF(J)-FUN(I)*RH*SH*DET*GAMA 30 CONTINUE

98 CONTINUE

99 CONTINUE

CALL ASLOAD (R)

10 CONTINUE

RETURN

END

C *********************************************

SUBROUTINE FACER(iew,NSE,R,MEL,COOR,JR,WG)

DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),wg(*)

* ,XYZ(2,4),iew(*),PR(2)

COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN2/ N,MX,NH

COMMON /CMN3/ RF(8),SKE(8,8),NN(8)

COMMON /CMN4/ NEE,NME

COMMON /GAUSS/ RSTG(3),H(3)

H(1)=1.0

H(2)=1.0

RSTG(1)=-0.5773502691896260

RSTG(2)=-RSTG(1)

nwf=0

nnf=0

ir=wg(1)+0.1

if(ir.eq.1)nwf=1

if(ir.eq.2)nnf=1

DO 510 IE=1,NSE

DO I=1,8

RF(I)=0.00

ENDDO

nee=iew(ie)

DO 575 K=1,4

IEK=MEL(K,NEE)

DO 595 M=1,2

JJ=2*(K-1)+M

NN(JJ)=JR(M,IEK)

595 XYZ(M,K)=COOR(M,IEK)

575 CONTINUE

IF(NWF.EQ.1) then

GAMA=WG(2)

Z0=WG(3)

NSU=WG(4)+0.1

CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,1)

endif

IF(NNF.EQ.1) then

q=WG(2)

NSU=WG(4)+0.1

do j=1,2

PR(J)=q

enddo

CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,2)

endif

CALL ASLOAD (R)

510 CONTINUE

RETURN

END

C *********************************************

SUBROUTINE SURLOD (NSU,XYZ,PR,Z0,GAMA,NSI)

DIMENSION XYZ(2,*),RST(3),PR(2),KCRD(4),KFACE(2,4), & FVAL(4),NODES(2),FACT(4)

COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN2/ N,MX,NH

COMMON /CMN3/ RF(8),SKE(8,8),NN(8)

COMMON /CMN4/ NEE,NME

COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)

COMMON /GAUSS/ RSTG(3),H(3)

DATA KCRD/1,1,2,2/

DATA KFACE/1, 4,

* 2, 3,

* 1, 2,

* 4, 3/

DATA FVAL/-1.00,1.00,-1.00,1.00/ FACT(1)=1.0

FACT(2)=-1.0

FACT(3)=-1.0

FACT(4)=1.0

FACTNUS=FACT(NSU)

DO I=1,2

J=KFACE(I,NSU)

NODES(I)=J

ENDDO

IF (NSI.EQ.1) THEN

DO I=1,2

J=NODES(I)

Z=Z0-XYZ(2,J)

PR(I)=0.00

IF (Z.GT.0.00) PR(I)=Z*GAMA ENDDO

ENDIF

ML=KCRD(NSU)

IF(ML.EQ.1)MM=2

IF(ML.EQ.2)MM=1

RST(ML)=FVAL(NSU)

DO 70 LX=1,2

RST(MM)=RSTG(LX)

CALL FUN8 (XYZ,RST(1),RST(2),DET)

PXYZ=0.00

DO 25 I=1,2

J=NODES(I)

PXYZ=PXYZ+FUN(J)*PR(I)

25 CONTINUE

A1=XJAC(MM,2)

A2=-XJAC(MM,1)

30 DO 60 I=1,2

J=NODES(I)

K2=2*J

K1=K2-1

Q=PXYZ*FUN(J)*H(LX)*FACTNUS

RF(K1)=RF(K1)+Q*A1

RF(K2)=RF(K2)+Q*A2

60 CONTINUE

70 CONTINUE

RETURN

END

C *********************************************

SUBROUTINE ASLOAD (R)

DIMENSION R(*)

COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN3/ RF(8),SKE(8,8),NN(8)

DO 20 I=1,8

L=NN(I)

IF (L.EQ.0) GO TO 20

R(L)=R(L)+RF(I)

20 CONTINUE

RETURN

END

C ***********************************************

SUBROUTINE DECOP (SK,MA)

DIMENSION SK(*),MA(*)

COMMON /CMN2/ N,MX,NH

DO 50 I=2,N

L=I-MA(I)+MA(I-1)+1

K=I-1

L1=L+1

IF (L1.GT.K) GO TO 30

DO 20 J=L1,K

IJ=MA(I)-I+J

M=J-MA(J)+MA(J-1)+1

IF (L.GT.M) M=L

MP=J-1

IF (M.GT.MP) GO TO 20

DO 10 LP=M,MP

IP=MA(I)-I+LP

JP=MA(J)-J+LP

SK(IJ)=SK(IJ)-SK(IP)*SK(JP)

10 CONTINUE

20 CONTINUE

30 IF (L.GT.K) GO TO 50

DO 40 LP=L,K

IP=MA(I)-I+LP

LPP=MA(LP)

SK(IP)=SK(IP)/SK(LPP)

II=MA(I)

SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP) 40 CONTINUE

50 CONTINUE

等参单元

5．等参单元本章包括以下内容： 5.1等参单元的基本概念 5.2四边形八节点等参单元 5.3等参单元的单元分析 5.4六面体等参单元 5.1等参单元的基本概念在进行有限元分析时，单元离散化会带来计算误差，主要采用两种方法来降低单元离散化产生的误差：1）提高单元划分的密度，被称为h 方法（h-method ）；2）提高单元位移函数多项式的阶次，被称为p 方法（p-method ）。在平面问题的有限单元中，我们可以选择四结点的矩形单元，如图5-1所示，该矩形单元在x 及y 方向的边长分别为2a 和2b 。图5-1 四结点矩形单元同第三章的方法类似，将单元的位移模式选为， xy a y a x a a u 4321+++= xy a y a x a a v 8765+++= （5-1）可得到， p p m m j j i i u N u N u N u N u +++= p p m m j j i i v N v N v N v N v +++= （5-2）形态函数为， )1)(1(41b y a x N i --= )1)(1(41b y a x N j - + =

)1)(1(41b y a x N m ++= )1)(1(4 1b y a x N p + - = （5-3）上述单元位移模式满足位移模式选择的基本要求： 1）反映了单元的刚体位移和常应变， 2）单元在公共边界上位移连续。在矩形单元的边界上，坐标x 和y 的其中一个取常量，因此在边界上位移是线性分布的，由两个结点上的位移确定。与三结点三角形单元相比，四结点矩形单元的位移模式是坐标的二次函数，能够提高计算精度，但也有显著的缺点，两种单元的比较如下。表5-1 三结点三角形单元与四结点矩形单元比较如果任意形状的四边形四结点单元采用矩形单元的位移模式，则在公共边界上不满足位移连续性条件。为了既能得到较高的计算精度，又能适应复杂的边界形状，可以采用坐标变换。图5-2任意四结点四边形单元图5-3四结点正方形单元在图5-2所示的任意四边形单元上，用等分四条边的两族直线分割四边形，以两族直线的中心为原点，建立局部坐标系),(ηξ，沿ξ及η增大的方向作为ξ轴和η轴，并令四条边上的ξ及η值分别为1±。为了求出位移模式，以及局部坐标与整体坐标之间的变换式，在局部坐标系中定义一个四结点正方形单元，如图5-3所示。参照矩形单元，四结点正方形单元的位移模式为， 44332211u N u N u N u N u +++= 44332211v N v N v N v N v +++= （5-4）

平面四节点等参单元matlab实现

计算力学报告平面四节点等参单元学生姓名：朱彬学号：20100311

一、问题描述及分析在无限大平面内有一个小圆孔。孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=0.5mm，在远离大圆边界的地方模型是比较精确的。由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。二、有限元划分描述在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。具体操作如下：打开ansys，在单元类型中选择solid->Quad 4 node 182单元；建立类半径为0.5外半径为10的圆环；使用mashtool中的智能划分和将单元退化成三角形单元；使用工具栏中List中的Nodes和Elements 选项将节点和单元数据导出并导入Excle中，总共得到了207个单元和229个节点。如下图：

图1 三、有限元程序及求解程序求解使用了matlab语言。具体如下：程序： clc clear E=2e11; %弹性模量 NU=0.3; %泊松比 t=0.1; %厚度 X=xlsread('D:\data','nodes'); %读取节点坐标 elem=xlsread('D:\data','elements'); %读取单元编号 w=[1,2,3,4,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28]; %有位移

约束的节点 n=size(X,1); %节点数 m=size(elem,1); %单元数 K=zeros(2*n); %初始总体刚度矩阵 for i=1:m syms Ks Et x y I1 I2 a b B; %定义可能存在的变量 e=[1,1;-1,1;-1,-1;1,-1]; for j=1:4 %形函数 N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et); end x=0;y=0; for j=1:4 %标准母单元映射到真实单元x=x+N(j)*X(elem(i,j),1); y=y+N(j)*X(elem(i,j),2); end J1=jacobian([x;y],[Ks;Et]); %雅克比矩阵及其转置 J=J1'; for j=1:4 I1=diff(N(j),Ks); %形函数分别对Ks和Et的偏导数 I2=diff(N(j),Et); C=(J^-1)*[I1;I2]; a=C(1); %形函数对x，y的偏导数 b=C(2); B(1,2*j-1)=a; %组成B阵 B(1,2*j)=0; B(2,2*j-1)=0; B(2,2*j)=b; B(3,2*j-1)=b; B(3,2*j)=a; end D=(E/(1-NU*NU))*[1,NU,0;NU,1,0;0,0,(1-NU)/2]; %D阵 k=zeros(8,8); Kss=[-0.906179,-0.538469,0,0.538469,0.906179]; %5*5高斯积分点 Ett=[-0.906179,-0.538469,0,0.538469,0.906179]; H=[0.236926,0.478628,0.568888,0.478628,0.236926];%高斯积分权系数 for j=1:5 %高斯积分求单元刚度阵 for l=1:5 Ks=Kss(j);Et=Ett(l); B=subs(B); J=subs(J);

平面四节点等参单元分析程序

变分原理与有限元大作业平面四节点等参单元分析程序；姓名：潘清【学号：SQ 完成时间：2011-4-26 | 一、概述

! 通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。 C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。二、编程思想 ; 本程序采用C语言编程，编制平面四边形四节点等参元程序，用以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移，然后求中间高斯点的应力，最后用绕节点平均法讲单元应力等效到节点上，再将结果写到tecplot文件中。在有限元程序中，变量数据需赋值的可分为节点信息，单元信息，载荷信息等。对于一个节点来说，需以下信息：节点编号（整型），节点坐标（实型），节点已知位移（实型），节点载荷（实型），边界条件（实型）等。同样，对于一个单元来说，需以下信息：单元的

平面四边形四节点等参单元Fortran源程序

C ************************************************ C * FINITE ELEMENT PROGRAM * C * FOR Two DIMENSIONAL ELASticity PROBLEM * C * WITH 4 NODE * C ************************************************ PROGRAM ELASTICITY character*32 dat,cch DIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200), & WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) WRITE(*,*)'PLEASE ENTER INPUT FILE NAME' READ(*,'(A)')DAT OPEN(4,FILE=dat,STATUS='OLD') OPEN(7,FILE='OUT',STATUS='UNKNOWN') READ(4,*)NP,NE,NM,NR WRITE(7,'(A,I6)')'NUMBER OF NODE---------------------NP=',np WRITE(7,'(A,I6)')'NUMBER OF ELEMENT------------------NE=',ne WRITE(7,'(A,I6)')'NUMBER OF MATERIAL-----------------NM=',nm WRITE(7,'(A,I6)')'NUMBER OF surporting---------------NC=',Nr CALL INPUT (JR,COOR,AE,MEL) CALL CBAND (MA,JR,MEL) DO I=1,NH SK(I)= enddo CALL SK0(SK,MEL,COOR,JR,MA,AE) do I=1,N R(I)= enddo pause 'aaa' stop READ(4,*)NCP,NBE,iz WRITE(*,'(5i8)')NCP,NBE,iz WRITE(7,'(5i8)')NCP,NBE,iz IF CONCR(NCP,R,JR) IF CALL BODYR(NBE,R,MEL,COOR,JR,AE) IF do jj=1,iz READ (4,*)Js,nse,(WG(I),I=1,4) read(4,*)(iew(m),m=1,nse) CALL FACER(iew,NSE,R,MEL,COOR,JR,WG) enddo endif CALL DECOP (SK,MA)

平面四边形4结点等参有限单元法

有限元程序设计平面四边形4结点等参有限单元法程序设计 1、程序功能及特点

b.前处理采用网格自动划分技术，自动生成单元及结点信息。 b.能计算受集中力、自重体力、分布面力和静水压力的作用。 c.计算结点的位移和单元中心点的应力分量及其主应力。 d.后处理采取整体应力磨平求得各个结点的应力分量。 e.算例计算结果与ANSYS计算结果比较，并给出误差分析。 f.程序采用Visual Fortran 5.0编制而成。 2、程序流程及图框图2-１程序流程图

图2-２子程序框图其中，各子程序的主要功能为： INPUT――输入原始数据 HUAFEN――自动网格划分，形成COOR(2,NP),X,Y的坐标值与单元信息CBAND――形成主元素序号指示矩阵MA（*） SKO――形成整体刚度矩阵[K] CONCR――计算集中力引起的等效结点荷载{R}e BODYR――计算自重体力引起的等效结点荷载{R}e FACER――计算分布面力引起的等效结点荷载{R}e DECOP――支配方程LU三角分解 FOBA――LU分解直接解法中的回代过程 OUTDISP――输出结点位移分量 STRESS――计算单元应力分量 OUTSTRE――输出单元应力分量 STIF――计算单元刚度矩阵 FDNX――计算形函数对整体坐标的导数 T i i y N x N ? ? ? ? ? ? ? ? ? ? ，= i1,2,3,4。 FUN8――计算形函数及雅可比矩阵[J] SFUN ――应力磨平-单元下的‘K’=NCN‘ SCN――应力磨平-单元下的右端项系数‘CN‘ SUMSKN――应力磨平-单元下的右端项集成到总体的‘P‘ SUMSTRS――应力磨平-单元下的集成到总体的‘K‘ GAUSTRSS――高斯消元求磨平后的应力 3、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。在运行程序之前，根据程序中INPUT需要的数据输入建立一个存放原始数据的文件，这个文件的名

平面四节点等参单元分析程序

变分原理与有限元大作业平面四节点等参单元分析程序姓名：潘清学号：SQ10018014033 完成时间：2011-4-26

一、概述通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。 C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。二、编程思想本程序采用C语言编程，编制平面四边形四节点等参元程序，用以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移，然后求中间高斯点的应力，最后用绕节点平均法讲单元应力等效到节点上，再将结果写到tecplot文件中。在有限元程序中，变量数据需赋值的可分为节点信息，单元信息，载荷信息等。对于一个节点来说，需以下信息：节点编号（整型），节点坐标（实型），节点已知位移（实型），节点载荷（实型），边界条件（实型）等。同样，对于一个单元来说，需以下信息：单元的节点联接信息（整型），材料信息（弹性模量，泊松比等）（实型）等。

八节点等参单元平面有限元

八节点等参单元平面有限元分析程序土木工程学院

目录 1.概述 (1) 2.编程思想 (2) 1 (2) 2 (2) 2.1.八节点矩形单元介绍 (2) 1 (5) 2 (5) 2.1 (5) 2.2.有限元分析的模块组织 (5) 3.程序变量及函数说明 (6) 3. (6) 3 (6) 3.1.主要变量说明： (6) 3.1. (7) 3.2.主要函数说明 (7) 4.程序流程图 (8) 5.程序应用与ANSYS分析的比较 (9) 4 (9) 5 (9) 5.1.问题说明 (9) 5.2.ANSYS分析结果 (10) 5.3.自编程序分析结果 (12) 5.4.结果分析比较 (12) 参考文献 (15) 附录程序源代码 (16)

《计算力学》课程大作业 1.概述通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。 C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。

平面四边形四节点等参单元Fortran源程序复习过程

平面四边形四节点等参单元F o r t r a n源程序

C ************************************************ C * FINITE ELEMENT PROGRAM * C * FOR Two DIMENSIONAL ELASticity PROBLEM * C * WITH 4 NODE * C ************************************************ PROGRAM ELASTICITY character*32 dat,cch DIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200), & WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) WRITE(*,*)'PLEASE ENTER INPUT FILE NAME' READ(*,'(A)')DAT OPEN(4,FILE=dat,STATUS='OLD') OPEN(7,FILE='OUT',STATUS='UNKNOWN') READ(4,*)NP,NE,NM,NR WRITE(7,'(A,I6)')'NUMBER OF NODE---------------------NP=',np WRITE(7,'(A,I6)')'NUMBER OF ELEMENT------------------NE=',ne WRITE(7,'(A,I6)')'NUMBER OF MATERIAL-----------------NM=',nm WRITE(7,'(A,I6)')'NUMBER OF surporting---------------NC=',Nr CALL INPUT (JR,COOR,AE,MEL) CALL CBAND (MA,JR,MEL) DO I=1,NH SK(I)=0.0 enddo CALL SK0(SK,MEL,COOR,JR,MA,AE) do I=1,N R(I)=0.0 enddo pause 'aaa' stop READ(4,*)NCP,NBE,iz WRITE(*,'(5i8)')NCP,NBE,iz WRITE(7,'(5i8)')NCP,NBE,iz IF(NCP.GT.0)CALL CONCR(NCP,R,JR) IF(NBE.GT.0) CALL BODYR(NBE,R,MEL,COOR,JR,AE) IF(iz.GT.0)then do jj=1,iz READ (4,*)Js,nse,(WG(I),I=1,4) read(4,*)(iew(m),m=1,nse) CALL FACER(iew,NSE,R,MEL,COOR,JR,WG) enddo endif

平面四节点等参单元分析程序

\ 变分原理与有限元大作业平面四节点等参单元分析程序 . 姓名：潘清学号：SQ 完成时间：2011-4-26 ：一、概述

通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。 C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。二、，三、编程思想本程序采用C语言编程，编制平面四边形四节点等参元程序，用以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移，然后求中间高斯点的应力，最后用绕节点平均法讲单元应力等效到节点上，再将结果写到tecplot文件中。在有限元程序中，变量数据需赋值的可分为节点信息，单元信息，载荷信息等。对于一个节点来说，需以下信息：节点编号（整型），节点坐标（实型），节点已知位移（实型），节点载荷（实型），边界条件（实型）等。同样，对于一个单元来说，需以下信息：单元的

四边形四节点等参元matlab程序

悬臂钢梁，尺寸如图一所示；v=0.3。h=1，E=2.1e11. 图一悬臂钢梁图二单元划分与结点编号

附录： %---------------四边形四节点等参元matlab计算程序---------------------------- % 2013年 % 13级建筑与土木工程Bruce % B15-405 % 本程序只输出了节点位移、单元中心点的应力 %******************************************************************* % 变量说明 % E v h % 弹性模量泊松比厚度 % NPOIN NELEM NVFIX NNODE NFPOIN % 总结点数, 单元数, 约束结点个数, 单元节点数,受力结点数 % COORD LNODS % 结构节点整体坐标数组, 单元定义数组, % FPOIN FORCE FIXED % 结点力数组，总体荷载向量, 约束信息数组 % HK DISP % 总体刚度矩阵,结点位移向量 %****************************** clear all format short e FP1=fopen(' sjd.txt','rt'); %打开数据文件 %%读入控制数据 E=fscanf(FP1,'%f',1); %弹性模量 v=fscanf(FP1,'%f',1); % 泊松比 h=fscanf(FP1,'%f',1); %厚度 NELEM=fscanf(FP1,'%d',1); %单元数 NPOIN=fscanf(FP1,'%d',1); % 总结点数 NNODE=fscanf(FP1,'%d',1); %单元节点数 NFPOIN=fscanf(FP1,'%d',1); %受力结点数 NVFIX=fscanf(FP1,'%d',1); %约束结点个数 LNODS=fscanf(FP1,'%f',[NNODE,NELEM])'; % 单元定义：单元结点号（逆时针）COORD=fscanf(FP1,'%f',[2,NPOIN])'; % 结点号x,y坐标(整体坐标下) FPOIN=fscanf(FP1,'%f',[3,NFPOIN])'; % 节点力：结点号、X方向力(向右正)，Y方向力(向上正) FIXED=fscanf(FP1,'%d',[3,NVFIX])'; %约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号

平面四边形四节点等参单元Fortran源程序

等参单元

平面四节点等参单元matlab实现

平面四节点等参单元分析程序

最新平面四边形4结点等参有限单元法

平面四边形四节点等参单元Fortran源程序

平面四边形4结点等参有限单元法

平面四节点等参单元分析程序

八节点等参单元平面有限元

平面四边形四节点等参单元Fortran源程序复习过程

平面四节点等参单元分析程序

最新第五章-等参单元1

四边形四节点等参元matlab程序