Novembro Método dos Elementos Finitos9ªAula e 10ªAula

1

Sumário e Objectivos

Sumário: Programação do Método dos Elementos FinitosObjectivos da Aula: Apreensão da forma como se pode estruturar um programa de Elementos Finitos

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

2

Aplicações

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

3

TEXTO RECOMENDADO

Finite Element ProgramingAutores: Hinton, E. and Owen, D. R. J.Editora: Academic PressAno: 1977

O Programa descrito no texto está escrito em FORTRAN77

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

4

Estrutura do ProgramaSubrotinas Principais

Programa Principal - FEMPLANO

INPUT STIF B

PS

PB

LOAD B

PS

PB

FRONT STRE B

PS

PB

Loop número de Problemas

Loop nº de Casos de Carga

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

5

Descrição

FEMPLANO é o programa principal que recorre às subrotinas:

INPUT- Leitura de dados Geométricos e materiais

Capítulo 3 do Texto Recomendado

STIF(B,PS,PB) – Subrotina que procede ao cálculo das Matrizes de Rigidez dos Elementos.

STIFB Capítulo 4 do Texto Recomendado

STIFPS Capítulo 5 do Texto Recomendado

STIFPS Capítulo 6 do Texto Recomendado

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

6

Descrição

LOAD(B,PS,PB) – Subrotina que procede ao Cálculo dos Vectores de Forças Nodais dos Elementos.

LOADB Capítulo 4 do Texto RecomendadoLOADPS e LOADPB Capítulo 7 do Texto

RecomendadoFRONT – Procede à Assemblagem das Matrizes e Vectores

Elementares e resolve o Sistema de Equações por forma a permitir a obtenção dos deslocamentos.

Capítulo 8 do Texto Recomendado STRE(B,PS,PB) – Subrotina que procede ao cálculo das

Tensões ao nível dos Elementos.

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

7

Estrutura do ProgramaSubrotina INPUT

INPUT

CHECK1 NODEXY GAUSSQ CHECK2

ECHO

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

8

Subrotinas AuxiliaresSubrotina NODEXY

Determina as coordenadas do nó j a partir das coordenadas dos nós i e k

i j k

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

9

Subrotinas CHECK1 CHECK2 e GAUSSQ

As subrotinas CHECK1 e CHECK2 verificam se existem erros detectáveis nos dados lidos na Subrotina INPUT. Assim que um erro é detectado o programa aborta.

Os dados não lidos após a detecção de um erro são escritos recoorendo à subrotina ECHO.

A subrotina GAUSSQ estabelece as posições e os pesos associados aos pontos de integração a serem usados na integração numérica de Gauss.

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

10

Estrutura do ProgramaSubrotina STIF(B,PS,PB)

STIF (B,PS,PB)

MOD

(B, PS, PB)

SFR

(1,2)

JACOB

(1,2)

BMAT(B, PS, PB)

DBE

A subrotina STIFB constrói a matriz de rigidez para um Elemento de viga, A subrotina STIFPSa matriz de rigidez para Elementos planos e a subrotina STIFPB constrói a Matriz de Rigidez para elemento de placa de Mindlin. dV∫

e

Te V

κ = B DB

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

11

SUBROUTINE STIFPSINCLUDE 'FEMPLANO.COM'DIMENSION ESTIF(16,16)

CC*** LOOP OVER EACH ELEMENTC

DO 70 IELEM=1,NELEMLPROP=MATNO(IELEM)

CC*** EVALUATE THE COORDINATES OF THE ELEMENT NODAL POINTSC

DO 10 INODE=1,NNODELNODE=LNODS(IELEM,INODE)DO 10 IDIME=1,NDIME

10 ELCOD(IDIME,INODE)=COORD(LNODE,IDIME)CC*** EVALUATE THE D-MATRIXC

CALL MODPS(LPROP)THICK=PROPS(LPROP,3)

CC*** INITIALIZE THE ELEMENT STIFFNESS MATRIXC

DO 20 IEVAB=1,NEVABDO 20 JEVAB=1,NEVAB

20 ESTIF(IEVAB,JEVAB)=0.0KGASP=0

CC*** ENTER LOOPS FOR AREA NUMERICAL INTEGRATIONCDO 50 IGAUS=1,NGAUS

DO 50 JGAUS=1,NGAUSKGASP=KGASP+1EXISP=POSGP(IGAUS)ETASP=POSGP(JGAUS)

CC*** EVALUATE THE SHAPE FUNCTIONS,ELEMENTAL VOLUME,ETC.C

CALL SFR2(EXISP,ETASP)CALL JACOB2(IELEM,DJACB,KGASP)DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS)IF(THICK.NE.0.0) DVOLU=DVOLU*THICK

CC*** EVALUATE THE B AND DB MATRICESC

CALL BMATPSCALL DBE

CC*** CALCULATE THE ELEMENT STIFFNESSESC

DO 30 IEVAB=1,NEVABDO 30 JEVAB=IEVAB,NEVABDO 30 ISTRE=1,NSTRE

30 ESTIF(IEVAB,JEVAB)=ESTIF(IEVAB,JEVAB)+BMATX(ISTRE,IEVAB)*. DBMAT(ISTRE,JEVAB)*DVOLU

CC*** STORE THE COMPONENTS OF THE DB MATRIX FOR THE ELEMENTC

DO 40 ISTRE=1,NSTREDO 40 IEVAB=1,NEVAB

40 SMATX(ISTRE,IEVAB,KGASP)=DBMAT(ISTRE,IEVAB)50 CONTINUE

CC*** CONSTRUCT THE LOWER TRIANGLE OF THE STIFFNESS MATRIXC

DO 60 IEVAB=1,NEVABDO 60 JEVAB=1,NEVAB

60 ESTIF(JEVAB,IEVAB)=ESTIF(IEVAB,JEVAB)CC*** STORE THE STIFFNESS MATRIX,STRESS MATRIX AND SAMPLING POINTC COORDINATES FOR EACH ELEMENT ON DISC FILEC

WRITE(91) ESTIFWRITE(93) SMATX,GPCOD

70 CONTINUERETURNEND

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

12

Subrotina MOD(B, PS, PB)

( )2

12

1 0E 1 0

10 0 −ν

⎡ ⎤ν⎢ ⎥⎢ ⎥= ν

−ν ⎢ ⎥⎢ ⎥⎣ ⎦

D

A subrotina MODB(LPROP) calcula a matriz das constantes elásticas para o Elemento de Viga plano. D(1,1)=EI, D(2,2)=S (Rigidez ao corte). A subrotina MODPS(LPROP) calcula as matrizes das constantes elásticas para o Elemento 2D para estados planos de Tensão e Deformação. DMATX(I,J), I=1,3 e J=1,3

1 01

(1 ) 1 0(1 )(1 2 ) 1

1 20 02(1 )

⎡ ⎤⎢ ⎥−⎢ ⎥

− ⎢ ⎥= ⎢ ⎥+ − −⎢ ⎥

−⎢ ⎥⎢ ⎥−⎣ ⎦

D E

νν

ν νν ν ν

νν

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

13

MODPS(LPROP)

SUBROUTINE MODPS(LPROP)

INCLUDE 'FEMPLANO.COM'

YOUNG=PROPS(LPROP,1)POISS=PROPS(LPROP,2)DO 10 ISTRE=1,NSTREDO 10 JSTRE=1,NSTREDMATX(ISTRE,JSTRE)=0.0

10 CONTINUEIF(NTYPE.NE.1) GO TO 20

CC*** D MATRIX FOR PLANE STRESS

CASEC

CONST=YOUNG/(1.0-POISS*POISS)DMATX(1,1)=CONSTDMATX(2,2)=CONSTDMATX(1,2)=CONST*POISSDMATX(2,1)=CONST*POISSDMATX(3,3)=(1.0-POISS)*CONST/2.0GO TO 30

20 IF(NTYPE.NE.2) GO TO 30CC*** D MATRIX FOR PLANE STRAIN CASEC

CONST=YOUNG*(1.0-POISS)/((1.0+POISS)*(1.0-2.0*POISS))

DMATX(1,1)=CONSTDMATX(2,2)=CONSTDMATX(1,2)=CONST*POISS/(1.0-POISS)DMATX(2,1)=CONST*POISS/(1.0-POISS)DMATX(3,3)=CONST*(1.0-2.0*POISS)/(2.0*(1.0-

POISS))30 CONTINUE

RETURNEND

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

14

Subrotina MODPBDMATX(I,J), I=1,5 e J=1,5

⎡ ⎤= ⎢ ⎥⎣ ⎦

f

s

D 0D0 D

1 0D 1 0

0 0 (1 ) 2

ν⎡ ⎤⎢ ⎥= ν⎢ ⎥⎢ ⎥− ν⎣ ⎦

fD

( )1 0Ee0 12 1⎡ ⎤

= κ ⎢ ⎥+ ν ⎣ ⎦sD

A subrotina MODPB(LPROP) calcula a matriz das constantes elásticas para o Elemento o Elemento de Placa baseado na Teoria de Mindlin

3

2

EtD12(1 )

=−ν

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

15

Subrotina SFR1 e SFR2

A subrotina SFR1 determina as funções de forma e as suas derivadas para um elemento de viga isoparamétrico de 3 nós.

A subrotina SFR2 determina as funções de forma e as suas derivadas para um elemento plano (2D) isoparamétrico da família de Serindipity.

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

16

Subrotina SFR2

SUBROUTINE SFR2(S,T)CC*** CALCULATES SHAPE FUNCTIONS AND THEIR DERIVATIVES FOR 2D ELEMENTSC

INCLUDE 'FEMPLANO.COM'S2=S*2.0T2=T*2.0SS=S*STT=T*TST=S*TSST=S*S*TSTT=S*T*TST2=S*T*2.0

CC*** SHAPE FUNCTIONSC

SHAPE(1)=(-1.0+ST+SS+TT-SST-STT)/4.0SHAPE(2)=(1.0-T-SS+SST)/2.0SHAPE(3)=(-1.0-ST+SS+TT-SST+STT)/4.0SHAPE(4)=(1.0+S-TT-STT)/2.0

SHAPE(5)=(-1.0+ST+SS+TT+SST+STT)/4.0SHAPE(6)=(1.0+T-SS-SST)/2.0SHAPE(7)=(-1.0-ST+SS+TT+SST-STT)/4.0SHAPE(8)=(1.0-S-TT+STT)/2.0

CC*** SHAPE FUNCTION DERIVATIVESC

DERIV(1,1)=(T+S2-ST2-TT)/4.0DERIV(1,2)=-S+STDERIV(1,3)=(-T+S2-ST2+TT)/4.0DERIV(1,4)=(1.0-TT)/2.0DERIV(1,5)=(T+S2+ST2+TT)/4.0DERIV(1,6)=-S-STDERIV(1,7)=(-T+S2+ST2-TT)/4.0DERIV(1,8)=(-1.0+TT)/2.0DERIV(2,1)=(S+T2-SS-ST2)/4.0DERIV(2,2)=(-1.0+SS)/2.0DERIV(2,3)=(-S+T2-SS+ST2)/4.0DERIV(2,4)=-T-STDERIV(2,5)=(S+T2+SS+ST2)/4.0DERIV(2,6)=(1.0-SS)/2.0DERIV(2,7)=(-S+T2+SS-ST2)/4.0DERIV(2,8)=-T+STRETURNEND

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

17

Subrotina JACOB (1,2)

A subrotina JACOB 1 determina as coordenadas dos pontos de Gauss, o Jacobiano e as derivadas cartesianas das funções de forma para um elemento de viga.

A subrotina JACOB 2 determina as coordenadas dos pontos de Gauss, a matriz Jacobiana, a inversa da Jacobiana e as derivadas cartesianas das funções de forma para um elemento plano 2D.

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

18

SUBROUTINE JACOB2

SUBROUTINE JACOB2(IELEM,DJACB,KGASP)CC*** CALCULATES COORDINATES OF GAUSS POINTS AND THE JACOBIAN MATRIXC AND ITS DETERMINANT AND THE INVERSE FOR 2D ELEMENTSC

INCLUDE 'FEMPLANO.COM’DIMENSION XJACM(2,2),XJACI(2,2)

CC*** CALCULATE COORDINATES OF SAMPLING POINTC

DO 10 IDIME=1,NDIMEGPCOD(IDIME,KGASP)=0.0DO 10 INODE=1,NNODE

GPCOD(IDIME,KGASP)=GPCOD(IDIME,KGASP)+ELCOD(IDIME,INODE).*SHAPE(INODE)

10 CONTINUECC*** CREATE JACOBIAN MATRIX XJACMC

DO 20 IDIME=1,NDIMEDO 20 JDIME=1,NDIMEXJACM(IDIME,JDIME)=0.0DO 20 INODE=1,NNODE

XJACM(IDIME,JDIME)=XJACM(IDIME,JDIME)+DERIV(IDIME,INODE)*.ELCOD(JDIME,INODE)

20 CONTINUECC*** CALCULATE DETERMINANT AND INVERSE OF JACOBIAN MATRIXC

DJACB=XJACM(1,1)*XJACM(2,2)-XJACM(1,2)*XJACM(2,1)IF(DJACB.GT.0.0) GO TO 30STOP

30 XJACI(1,1)=XJACM(2,2)/DJACBXJACI(2,2)=XJACM(1,1)/DJACBXJACI(1,2)=-XJACM(1,2)/DJACBXJACI(2,1)=-XJACM(2,1)/DJACB

CC*** CALCULATE CARTESIAN DERIVATIVESC

DO 40 IDIME=1,NDIMEDO 40 INODE=1,NNODECARTD(IDIME,INODE)=0.0DO 40 JDIME=1,NDIMECARTD(IDIME,INODE)=CARTD(IDIME,INODE)+XJACI(IDIME,JDIME)*.DERIV(JDIME,INODE)

40 CONTINUE900 FORMAT(//36H PROGRAM HALTED IN SUBROUTINE JACOB2/11X,

.22H ZERO OR NEGATIVE AREA/10X,16H ELEMENT NUMBER ,I5)RETURNEND

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

19

Subrotinas BMAT(B, PS, PB)

A subrotina BMATB determina a matriz das deformações Bpara o elemento de viga.

A subrotina BMATPS determina a matriz das deformações B para o elemento plano 2D

A subrotina BMATPB determina a matriz das deformações B para o elemento de placa.

= = =Lu LNd Bdε

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

20

Subrotina BMATPS

SUBROUTINE BMATPS

INCLUDE 'FEMPLANO.COM'

NGASH=0DO 10 INODE=1,NNODEMGASH=NGASH+1NGASH=MGASH+1BMATX(1,MGASH)=CARTD(1,INODE)BMATX(1,NGASH)=0.0BMATX(2,MGASH)=0.0BMATX(2,NGASH)=CARTD(2,INODE)BMATX(3,MGASH)=CARTD(2,INODE)BMATX(3,NGASH)=CARTD(1,INODE)

10 CONTINUERETURNEND

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

21

Subrotina DBE

A subrotina DBE multiplica a matriz D pela matriz B

SUBROUTINE DBEINCLUDE 'FEMPLANO.COM'

CC*** CALCULATE D X BC

DO 10 ISTRE=1,NSTREDO 10 IEVAB=1,NEVABDBMAT(ISTRE,IEVAB)=0.0DO 10 JSTRE=1,NSTREDBMAT(ISTRE,IEVAB)=DBMAT(ISTRE,IEVAB)+

.DMATX(ISTRE,JSTRE)*BMATX(JSTRE,IEVAB)10 CONTINUE

RETURNEND

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

22

Subrotinas LOAD(B,PS,PB)

A subrotina LOADB determina as nodais equivalentes a a uma carga distribuída no elemento de viga.

A subrotina LOADPS determina as cargas nodais equivalentes a cargas gravíticas, a cargas distribuídas sobre os lados do elemento, a cargas devidas à temperatura e a cargas nodais no elemento.

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

23

Subrotinas LOAD(B,PS,PB)

LOAD (B,PS,PB)

SFR

(1,2)

JACOB

(1,2)

xei i

y

bf N dxdy

b⎧ ⎫

= − ⎨ ⎬⎩ ⎭

∫

Para mais detalhes ver: Finite Element ProgramingAutores: Hinton, E. and Owen, D. R. J.,Editora: Academic Press,Ano: 1977

Força Gravítica

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

24

Subrotina FRONT

A matriz de Rigidez Global e o Vector Força Global são obtidos e procede-se à resolução do sistema de equações obtido recorrendo ao método de eliminação de Gauss e substituição.

Novembro Método dos Elementos Finitos9ªAula e 10ªAula

25

Subrotinas STRE (B,PS,PB)

Determina as tensões nos pontos de GAUSS a partir dos deslocamentos nodais conhecidos e determinados por resolução do sistema de equações anteriormente obtido