of 25/25
Introdução ao Método dos Elementos Finitos J. A. J. Avila Departamento de Matemática e Estatística - DEMAT, UFSJ 36307 904, São João Del - Rei, MG [email protected] RESUMO Neste minicurso aprenderemos a importância do Método dos Elementos Finitos em relação aos outros métodos numéricos e sua aplicabilidade em problemas da ciência e da engenharia. Em particular, resolveremos numericamente, pelo Método dos Elementos Finitos Galerkin, a Equação de Condução do Calor 2D numa placa quadrada. Serão mostrados resultados numéricos da distribuição de temperatura durante o regime transiente até o regime estacionário. Palavras-chave: Método de Galerkin, Método dos Elementos Finitos, Equação de Condução do Calor. ERMAC 2010: I ENCONTRO REGIONAL DE MATEMÁTICA APLICADA E COMPUTACIONAL 11 - 13 de Novembro de 2010, São João del-Rei, MG; pg 65 - 89 65

Introdução ao Método dos Elementos Finitos - ufsj.edu.br§ão ao Método dos Elementos Finitos . J. A. J. Avila . Departamento de Matemática e Estatística - DEMAT, UFSJ . 36307

  • View
    215

  • Download
    1

Embed Size (px)

Text of Introdução ao Método dos Elementos Finitos - ufsj.edu.br§ão ao Método dos Elementos Finitos ....

  • Introduo ao Mtodo dos Elementos Finitos

    J. A. J. Avila

    Departamento de Matemtica e Estatstica - DEMAT, UFSJ 36307 904, So Joo Del - Rei, MG

    [email protected]

    RESUMO

    Neste minicurso aprenderemos a importncia do Mtodo dos Elementos Finitos

    em relao aos outros mtodos numricos e sua aplicabilidade em problemas

    da cincia e da engenharia. Em particular, resolveremos numericamente, pelo

    Mtodo dos Elementos Finitos Galerkin, a Equao de Conduo do Calor 2D

    numa placa quadrada. Sero mostrados resultados numricos da distribuio

    de temperatura durante o regime transiente at o regime estacionrio.

    Palavras-chave: Mtodo de Galerkin, Mtodo dos Elementos Finitos, Equao

    de Conduo do Calor.

    ERMAC 2010: I ENCONTRO REGIONAL DE MATEMTICA APLICADA E COMPUTACIONAL

    11 - 13 de Novembro de 2010, So Joo del-Rei, MG; pg 65 - 89 65

  • 1. INTRODUO

    O estudo de muitos fenmenos fsicos que acontecem na engenharia, biologia, oceanografia,

    astronomia, cosmologia, etc., usam domnios que so geometricamente muito complicados e

    difceis de desenhar, o qual dificulta sua resoluo. O Mtodo de Diferenas Finitas, num

    primeiro momento, trataria de resolver-las com o uso de transformaes conformes ou outras

    transformaes, porm, isto no fcil, pois, envolve sistemas de equaes diferenciais

    parciais elpticas. assim que uma nova tcnica potentssima chamada Mtodo dos Elementos

    Finitos nos ajuda a resolver estes tipos de problemas. O objetivo deste minicurso conhecer e

    saber aplicar o Mtodo dos Elementos Finitos Galerkin na soluo numrica de Equaes

    Diferenciais Parciais. Para compreender este mtodo resolveremos numericamente, por

    exemplo, a equao de conduo do calor 2D numa placa quadrada unitria com condio

    inicial e de contorno. Comearemos com a formulao clssica da equao de conduo do

    calor, logo passaremos formulao integral, aproximaremos esta ltima formulao pelo

    Mtodo de Galerkin e discretizaremos o domnio pelo Mtodo de Elementos Finitos.

    Resultados numricos sero apresentados.

    66

  • 2. EQUAO DE CONDUO DO CALOR EM 2D

    Sabe-se que a transferncia de calor o transporte de energia num corpo material devido s

    diferenas de temperatura, ou seja, sempre que existir uma diferena de temperatura em um

    meio ou entre meios ocorrer transferncia de calor. Este fenmeno acontece em trs

    mecanismos diferentes:

    Conduo

    Conveco

    Radiao

    Neste trabalho estudaremos a transferncia do calor por conduo, que a troca de energia

    entre as partes de um meio continuo que, estando em diferentes temperaturas, transferem

    energia trmica pela transferncia de energia cintica entre as partculas individuais ou grupo

    de partculas, no nvel atmico.

    2.1 Propriedades Fsicas

    Uma propriedade fsica do material (ou meio onde ocorre a conduo) se chama

    condutividade trmica, k , que depende da natureza do material. Neste trabalho assumiremos

    que o material isotrpico.

    A difusividade trmica do metal definida por:

    k

    c

    (1)

    onde k a condutividade trmica, a densidade e c o calor especfico. As unidades de

    medidas para k , e c so

    3W/ m , kg/m , J/ kg o ok C c C (2)

    Segundo a Eq.(1) e Eq. (2) a unidade de medida da difusividade trmica ,

    22 2

    3

    W/ m J/s mWm m

    J J sJ/ kg kg/m

    o

    o

    Ck

    c C

    (3)

    Em forma geral sabemos que k depende da temperatura. Neste trabalho o considerado

    constante. Em HOLMAN (1986) podemos encontrar valores das propriedades de alguns metais,

    veja Tabela 1.

    67

  • Tabela 1. Propriedades de alguns metais a 20oC

    Metal W/ m ok C 3kg/m

    J/ kg oc C 2 m /s

    Prata 41,6563 10

    Ouro 41,27 10

    Cobre 386,0 8,954 383,1 41,1234 10

    Alumnio 58,418 10

    Ao 52,3 10

    2.2 Formulao Diferencial

    O domnio, 2R , para a Equao de Conduo do Calor definido como,

    2, : 0 1, 0 1x y R x y (4)

    onde pode ser, por exemplo, uma placa quadrada de um certo metal. Neste trabalho

    ser uma placa de cobre de 21 [m ] .

    O contorno de definido por,

    4

    1

    ii

    (5)

    Denote-se a distribuio transiente de temperatura pela funo escalar , ,u u x y t .

    A formulao matemtica do problema de Conduo do Calor 2D e sem fonte

    1 2 4 3

    , , , 0,

    , ,0 0, ,

    100 , 500 , 0,o o

    uu x y t T

    t

    u x y x y

    u u u C u C t T

    (6)

    Segundo a forma dada da Eq. (6) chamada de Problema de Valor Inicial e de Contorno PVIC.

    A temperatura nos cantos superiores da placa deve satisfazer:

    0,1, 1,1, 100 , 0,ou t u t C t T (7)

    Na Figura 1 temos o domnio com suas respectivas condies de contorno.

    68

  • Figura 2.1. O domnio e suas condies de contorno.

    2.3 Formulao Integral

    Na primeira equao de Eq. (6) multipliquemos por ,x y e integremos em ,

    , ,u

    dA u dA x yt

    (8)

    Pela Eq. (1),

    , ,u k

    dA u dA x yt c

    (9)

    , ,u

    c dA k u dA x yt

    (10)

    , ,u u

    c dA k ds k u dA x yt

    n

    (11)

    ou equivalentemente,

    , ,u u

    c dA k u dA k ds x yt

    n

    (12)

    onde n o vetor normal curva . Eq. (12) a Formulao Integral da Equao de Conduo

    do Calor 2D.

    69

  • 3. MTODO DOS ELEMENTOS FINITOS GALERKIN

    No inteno de este minicurso abarcar todos os casos que pode ser estudado pelo Mtodo

    dos Elementos Finitos e sim elementos simples e com ordem de aproximao linear.

    Estudaremos os elementos triangulares de Lagrange, ou seja, com ns nos vrtices do

    tringulo. No MEF existem dois tipos de estudos o global e o local.

    3.1 Estudo Global

    quando se trabalha com todo o domnio do problema e com todos os ns que existem nele.

    Definamos a soluo aproximada global de u por

    1

    , , u ,nN

    i i

    i

    u x y t t x y

    (1)

    onde nN o nmero de ns da discretizao de e 1nN

    i i

    as funes bases globais. Como

    Eq. Erro! Fonte de referncia no encontrada. valida para qualquer ento valida para

    uma famlia finita j , 1, nj N , isto ,

    j j ju u

    c dA k u dA k dst

    n

    (2)

    desse modo a Eq. (2) , na verdade, um sistema de nN - equaes.

    A derivada temporal de u seria,

    1

    , ,u ,

    N

    i i

    i

    u x y tt x y

    t

    (3)

    e suas derivadas parciais,

    1 1

    , , , , , ,u , u

    N Ni i

    i i

    i i

    u x y t x y u x y t x yt t

    x x y y

    (4)

    Substituindo Eq. (1) em Eq. (2),

    1 1 1

    u u u , 1,n n nN N N

    ii i j i i j j i n

    i i i

    c t dA k t dA k t ds j N

    n

    (5)

    ou equivalentemente,

    1 1 1

    u u u , 1,n n nN N N

    ii i j i i j i j n

    i i i

    t c dA t k dA t k ds j N

    n

    70

  • 1 1 1

    u u u , 1,n n nN N N

    ii j i i j i j i n

    i i i

    c dA t k dA t k ds t j N

    n

    Assim, temos

    1 1

    u u 0, 1,n nN N

    ii j j i j j j n

    j j

    c dA t k dA k ds t i N

    n (6)

    Sejam as matrizes globais,

    n n

    ij i jN NM m c dA

    (7)

    1 1 + n n

    j ji i

    ij i jN NK k k dA k dA

    x x y y

    (8)

    2 2n n

    i

    ij jN NK k k ds

    n (9)

    onde M a matriz massa e 1 2K K K a matriz rigidez.

    Logo,

    1 2u u 0ij j ij ij jm k k (10)

    ou equivalentemente,

    u u

    u 0 0

    M K f

    (11)

    onde , 0f f x y o termo fonte. Podemos observar que Eq. (11) representa um

    sistema linear de EDO.

    3.2 Estudo Local

    Chamado, tambm, estudo Elementar e quando se trabalha com um elemento arbitrrio do

    domnio e com todos os ns que existem nele.

    Na Figura 1 temos um elemento arbitrrio e de cujos ns so,

    1 1 1 2 2 2 3 3 3, , , , ,n x y n x y n x y (12)

    71

  • Figura 1. Elemento arbitrrio e com seus respectivos ns.

    Defina-se a soluo aproximada elementar de eu por

    3

    1

    , , u ,e ei ii

    u x y t t w x y

    (13)

    onde as bases locais lineares para cada elemento e so:

    1 2 3 3 2 2 3 3 2

    2 3 1 1 3 3 1 1 3

    3 1 2 2 1 1 2 2 1

    1,

    2

    1,

    2

    1,

    2

    e

    e

    e

    w x y x y x y y y x x x yA

    w x y x y x y y y x x x yA

    w x y x y x y y y x x x yA

    (14)

    onde eA a rea de cada elemento e .

    Para cada elemento arbitrrio, e , definamos

    1 3 2 1 2 3 1 2 3 3 2

    2 1 3 2 3 1 2 3 1 1 3

    3 2 1 3 1 2 3 1 2 2 1

    a x x b y y d x y x y

    a x x b y y d x y x y

    a x x b y y d x y x y

    (15)

    O Jacobiano de um elemento arbitrrio, e , definido por

    3 2 2 3eJ a b a b (16)

    e 2e eJ A .

    Substituindo equaes (15) e (16) em Eq. (14) temos que as bases locais ficam

    1

    , , 1,2,3i i i iew x y d b x a y iJ (17)

    Com a nica finalidade de simplificar os clculos de integrais sobre e , transformaremos o

    tringulo e num triangulo e centrado na origem, tal como mostra a Figura 3.

    72

  • Figura 2. Transformao do triangulo e no triangulo e .

    Considerando esta transformao, temos

    1

    2

    3

    , 1

    ,

    ,

    w

    w

    w

    (18)

    Fazendo alguns clculos encontramos as coordenadas de rea (ou coordenadas naturais),

    1 2 3, , dadas por

    1

    , 1,2,3i i i ie d b x a y iJ (19)

    De aqui temos que as bases locais coincidem com as coordenadas de rea,

    1 1 2 3 1

    2 1 2 3 2

    3 1 2 3 3

    , ,

    , ,

    , ,

    w

    w

    w

    (20)

    Desse modo, o clculo da integral de uma funo ,f f x y sobre um elemento qualquer,

    e , seria

    1 2 3 1 2 , , , , e e

    e e ef x y dxdy J f d d J f d d (21)

    As derivadas das coordenadas de rea seriam

    1 1 1 1

    2 2 2 2

    3 3 3 3

    e e

    e e

    e e

    b a

    x J y J

    b a

    x J y J

    b a

    x J y J

    (22)

    73

  • e para calcular a integral de ,f x y

    x

    sobre e , temos

    1 2 31 2

    1 2 3

    1 2 3 1 21 2 3

    +

    +

    e

    e e

    e

    f f f fdxdy J d d

    x x x x

    f f fb b b d d

    (23)

    De forma anloga para a seguinte integral

    1 2 3 1 21 2 3

    + e e

    f f f fdxdy a a a d d

    y

    (24)

    A frmula para calcular a integral sobre um elemento e dada por

    1 2 3 1 2 0

    ! ! !, , ,

    2 !

    p q r

    e

    p q rd d p q r

    p q r

    (25)

    Seja ei um segmento do contorno de e . Definimos o comprimento de

    e

    i por iL ,

    2 2 2 2 2 2

    1 3 3 2 1 1 3 2 2, , L a b L a b L a b (26)

    A frmula para calcular a integral de linha sobre o segmento ei dada por

    1 2 1 2 0

    ! ! , ,

    1 !e ei i

    p q p q

    i i

    p qw w ds L ds L p q

    p q

    (27)

    Procedendo de forma anloga, como no caso global, chegamos ao seguinte sistema linear de

    EDO

    u u

    u 0 0

    e e e e e

    e

    M K f

    (28)

    onde 0ef o termo fonte.

    3.3 Clculo das Matrizes e Vetores Elementares

    Agora calcularemos as matrizes elementares: a matriz massa e a matriz de rigidez e os vetores

    elementares: vetor carga.

    3.3.1. Matriz Massa

    3 3

    e e

    ij i je

    M m c w w dA

    (29)

    74

  • 1 2 e e

    i j i j i je e e

    c w w dxdy cJ w w d d cJ d d

    Ento,

    3 3

    e e e

    ij ijM m cJ I (30)

    onde 1 2 ij i j

    eI d d calculada pela formula dada em Eq. (25), isto ,

    1 2

    1se

    12

    1se

    24

    ij i je

    i j

    I d d

    i j

    (31)

    3.3.2. Matriz Rigidez

    A Matriz Rigidez elementar dada por

    1 2e e eK K K (32)

    Clculo de 1eK

    1 13 3

    e eij xi xj yi yje

    K k k w w w w dA

    (33)

    ento

    1 1 2 eij xi xj yi yj xi xj yi yj ij ije e e

    k k w w w w dA k w w dA w w dA k I I (34)

    onde,

    1 ij xi xje

    I w w dA (35)

    e

    2 ij yi yje

    I w w dA (36)

    Calculando 1ijI temos que

    11 2 3 1 2 3

    1 2 3 1 2 3

    1

    j j ji i iij xi xj ee e

    I w w dA b b b b b b dAJ

    (37)

    1

    1

    jik se e

    k s

    k ik s jse e

    b b dAJ

    b b dAJ

    onde ij o delta de Kronecker. Observe que para k i e s j , temos

    75

  • 1

    1 1 1

    2

    e

    ij i j i j i je e eeI bb dA bb A bb

    J J J (38)

    Por tanto,

    11

    2ij i je

    I bbJ

    (39)

    De forma anloga temos para

    2

    1 1 1

    2

    e

    ij i j i j i je e eeI a a dA a a A a a

    J J J (40)

    ou seja,

    21

    2ij i je

    I a aJ

    (41)

    Substituindo equaes (39) e (41) em Eq. (34),

    1 13 3

    1 1

    2 2 2

    e e

    ij i j i j i j i je e e

    kK k k bb a a bb a a

    J J J

    (42)

    Clculo de 2eK

    Com ajuda da Eq. (12), calculemos a normal exterior em cada lado do elemento e , veja Figura

    4:

    1 2 1 2 1 2 1 3 3 1 3 3

    2 3 2 3 2 3 2 1 1 2 1 1

    3 1 3 1 3 1 3 2 2 3 2 2

    , , n ,

    , , n ,

    , , n ,

    v n n x x y y a b b a

    v n n x x y y a b b a

    v n n x x y y a b b a

    (43)

    Como as normais que precisamos devem ser unitrias, ento

    3 311 121 1

    n , n ,n n

    b a (44)

    1 121 222 2

    n , n ,n n

    b a (45)

    2 231 323 3

    n , n ,n n

    b a (46)

    76

  • Figura 3. Fronteira de um elemento arbitrrio.

    Seja n o vetor normal curva 3

    1

    e e

    i

    i

    . Ento

    2 2 n

    e

    e e i

    ij j ijN N

    wK k k w ds kII

    (47)

    onde

    1 2 3

    1 2 3

    nn

    n n n

    e e

    e e e

    i

    ij j i j

    i j i j i j

    wII w ds w w ds

    w w ds w w ds w w ds

    (48)

    Para calcular a matriz ijII devemos saber qual o lado do elemento que faz parte da fronteira

    do domnio .

    Suponhamos que 1e , ento

    1

    11 12n + ne

    aij

    ij xi yi j

    I

    II w w w ds

    Calculo de a

    ijI

    1 1

    1

    111 12 11 12

    1 111 12 11 12

    111 12

    1 11 1 12 1 11 1 12

    12 11 2 12 2 11 2 12

    n + n n + n

    1!n n n n

    1 1 !

    n n2

    n n n n 0

    n n n n 02

    0 0 0

    e e

    e

    a

    ij xi yi j i i je

    i i j i ie e

    i ie

    e

    LI w w w ds b a ds

    J

    L Lb a ds b a

    J J

    Lb a

    J

    b a b aL

    b a b aJ

    77

  • Suponhamos que 2

    e , ento

    2

    21 22n + ne

    bij

    ij xi yi j

    I

    II w w w ds

    Calculo de b

    ijI

    2

    221 22 21 22

    22 21 2 22 2 21 2 22

    3 21 3 22 3 21 3 22

    n + n n n2

    0 0 0

    0 n n n n2

    0 n n n n

    e

    b

    ij xi yi j i ie

    e

    LI w w w ds b a

    J

    Lb a b a

    Jb a b a

    Suponhamos que 3e , ento

    3

    31 32n + ne

    cij

    ij xi yi j

    I

    II w w w ds

    Calculo de c

    ijI

    2

    3

    31 32 31 32

    1 31 1 32 1 31 1 32

    3

    3 31 3 32 3 31 3 32

    n + n n n2

    n n 0 n n

    0 0 02

    n n 0 n n

    e

    c

    ij xi yi j i ie

    e

    LI w w w ds b a

    J

    b a b aL

    Jb a b a

    3.4 Domnio Discretizado

    Na Figura 4, observamos uma malha para . Esta malha nos fornece os seguintes dados:

    12, 11, 3, 8e n i c n iN N N N N N (49)

    Figura 4. Malha para o domnio.

    78

  • onde eN o Nmero de elementos, nN o Nmero de ns, iN o Nmero de incgnitas e

    cN o Nmero de ns no contorno.

    Por exemplo, para o elemento 1 , como mostra a Figura 5,

    Figura 5. Elemento nmero 1 da Malha.

    temos que para 9i , 3j e 7k definimos a seguintes matriz elementar e o vetor

    derivada elementar, respectivamente,

    (1) (1) (1) (1)

    99 93 97 9

    (1) (1) (1) (1) (1) (1) (1)

    39 33 37 33 3(1) (1) (1) (1)

    79 73 77 7

    u

    , u u

    u

    ij

    m m m

    m m m m m

    m m m

    (50)

    Na Fig. 5, a escolha do n i coincidir com o n 9 feita pelo software gerador de malha.

    De forma anloga, temos para

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

    99 93 97 99 93 97

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

    39 33 37 39 33 373 3 3 31(1) 1(1) 1(1) 2(1) 2(1) 2(1)

    79 73 77 79 73 77

    , ij ij

    k k k k k k

    k k k k k k k k k k

    k k k k k k

    (51)

    e

    (1) (1)

    9 9

    (1) (1) (1) (1)

    3 3

    (1) (1)

    7 7

    u

    u u ,

    u

    f

    f f

    f

    (52)

    Do mesmo modo faz para todos os elementos da malha.

    3.5 Ensamblagem

    Consiste em montar todas as eN matrizes elementares numa nica matriz chamada de matriz

    global cuja ordem n nN N .

    3.5.1 Montagem

    Uma maneira de visualizar todas as matrizes elementares numa nica matriz como

    mostramos a continuao

    79

  • 1(1)

    2(2)

    ( )

    1 11(1) 2(1)

    2 21(2) 2(2)

    1( ) 2( )

    u0 0

    0 0 u

    0 0 u

    u0 0 0 0

    0 0 0 0 u

    0 0 0 0 u

    ee

    e ee e

    N N

    N N N N

    m

    m

    m

    fk k

    k k f

    k k f

    (53)

    3.5.2 Obteno da Matriz e Vetor Global

    O gerador de malha nos fornece um arquivo de dados, assim como mostra a Tabela 1.

    Tabela 1

    Elemento i j k

    1 9 3 7

    2 6 3 9

    3 5 2 9

    4 9 2 6

    5 10 1 5

    6 8 1 10

    7 5 9 10

    8 7 4 11

    9 11 4 8

    10 11 9 7

    11 11 8 10

    12 9 11 10

    Como o 12eN , ento devemos elaborar 12 tabelas cada uma com sua respectiva matriz

    elementar. A seguir estas tabelas,

    1 1 2 3 4 5 6 7 8 9 10 11 2 1 2 3 4 5 6 7 8 9 10 11

    1 1

    2 2

    3 3

    4 4

    5 5

    6 6

    7 7

    8 8

    9 9

    10 10

    11 11

    80

  • 3 1 2 3 4 5 6 7 8 9 10 11 4 1 2 3 4 5 6 7 8 9 10 11

    1 1

    2 2

    3 3

    4 4

    5 5

    6 6

    7 7

    8 8

    9 9

    10 10

    11 11

    5 1 2 3 4 5 6 7 8 9 10 11 6 1 2 3 4 5 6 7 8 9 10 11

    1 1

    2 2

    3 3

    4 4

    5 5

    6 6

    7 7

    8 8

    9 9

    10 10

    11 11

    7 1 2 3 4 5 6 7 8 9 10 11 8 1 2 3 4 5 6 7 8 9 10 11

    1 1

    2 2

    3 3

    4 4

    5 5

    6 6

    7 7

    8 8

    9 9

    10 10

    11 11

    9 1 2 3 4 5 6 7 8 9 10 11 10 1 2 3 4 5 6 7 8 9 10 11

    1 1

    2 2

    3 3

    4 4

    5 5

    6 6

    7 7

    8 8

    9 9

    10 10

    11 11

    11 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11

    1 1

    2 2

    3 3

    4 4

    5 5

    6 6

    7 7

    8 8

    9 9

    10 10

    11 11

    81

  • Observe que a cor escura, em cada tabela, representa as componentes das matrizes

    elementares. O passo seguinte imaginar-se que estas tabelas representam matrizes de

    ordem 11 11 . Compreendido isso, agora s temos que somar estas 12 matrizes, resultando

    em uma nica matriz (ou tabela). Esta matriz chamada matriz global.

    Tabela 2. Matriz global de ordem 11 11 1 2 3 4 5 6 7 8 9 10 11

    1 511+611 0 0 0 515 0 0 618 0 51,10+61,10 0

    2 0 322+422 0 0 325 426 0 0 329+429 0 0

    3 0 0 133+233 0 0 236 137 0 139+239 0 0

    4 0 0 0 844+944 0 0 847 948 0 0 84,11+94,11

    5 551 352 0 0 355+555 +755

    0 0 0 359+759 55,10+75,10 0

    6 0 462 263 0 0 266+466 0 0 269+469 0 0

    7 0 0 173 874 0 0 177+877 +1077

    0 179+1079 0 87,11+107,11

    8 681 0 0 984 0 0 0 688+988 +1188

    0 68,10+118,10 98,11+118,11

    9 0 392+492 193+293 0 395+795 296+496 197+1097 0 199+299 +399+499 +799+1099 +1299

    79,10+129,10 109,11+129,11

    10 510,1+610,1 0 0 0 510,5 +710,5

    0 0 610,8 +1110,8

    710,9+ 1210,9

    510,10+ 610,10+ 710,10+ 1110,10+ 1210,10

    1110,11 +1210,11

    11 0 0 0 811,4+ 911,4

    0 0 811,7+ 1011,7

    911,8+ 1111,8

    1011,9+ 1211,9

    1111,10 1211,10

    811,11+911,11 +1011,11+ 1111,11+1211,11

    A estrutura da matriz global tem a seguinte forma

    1 2 3 4 5 6 7 8 9 10 11

    1

    2

    3

    4

    5

    6

    7 8 9 10 11

    onde os quadradinhos brancos representam os zeros da matriz. De forma anloga, devemos

    elaborar 12 colunas cada uma com seu respectivo vetor elementar. A seguir estas colunas,

    1 2 3 4 5 6 7 8 9 10 11 12

    1

    2

    3

    4

    5

    6

    7 8 9 10 11

    82

  • Tabela 3. Vetor global de ordem 11 1

    1 51+61 2 32+42

    3 13+23

    4 84+94

    5 35+55+75 6 26+46

    7 17+87+107

    8 68+98+118

    9 19+29+39+79+109+129

    10 510+610+710+1110+1210

    11 811+9110+1011+1111+1211

    Agora expressemos em forma matricial

    11 15 18 1 10

    22 25 26 29

    33 36 37 39

    44 47 48 4 11

    51 52 55 59 5 10

    62 63 66 69

    73 74 77 79 7 10

    81 84 88 8 10 8 11

    92 93 95 96 97 99 9 10 9 1

    0 0 0 0 0 0 0

    0 0 0 0 0 0 0

    0 0 0 0 0 0 0

    0 0 0 0 0 0 0

    0 0 0 0 0 0

    0 0 0 0 0 0 0

    0 0 0 0 0 0

    0 0 0 0 0 0

    0 0 0

    m m m m

    m m m m

    m m m m

    m m m m

    m m m m m

    m m m m

    m m m m m

    m m m m m

    m m m m m m m m

    1

    2

    3

    4

    5

    6

    7

    8

    1 9

    10 1 10 5 10 8 10 9 10 10 10 11 10

    11 4 11 7 11 8 11 9 11 10 11 11 11

    u

    u

    u

    u

    u

    u

    u

    u

    u

    0 0 0 0 0 u

    0 0 0 0 0 u

    m m m m m m

    m m m m m m

    +

    11 15 18 1 10

    22 25 26 29

    33 36 37 39

    44 47 48 4 11

    51 52 55 59 5 10

    62 63 66 69

    73 74 77 79 7 10

    81 84 88 8 10 8 11

    92 93 95 96 97 99 9 10 9 1

    0 0 0 0 0 0 0

    0 0 0 0 0 0 0

    0 0 0 0 0 0 0

    0 0 0 0 0 0 0

    0 0 0 0 0 0

    0 0 0 0 0 0 0

    0 0 0 0 0 0

    0 0 0 0 0 0

    0 0 0

    k k k k

    k k k k

    k k k k

    k k k k

    k k k k k

    k k k k

    k k k k k

    k k k k k

    k k k k k k k k

    1 1

    2 2

    3 3

    4 4

    5 5

    6 6

    7 7

    8 8

    1 9 9

    10 1 10 5 10 8 10 9 10 10 10 11 10 10

    11 4 11 7 11 8 11 9 11 10 11 11 11 11

    u

    u

    u

    u

    u

    u

    u

    u

    u

    0 0 0 0 0 u

    0 0 0 0 0 u

    f

    f

    f

    f

    f

    f

    f

    f

    f

    k k k k k k f

    k k k k k k f

    1 2u ug g g g g gM K K f (54)

    83

  • u ug g g g gM K f (55)

    onde 0gf e 1 2g g gK K K .

    3.5.3 Obteno da Matriz e Vetor a Resolver

    Anlise na matriz M . Usando as condies de contorno

    11m 0 0 0 15m 0 0 18m 0 1 10m 0

    0 22m 0 0 25m 26m 0 0 29m 0 0

    0 0 33m 0 0 36m 37m 0 39m 0 0

    0 0 0 44m 0 0 47m 48m 0 0 4 11m

    51m 52m 0 0 55m 0 0 0 59m 5 10m 0

    0 62m 63m 0 0 66m 0 0 69m 0 0

    0 0 73m 74m 0 0 77m 0 79m 0 7 10m

    81m 0 0 84m 0 0 0 88m 0 8 10m 8 11m

    92 93 95 96 97

    10 1 10 5 1

    1

    99 9 10 9 11

    10 9 10 10 10 11

    11 9 11 1

    0 8

    11 4 11 7 11 8 0 11 11

    0 0 0

    0 0 0 0 0

    0

    u

    0 0 0 0

    m m m m m

    m m m

    m

    m m m

    m m

    m m m

    m m m

    2u

    3u

    4u

    5u

    6u

    7u

    8u

    9

    10

    11

    u

    u

    u

    De forma anloga para a matriz K . Desse modo temos

    99 9 10 9 11 9 99 9 10 9 11 9 9

    10 9 10 10 10 11 10 10 9 10 10 10 11 10 10

    11 9 11 10 11 11 11 11 9 11 10 11 11 11 11

    u u

    u u

    u u

    i iN N

    m m m k k k f

    m m m k k k f

    m m m k k k f

    92 93 95 96 97

    10 1 10 5 10 8

    11 4 11 7 11 8

    92 93 95

    1

    2

    3

    4

    5

    6

    7

    9

    8

    u

    u

    u

    u

    u

    u

    u

    u

    0 0 0

    0 0 0 0 0

    0 0 0 0 0

    0 0

    i cM N N

    m m m m m

    m m m

    m m m

    k k k k

    6 97

    10 1 10 5 10 8

    11 4 11 7 11 8

    1

    2

    3

    4

    5

    6

    7

    8

    u

    0

    0 0 0 0 0

    0

    u

    u

    u

    u

    u

    u

    0 0

    u

    0 0

    i cK N N

    k

    k k k

    k k k

    ou equivalentemente

    84

  • u u u uM K f M K (56)

    onde 0f . Assim,

    u u fM K (57)

    onde f u uf M K

    3.6 Soluo do Sistema de EDO

    Malha temporal

    O tamanho de passo temporal definido por:

    T

    tN

    (58)

    e os ns temporais, por

    0 , 1,nt t n t n N (59)

    onde 0 0t o tempo inicial.

    Condio Inicial

    0 0u t (60)

    Regra do Mdio ponto Generalizado

    De Eq. (60), temos que

    0 0u (61)

    Faa para 1,n N

    1 1 1 u f 1 fu n nn nM t K M t K t (62)

    Observe que Eq. (62) representa um Sistema de Equaes Lineares, veja LEWIS et. al (2004).

    Ax b (63)

    onde

    1 1u

    1 u f 1 f

    n

    n n n

    A M t K

    x

    b M t K t

    (64)

    Para resolver o Sistema Linear em Eq. (63) usamos o mtodo SOR. Lembrando que 0 1 .

    Neste trabalho usamos 0,5 .

    85

  • 4. RESULTADOS NUMRICOS

    O cdigo computacional foi implementado na Linguagem Fortran (Compaq Visual Fortran). Os

    programas foram rodados num laptop T6500 Processador Intel CoreTM 2 Duo. Para a gerao

    da malha se uso o software (livre) Gmsh 2.5.0 e para ps-processamento se uso o software

    (comercial) Tecplot 360.

    4.1 Malha Estruturada e no Estruturada

    As Figuras 4.1 e 4.2 mostram dois tipos de malhas que foram geradas com o software Gmsh,

    uma Malha Estruturada e uma Malha no Estruturada.

    Figura 4.1 Malha Estruturada. Figura 4.2 Malha no Estruturada

    A Tabela 4.1 mostra o Nmero de elemento e o Nmero de ns das duas malhas.

    Tabela 4.1 Elementos e ns das malhas

    Tipo de Malha eN nN

    Estruturada 1024 545

    No Estruturada 824 449

    4.2 Resultados Numricos

    O Nmero de incgnitas que resultou da malha estruturada foi 481iN e da malha no

    estruturada foi 377iN . O tempo CPU aps atingir o estado permanente mostrado na

    Tabela 4.2.

    Tabela 4.2 Nmero de incgnitas e tempo CPU.

    Tipo de Malha iN Tempo CPU em [s]

    Estruturada 481 22

    No Estruturada 377 14

    86

  • As Figuras 4.3 e 4.4 mostram a convergncia ao longo do tempo da distribuio de

    temperatura em ambas as malhas. Usou-se um 35 10t e se atingiu o estado estacionrio

    em 2 [s] com 400N . Tambm, pode-se observar que a distribuio de temperatura na

    placa quadrada comea a diminuir uniformemente desde a parte superior at um pouco mais

    da metade da placa em forma de elipses concntricas. Por exemplo, o valor da temperatura no

    centro da placa, isto , em 0.5,0.5 de 195,84[ ]oC .

    Figura 4.3 Distribuio de temperatura na malha estruturada.

    Figura 4.4 Distribuio de temperatura na malha no estruturada.

    A seguir mostraremos alguns quadros, em instantes de tempo diferentes, do estado transiente

    da conduo do calor at atingir o estado estacionrio.

    0 [ ]t s

    87

  • 0,5 [ ]t s

    1 [ ]t s

    2 [ ]t s

    Figura 4.5 Distribuio de temperatura em alguns instantes de tempo, desde o tempo inicial

    at o tempo em que atingiu o estado estacionrio.

    Desse modo conclumos a simulao numrica da conduo do calor numa placa quadrada.

    88

  • 5. REFERNCIAS BIBLIOGRFICAS

    HOLMAN J. P. Heat Transfer. 6a Edio. McGraw-Hill Book Company, New York, 1986.

    LEWIS Roland W, NITHIARASU Perumal e SEETHARAMU, Kankanhally N. Fundamentals of the

    Finite Element Method for Heat and Fluid Flow. John Wiley and Sons, 2004

    89