Fortran Library for FEM Programs Development

Fortran library for FEM programs development

Zitny R., CTU FME Prague, March 2001

Introduction

The following examples make use the FRONTAL METHOD for assembly and simultaneous solution of equations derived by finite element methods. The frontal method is implemented in subroutines of more than 15 years old fortran library INTERFOR[1], designed only for the FEM principles teaching, and not for development of efficient and large scale programs. Primary aim is simplicity and transparency of programs.

XMKPE1Plane stress analysis

XMKPE3Flow in pipeline networks (laminar/turbulent)

XMKPE5Temperature field in a liquid (cylindrical system, prescribed velocity field)

XMKPE7Rotationally symmetric shell

XMKPE9Beams and trusses

XMKPE10Pressure thermal forming (waffle baking)

Specific features of programs

All data - matrices, (coordinates of nodes, connectivity matrix, nodal parameters, properties of elements, stiffness matrices and so on) are dynamically allocated in the COMMON // area. In most of the examples the following lay-out of matrices in the COMMON area will be used:

1.matrixdirectory of all matrices in the COMMON area (each column of this directory describes dimensions, mapping function, etc by 9 integer numbers)

2.matrixNOE connectivity matrix (nodes of elements)

3.matrixXYZ coordinates of nodes

4.matrixIU indices of DOF corresponding to nodes

5.matrixU nodal parameters (vector of unknowns)

It is assumed that the vector U contains not only the unknown nodal parameters, but also the parameters prescribed as a boundary condition. To distinguish between them the following convention is adopted: if IUn,i> 0 the i-th parameter of n-th node is the calculated unknown, while if IUn,i <0 the parameter is prescribed as a boundary condition. In the both cases the absolute value |IUn,i| is an index of parameter in the vector U.

The most important subroutines are

MEFRIN(noe,iu)

prefrontal algorithm, which identifies the last occurrence of a node in the connectivity matrix. NOE – index of connectivity matrix (usually the second matrix in COMMON area, therefore NOE=2, see above), and IU is the index of matrix containing indices of DOF (usually IU=4, see above). The purpose of prefront is to identify the moment when all equations concerning a certain node have been assembled and therefore the moment when the corresponding nodal parameters can be eliminated even if the assembly of global system of equations has not been completed. The last occurrence of a node in the list of elements is indicated by a negative sign of its index in the connectivity matrix.

MEFRON(noe,iu,u,nfront,nloc,subloc,iwrite)

frontal solver, assembling local matrices and right hand side vectors defined in the subroutine SUBLOC, and performs the elimination of unknowns as soon as the appropriate equation has been completed[2]. In simple words, the MEFRON routine performs assembly of element matrices and the solution of system simultaneously. The first three parameters of subroutine MEFRON are indices of matrices (usually 2,4,5), NFRONT is maximum width of front, i.e. maximum dimension of submatrix containing equations which can not be eliminated in advance, NLOC is the maximum dimension of local element matrix and IWRITE is the number of logical device used for storing temporary global matrix. If IWRITE=0 the matrix will be allocated in the COMMON area.

Subroutine SUBLOC defines a specific problem and must have the following parameters

SUBROUTINE subloc(ie,nloc,nu,a,r)

where input parameters are IE (index of element), NLOC (max.dimension of element matrix), and NU (number of nodes of element). Output parameters are A(NLOC x NLOC) element matrix and R(NLOC) vector of loads, i.e. the right hand side vector.

Auxiliary subroutines

MTOOLS (maxdir,subkey)

management of data in COMMON area – interactive definition and modification of matrices using MATRIX EDITOR. MAXDIR-size of directory, maximum number of matrices in the COMMON area. Specifying MAXDIR=0 defines in fact directory for 32 matrices (this is the smallest size of directory). SUBKEY is the name of user defined subroutine activated after pressing Enter, function keys, arrows etc:

SUBROUTINE SUBKEY(imat,none,none,none,key)

Only the first and the last parameters are important inputs: IMAT- index of active matrix, KEY identify the entry point (KEY=51, SUBKEY has been called after pressing arrow keys, Enter, F1,F2,..., =52 before return, =53 after entry into MTOOLS, =54 before any keyboard scanning, =55 after pressing F10).

Another auxiliary routines, see also examples

MTOBIR(file,ierr)- reading common// from FILE in a binary format

MINPR(ienv,label,nchar,irow,icol,value)- input real VALUE

MLOCR(imat,irow,icol, loc,ityp,length,istep)position LOC of a matrix entry in COMMON

SKBDR(inkey)read character from keyboard

Examples

All programs describing selected examples have the same basic structure

c Main program – reads COMMON, adjust input data, and calls frontal method routines

call mtobir('xmkpe.mat',iend)

call mtools(0,subkey)

call mefrin(2,4)

call mefron(2,4,5,50,4,subloc,0)

call mtools(0,subkey)

end

subroutine subloc(ie,nl,nu,a,r)

dimension a(nl,nl),r(nl)

c calculation of local matrix a and right hand side vector for the element IE

c data must be transferred from COMMON//

...

end

3.1.XMKPE1 – plane stress

Displacement ux, uy in the x ,y direction are primary unknown variables. In terms of ux, uy the internal energy of deformation, as well as potential energy of external forces have to be expressed. To do this, three component of deformation tensor and three component of stress tensor have to calculated as a function of nodal displacements for a general element (triangle, qudrilateral,... having arbitrary number of nodes and shape functions Ni(x,y))

Deformation tensor has in fact 4 components, however ezz has no contribution to energy

(3.1.1)

Stresses are related to deformation according to Hook's law for the case of plane stress (szz=0)

(3.1.2)

Energy of deformation

(3.1.3)

Stiffness matrix K is therefore

(3.1.4)

Kmatrix is composed of submatrices Kij of dimension 2 x 2 corresponding to the combination of nodes i and j

(3.1.5)

Potential energy of external singular forces acting in nodes is simply

(3.1.6)

Overall potential energy is therefore

(3.1.7)

and minimum of this potential energy (Lagrangian principle) is achieved at

(3.1.8)

From now on we have to take into account specific form of element. The XMKPE1 program can be designed for different form of elements, for example 3 and 6 nodes triangles, and 4 and 8 nodes quadrilaterals.

The simplest case is the 3-node triangle with linear shape functions; their derivatives are constant, therefore no numerical integration is necessary.

(3.1.9)

where Sis the area of triangle and

(3.1.10)

More complicated are isoparametric elements.

(3.1.11)

In terms of these function the coordinate transformation is defined as

(3.1.12)

Derivatives are transformed

(3.1.13)

and partial derivatives with respect x,y can be obtained by inversion of Jacobi transformation matrix

(3.1.14)

where det J , called Jacobian, is determinant of Jacobi matrix. The Jacobian transforms area of element in the (,) plane to the area in the (x,y) plane, which makes integration easier

(3.1.15)

These equations are used in the following program XMKPE1.for

c MKPE1 @Zitny 4/92

c Plane stress analysis using

c a) 3-nodal triangles,

c b) 4/8-nodal isoparametric quadrilateral elements

c Matrices:

c 2 - noe [connectivity matrix]

c 3 - xy [x,y]

c 4 - ipu [indices of nodal parameters]

c 5 - rpu [nodal parameters u,v]

c 6 -

c 7 - help [text]

common /prop/E,rmi /scom/iss(50)

external lokal,sext

character *4 inkey

equivalence (kbdr,iss(13))

E=2.1e11

rmi=0.3

call scomin(-1)

call sintro('xmkpe1.for')

call mtobir('xmkpe1.mat',iend)

10 call sclear

call minpr(0,' Modulus of elasticity E=',25,14,2,E)

call minpr(0,' Poisson constant mi=',25,15,2,rmi)

1call mtools(0,sext)

c IMI 2821 IMA 20000 MTOOLS

c | Directory of matrices |

c | 1 MAP=2 2 MAP=1 3 MAP=1 4 MAP=1 5 MAP=1 |

c | 1> 288 289> 1088 1089> 1888 1889> 2288 2289> 2688 |

c |-IIIIIII------IIII------RR------II------RR------|

c | IIIIIII IIII RR II RR |

c | IIIIIII IIII_ 4 RR II RR |

c | IIIIIII_ 9 IIII RR_ 9 II_ 9 RR_ 9 |

c | IIIIIII _ 4 RR II RR |

c | IIIIIII RR II RR |

c | _ 7 _ 2 _ 2 _ 2 |

c | |

c | 7 MAP=1 6 MAP=0 |

c | 2689> 2821 4 5 3 |

c |------TTTTTTTTTTT----IIRRRR------|

c | TTTTTTTTTTT IIRRRR |

c | Please prepare new data: TTTTTTTTTTT IIRRRR |

c | Created by COSMOS [ F10 ] TTTTTTTTTTT_ 7IIRRRR_ 9 |

c | or manually: TTTTTTTTTTT IIRRRR |

c | Select matrix and [ Enter ] TTTTTTTTTTT IIRRRR |

c | Start: TTTTTTTTTTT IIRRRR |

c | Computation - by [ F9 ] _ 19 _ 6 |

c F1-SED F2-Save F3-Load F4-Dim F5-Copy F6-Del F7-New F8-Zero F9-Exit F0-Ext

if(kbdr.eq.67)then

call mefrin(2,4)

c maximum dimension of global matrix 60x60, and local matrix 16x16 (for 8-node quadrilaterals)

call mefron(2,4,5,60,16,lokal,0)

endif

call grid(kbd)

if(kbd.eq.62) goto 10

if(kbd.ne.27) goto 1

call scwrit(16,2,28,23,' Do You want to Quit ? [Y/N]')

call skbdr(inkey)

if(inkey.ne.'Y '.and.inkey.ne.'y ')goto 1

999 call sclear

end

subroutine lokal(ie,nl,nod,a,b)

c Input: IE-index of element, NL-dimension, NOD-number of nodes

dimension a(nl,nl),b(nl),

/ x(8),y(8),g(2),w(2),f(8),fx(8),fy(8),u(8),v(8),inod(8)

common /prop/e,rmi

c data g,w,ngaus/-.77459666,0.,.77459666,

c /.5555555,.8888888,.555555,3/

data g,w,ngaus/-.57735,.57735,

/1.,1.,2/

call get(ie,nod,x,y,u,v,inod)

c Type of element according to number of nodes of element

if(nod.eq.3)then

c triangular element (S2 is twice the triangle area), see Eq.(3.1.9)

S2=x(2)*y(3)-x(3)*y(2)-x(1)*y(3)+x(3)*y(1)+x(1)*y(2)-x(2)*y(1)

E1=S2/2.*E/(1-rmi**2)

do i=1,3

i1=mod(i,3)+1

i2=mod(i1,3)+1

fx(i)=(y(i1)-y(i2))/S2

fy(i)=(x(i2)-x(i1))/S2

enddo

do i=1,3

do j=1,3

c stiffness matrix according to Eq.(3.1.5)

a(2*i-1,2*j-1)=E1*(fx(i)*fx(j)+(1-rmi)/2.*fy(i)*fy(j))

a(2*i,2*j)=E1*(fy(i)*fy(j)+(1-rmi)/2.*fx(i)*fx(j))

a(2*i,2*j-1)=E1*(rmi*fy(i)*fx(j)+(1-rmi)/2.*fx(i)*fy(j))

a(2*i-1,2*j)=E1*(rmi*fx(i)*fy(j)+(1-rmi)/2.*fy(i)*fx(j))

enddo

elseif(nod.eq.4.or.nod.eq.8)then

c Isoparametric 4 or 8-node element. Gaussian integration 2 x 2

do i=1,nl

do j=1,nl

a(i,j)=0.

enddo

b(i)=0

enddo

do ig=1,ngaus

do jg=1,ngaus

call fdf(nod,x,y,g(ig),g(jg),f,fx,fy,det)

E1=det*w(ig)*w(jg)*E/(1-rmi**2)

do i=1,nod

i1=2*i-1

i2=2*i

do j=1,nod

j1=2*j-1

j2=2*j

a(i1,j1)=a(i1,j1)+E1*(fx(i)*fx(j)+(1-rmi)/2.*fy(i)*fy(j))

a(i2,j2)=a(i2,j2)+E1*(fy(i)*fy(j)+(1-rmi)/2.*fx(i)*fx(j))

a(i2,j1)=a(i2,j1)+E1*(rmi*fy(i)*fx(j)+(1-rmi)/2.*fx(i)*fy(j))

a(i1,j2)=a(i1,j2)+E1*(rmi*fx(i)*fy(j)+(1-rmi)/2.*fy(i)*fx(j))

enddo

endif

end

subroutine fdf(nod,x,y,rksi,rny,f,fx,fy,det)

c The subroutine calculates values [f] and first derivatives [fx,fy]

c of all shape functions (no. of shape functions = no. of nodal points)

c in point rksi,rny.

dimension x(*),y(*),f(*),fx(*),fy(*),dfksi(8),dfny(8),gx(8),gy(8)

data gx,gy/-1.,1.,1.,-1.,0.,1.,0.,-1.,-1.,-1.,1.,1.,-1.,0.,1.,0./

dxksi=0

dxny =0

dyksi=0

dyny =0

do 1 i=1,nod

xx=gx(i)*rksi

yy=gy(i)*rny

if(i.le.4.and.nod.eq.4)then

f(i)=(1+yy)*(1+xx)/4

dfksi(i)=gx(i)/4*(1+yy)

dfny (i)=gy(i)/4*(1+xx)

elseif(i.le.4)then

c Shape functions according to Eq.(3.1.11)

f(i)=(1+yy)*(1+xx)*(xx+yy-1)/4

dfksi(i)=gx(i)/4*(1+yy)*(2*xx+yy)

dfny (i)=gy(i)/4*(1+xx)*(2*yy+xx)

elseif(i.eq.5.or.i.eq.7)then

f(i)=(1-rksi**2)*(1+yy)/2

dfksi(i)=-rksi*(1+yy)

dfny(i)=gy(i)*(1-rksi**2)/2

else

f(i)=(1-rny**2)*(1+xx)/2

dfksi(i)=gx(i)*(1-rny**2)/2

dfny(i)=-rny*(1+xx)

endif

dxksi=dxksi+x(i)*dfksi(i)

dxny =dxny +x(i)*dfny (i)

dyksi=dyksi+y(i)*dfksi(i)

1dyny =dyny +y(i)*dfny (i)

c Jacobian, see Eq.(3.1.13)

det = dxksi*dyny-dxny*dyksi

c Transformed derivatives, Eq.(3.1.14)

do 2 i=1,nod

fx(i)=(dyny*dfksi(i)-dyksi*dfny(i))/det

2 fy(i)=(dxksi*dfny(i)-dxny*dfksi(i))/det

end

subroutine get(ie,nod,x,y,u,v,inod)

c Extract data for element IE from COMMON area. Result: NOD-number of nodes, X,Y-coordinates,

c U,V-displacements, INOD-indices of nodes.

common z(20000)

dimension x(8),y(8),u(8),v(8),inod(8),iz(20000)

equivalence (z,iz)

c MLOCR is a library routine, identifying properties of element in the row IE and column 1

c of the connectivity matrix (number 2): LOC is position in COMMON, IT type, NOD-length of row

c (i.e. number of nodes of element IE), and IDIF is a "distance" between neighbouring elements

c in the matrix row.

call mlocr(2,ie,1,loc,it,nod,idif)

do i=1,nod

indu=iabs(iz(loc+(i-1)*idif))

inod(i)=indu

call mlocr(5,indu,1,locp,it,npu,idpu)

u(i)=z(locp)

v(i)=z(locp+idpu)

call mlocr(3,indu,1,locp,it,npu,idpu)

x(i)=z(locp)

y(i)=z(locp+idpu)

enddo

end

The following Fig. is an example of application of the three types of finite elements simultaneously: triangles 1,2,3,4, quadrilaterals with 4 nodes 5,6, and one 8-node element number 7. Graphics implemented in XMKPE1 draws actual shape of isoparametric elements with sides defined by transformation (3.1.12).

3.2.XMKPE3 – flow in a pipeline

The program calculates pressures at nodes of pipeline networks. Network is defined as a system of pipes (finite elements with two nodes). Each element, pipe, is characterised by volumetric flow-rate and radius of pipe.

Flow-rate in a pipe connecting nodes i and j is related to the pressure difference

(3.2.1)

Sum of oriented flow-rates in any node must be zero – this is the continuity equation written for example for node i

(3.2.2)

Overall system is described by equations

(3.2.3)

where Q1, Q2, are non-zero only if a given node is a source, e.g. a pump supplying liquid into the pipeline. The global matrix is assembly of the local element matrices

(3.2.4)

where

,for laminar flow Re<2300(3.2.5)

, for 2300<Re<105 (Blasius).(3.2.6)

The problem would be simple for laminar flows because then the coefficients are constant, and the system is linear. In the program XMKPE3 the coefficients ij are calculated from (3.2.5) or (3.2.6) using values of pressures from previous iteration. The decision whether the flow is laminar or turbulent is based upon Re calculated from assumption that the flow is laminar,

.(3.2.7)

The nonlinear system of equations is solved iteratively, using the so called flip-flap technique for storing results of iterations.

c MKPE3

c Vypocet tlaku v potrubni siti. Laminarni/turbulentni proudeni.

c Reseni soustavy nelinearnich rovnic metodou substituci.

c Spolecne parametry pro vsechny elementy: rho-hustota, mi-viskozita

c Vyznamy poli:

c 2 - noe [indexy uzlu elementu, matice ma dva sloupce pro 2 uzly elementu]

c 3 - pe [Re,Q-objemovy prutok,r-polomer trubky]

c 4 - xy [x,y]

c 5 - ipu [indexy parametru uzlu, zaporna hodnota znamena predepsany tlak]

c 6 - rpu1 [parametry uzlu-tlaky]

c 7 - rpu [zadavane hodnoty tlaku a prutoku]

c 8 - rpu2 [parametry uzlu-tlaky]

c 9 - rpu1 & rpu2 [soucasne zobrazeni dvou iteraci]

c 10 [textove pole hlavicek matic]

c Pozn. Pole 6 a 8 se vyuzivaji jako houpacka [stridave ukladani vysledku

c v jednotlivych iteracich]. Pred kazdou iteraci je ve vektoru 7 ulo-

c zeno zadani (zadavane tlaky nebo prutoky) a ve vektoru 6 nebo 8

c je vysledek predchozi iterace (tlaky ve vsech uzlech), ktere potre-

c bujeme pro odhad Re.

common /cradle/ip1,ip2,rho,ami

character inkey*4

external lokal,sext

data ifl/1/

rho=1000.

ami=1.

call scomin(-1)

call sintro('xmkpe3.for')

call mtobir('XMKPE3.MAT',iend)

call mtools(0,sext)

111 call sclear

call minpr(0,' Density [kg/m^3]=',18,16,2,rho)

call minpr(0,' Viscosity [Pa.s]=',18,17,2,ami)

call mefrin(2,5)

1 ip1=7+ifl

ip2=7-ifl

ifl=-ifl

c presun vektoru zadavanych tlaku do vektoru IP1

call move(ip1,7)

call mefron(2,5,ip1,50,2,lokal,0)

call graf(kbd)

c F1-zoom F2-zoomout F3-tools F4-list F5-visc [Esc]-Exit

if(kbd.eq.61)then

call mtools(0,sext)

elseif(kbd.eq.62)then

call sintro('xmkpe3.for')

elseif(kbd.eq.63)then

goto111

elseif(kbd.eq.27)then

call scwrit(16,10,28,23,' Do You want to quit? (Y/N) ')

call skbdr(inkey)

if(inkey.eq.'Y '.or.inkey.eq.'y ') goto999

endif

goto 1

999 call sclear

end

subroutine lokal(ie,nl,nn,a,b)

c Lokalni matice prutoku a(2,2) a vektor prutoku v uzlech b(2)

dimension a(nl,nl),b(nl),x(2),y(2),pres(2),ibound(2)

call gett3(ie,fi,x,y,pres,ibound,Re)

a(1,1)=fi

a(2,2)=fi

a(1,2)=-fi

a(2,1)=-fi

b(1)=0.

b(2)=0.

end

subroutine gett3(ie,fi,x,y,pres,ibound,re)

common z(10000) /cradle/ip1,ip2,rho,ami

dimension iz(10000),x(2),y(2),pres(2),ibound(2)

equivalence (iz,z)

c stanoveni polomeru trubky elementu IE

call mlocr(3,ie,1,loc,it,npe,idif)

r=z(loc+2*idif)

drsnost=z(loc+3*idif)

call mlocr(2,ie,1,locu,it,nue,idifu)

iu1=iabs(iz(locu))

iu2=iabs(iz(locu+idifu))

c zjisteni delky elementu IE ze souradnic x,y uzlu 1,2

call mlocr(4,iu1,1,lx1,idum,nxy,it1)

call mlocr(4,iu2,1,lx2,idum,nxy,it2)

x(1)=z(lx1)

x(2)=z(lx2)

y(1)=z(lx1+it1)

y(2)=z(lx2+it2)

h=sqrt((x(1)-x(2))**2+(y(1)-y(2))**2)

c indexy parametru uzlu ibound(i) elementu IE

c a hodnoty tlaku pres(i)

call mloc(5,iu1,1,lu1,it)

call mloc(5,iu2,1,lu2,it)

ibound(1)=iz(lu1)

ibound(2)=iz(lu2)

call mloc(ip2,iabs(iz(lu1)),1,lr1,it)

call mloc(ip2,iabs(iz(lu2)),1,lr2,it)

pres(1)=z(lr1)

pres(2)=z(lr2)

dp=abs(z(lr1)-z(lr2))

Re=r**3*dp*rho/(4*ami**2*h)

if(Re.lt.2300.)then

c Hagen Poisseuile Lamba=64/Re

fi=3.141*r**4/(8*ami*h)

else

c Blasius Lambda=0.316/Re^.25

fi=(.558/h)**(4./7.)*(2*r)**(19./7.)/

/ ((rho*dp)**(3./7.)*ami**(1./7.))

endif

c Ulozeni vypocteneho prutoku Q jako druheho parametru elementu

z(loc+idif)=fi*dp

c Ulozeni vypocteneho Re jako prvniho parametru elementu

z(loc)=2*z(loc+idif)*rho/(3.141*r*ami)

end

3.3.XMKPE5 – Heat transfer in a pipe (velocity profile prescribed)

The program solves the Fourier Kirchhoff equation written in cylindrical coordinate system, assuming steady state and prescribed velocity field:

(3.3.1)

Applying weighted residual method (w(x,r)) and Green's theorem we arrive to the weak integral formulation

(3.3.2)

The first term is zero if temperatures are specified at inlet (then w=0) or if the temperature profile is fully developed at outlet (dT/dx=0). The second term can be expressed using boundary equation of the third kind at wall, i.e.

(3.3.3)

Using approximations

(3.3.4)

the Eq(3.3.2) is reduced to the linear algebraic equation

(3.3.5)

c MKPE5 @Zitny 4/92

c Reseni Fourierovy Kirchhoffovy rovnice

c u.dT/dx+v.dT/dr= d(a.dT/dx)/dx +d(a.r.dT/dr)dr/r

c dvourozmerneho prenosu tepla v cylindr. s.s.

c Pouziti bilinearnich obdelnikovych elementu a Gaussovy integrace 2 x 2

c Vyznamy poli:

c 2 - noe [indexy uzlu elementu]

c 3 - xruv [x,r,u,v] v kazdem uzlu se zadavaji krome souradnic i rychlosti

c 4 - ipu [indexy parametru uzlu]

c 5 - rpu [parametry uzlu-teploty]

c 6 -

c 7 - help [textove pole]

common /prop/a /scom/iss(50)

external lokal,sext

character *4 inkey

equivalence (kbdr,iss(13))

a=1.

call scomin(-1)

call sintro('xmkpe5.for')

call mtobir('xmkpe5.mat',iend)

10 call sclear

call minpr(1,' Temperature conductivity a=',28,14,2,a)

1 call mtools(0,sext)

if(kbdr.eq.67)then

call mefrin(2,4)

call mefron(2,4,5,50,4,lokal,0)

endif

call grid(kbd)

if(kbd.eq.62) goto 10

if(kbd.ne.27) goto 1

call scwrit(16,2,28,23,' Do You want to Quit ? [Y/N]')

call skbdr(inkey)

if(inkey.ne.'Y '.and.inkey.ne.'y ')goto 1

999 call sclear

end

subroutine lokal(ie,nl,nn,ak,b)

dimension ak(nl,nl),b(nl),

/ r(8),x(8),u(8),v(8),g(3),w(3),f(8),fx(8),fr(8),temp(8),ibound(8)

common /prop/a

data g,w,ngaus/-.77459666,0.,.77459666,

/.5555555,.8888888,.555555,3/

call gett4(ie,nl,x,r,u,v,temp,ibound)

do 1 i=1,nl

b(i)=0.

do 1 j=1,nl

1 ak(i,j)=0.

do 2 ig=1,ngaus

do 2 jg=1,ngaus

call fdf(nl,x,r,g(ig),g(jg),f,fx,fr,det)

rr=0

uu=0

vv=0

do 20 i=1,nl

rr=rr+r(i)*f(i)

uu=uu+u(i)*f(i)

20 vv=vv+v(i)*f(i)

rr=rr*det*w(ig)*w(jg)

do 3 i=1,nl

do 3 j=1,nl

3 ak(i,j)=ak(i,j)+rr*(f(i)*(uu*fx(j)+vv*fr(j))+a*(fx(i)*fx(j)+

/fr(i)*fr(j)))

2 continue

end

subroutine fdf(nl,x,r,rksi,rny,f,fx,fr,det)

c vypocet hodnot [f] a derivaci [fx,fr] bazovych funkci v bode rksi,rny.

dimension x(*),r(*),f(*),fx(*),fr(*),dfksi(8),dfny(8),gx(8),gr(8)

data gx,gr/-1,1,1,-1,0,1,0,-1,-1,-1,1,1,-1,0,1,0/

dxksi=0

dxny =0

drksi=0

drny =0

do 1 i=1,nl

f(i)=(1+gr(i)*rny)*(1+gx(i)*rksi)/4

dfksi(i)=gx(i)/4*(1+gr(i)*rny)

dfny (i)=gr(i)/4*(1+gx(i)*rksi)

dxksi=dxksi+x(i)*dfksi(i)

dxny =dxny +x(i)*dfny (i)

drksi=drksi+r(i)*dfksi(i)

1 drny =drny +r(i)*dfny (i)

det = dxksi*drny-dxny*drksi

do 2 i=1,nl

fx(i)=(drny*dfksi(i)-drksi*dfny(i))/det

2 fr(i)=(dxksi*dfny(i)-dxny*dfksi(i))/det

end

subroutine gett4(ie,nl,x,r,u,v,temp,ibound)

c data transfer

common z(10000)

dimension x(8),r(8),u(8),v(8),temp(8),ibound(8),iz(10000)

equivalence (z,iz)

call mlocr(2,ie,1,loc,it,np,idif)

do 1 i=1,nl

indu=iabs(iz(loc+(i-1)*idif))

call mlocr(4,indu,1,locp,it,npu,idpu)

ibound(i)=iz(locp)

call mlocr(5,iabs(iz(locp)),1,locp,it,npu,idpu)

temp(i)=z(locp)

call mlocr(3,indu,1,locp,it,npu,idpu)

x(i)=z(locp)

r(i)=z(locp+idpu)

u(i)=z(locp+2*idpu)

1 v(i)=z(locp+3*idpu)

end

3.4.XMKPE9 – Rotationally symmetric shell

This program implements finite elements (2-nodes) for rotationally symmetric shells, described in textbook Schneider P., Vykutil J.: Stavba chemickych zarizeni II. VUT Brno, 1990.

Program calculates primary variables (3 DOF in a node):

ux - axial displacement

ur -radial displacement

-rotation

and derived stresses and moments

Nmmeridian force [N/m]

Nccircumferential force [N/m]

Qtransversal force in the normal direction [N/m]

Mmmeridian moment [N]

Mccircumferential moment

c MKPE7 @Zitny 4/92

c Rotacne symetricka skorepina zatizena vnitrnim pretlakem

c Pouzit dvouuzlovy prvek, popisovany v (opraveny 2 chyby v mat.tuh.):

c Schneider P., Vykutil J.: Stavba chemickych zarizeni II, Brno 1990

c Vyznamy poli:

c 2 - noe [indexy uzlu elementu]

c 3 - rx [r,x]

c 4 - ipu [indexy parametru uzlu-matice se 3-mi sloupci u,w,beta]

c 5 - rpu [parametry uzlu- u,w,beta]

c 6 - pe [parametry elementu-p/tlak/,t/tloustka/,

c N1,N2,M1,M2,Q - vnitrni sily]

c 7 - help [textove pole]

common /prop/E,rmi /scom/iss(50)

external lokal,sext

character inkey*4

equivalence (ict,iss(3)),(kbdr,iss(13))

E=2.1e11

rmi=0.3