of 20 /20
ELEMENTOS FINITOS Resumo Será apresentado um desenvolvimento matemático para o problema da equação diferencial  proposta básica, incluindo a estratégia para desenvolvimento do código computacional e validação. A seguir, nos demais projetos (proposta avançada e treliça 2D) serão apresentados apenas os códigos, levando em conta que a ideia para o desenvolvimento matemático da solução equação diferencial pelo método dos elementos finitos é similar. 1. Equação diferencial  proposta básica 1. Descrição do problema Será resolvido o seguinte problema:  Ou seja, uma EDO de ordem 2 da forma descrita acima, onde os dados de entrada serão os valores dos coeficientes A, B e C (números reais), através de um código onde será utilizaa a formulação fraca para solução em MEF por ‘n’ elementos 1D com função de aproximação linear ou quadrática. Posteriormen te será realizada a comparação entre o método aproximado e a s olução exata para validação da solução aproximada. Desenho elemento linear  nós 1 ,2 e 3 2. Desenvolvimento e solução para elementos lineares e quadráticos com condições de contorno especificadas 2.1  Desenvolvimento matemático  forma fraca da declaração integral A função y(x) será aproximada por uma função u , sendo assim, um resíduo R surgirá por conta da aproximação realizada: u u  Resta-se ainda esclarecer qual função aproximada, u(x), será utilizada em questão. Para solução pelo método dos elementos finitos, tal função de aproximação deverá ter o valor da função exata y(x i) nos pontos (nós) do elemento, e, entre estes, será realizada uma aproximação por uma função polinomial. Ou seja: u   Onde as funções N1, N2 e N3 são as chamadas “funções e forma” do elemento de tal maneira que: u u u u   

Elementos finitos: estudo de caso

Embed Size (px)

DESCRIPTION

elementos finitos 1. Elementos lineares.Solução de EDO por elementos finitos e estudo de caso.Comparações com solução analítica. Utilização de matlab e mathcad.

Text of Elementos finitos: estudo de caso

  • ELEMENTOS FINITOS

    Resumo

    Ser apresentado um desenvolvimento matemtico para o problema da equao

    diferencial proposta bsica, incluindo a estratgia para desenvolvimento do cdigo

    computacional e validao. A seguir, nos demais projetos (proposta avanada e trelia 2D)

    sero apresentados apenas os cdigos, levando em conta que a ideia para o

    desenvolvimento matemtico da soluo equao diferencial pelo mtodo dos elementos

    finitos similar.

    1. Equao diferencial proposta bsica

    1. Descrio do problema

    Ser resolvido o seguinte problema:

    Ou seja, uma EDO de ordem 2 da forma descrita acima, onde os dados de entrada

    sero os valores dos coeficientes A, B e C (nmeros reais), atravs de um cdigo onde ser

    utiliza a a formulao fraca para soluo em MEF por n elementos 1D com funo de

    aproximao linear ou quadrtica. Posteriormente ser realizada a comparao entre

    o mtodo aproximado e a soluo exata para validao da soluo aproximada.

    Desenho elemento linear ns 1 ,2 e 3

    2. Desenvolvimento e soluo para elementos lineares e quadrticos com condies de

    contorno especificadas

    2.1 Desenvolvimento matemtico forma fraca da declarao integral

    A funo y(x) ser aproximada por uma funo u , sendo assim, um resduo R

    surgir por conta da aproximao realizada:

    u

    u

    Resta-se ainda esclarecer qual funo aproximada, u(x), ser utilizada em questo.

    Para soluo pelo mtodo dos elementos finitos, tal funo de aproximao dever ter o

    valor da funo exata y(xi) nos pontos (ns) do elemento, e, entre estes, ser realizada

    uma aproximao por uma funo polinomial. Ou seja:

    u

    Onde as funes N1, N2 e N3 so as chamadas funes e forma do elemento de tal

    maneira que:

    u u u u

  • Assim, para que isso ocorra, necessrio que:

    { se se

    Ou seja, , e

    .

    Generalizando, ( ) { se i se i

    Para que a funo de aproximao u(x) seja linear, temos que as funes de forma

    devero ser:

    Garantindo, desta maneira, que: e . E assim, a

    funo de aproximao u(x) ser ento:

    u

    Onde e correspondem aos valores da funo exata nos ns 1 e 2 do elemento, ou seja

    y( e y( , satisfazendo o descrito em (4). Algumas simplificaes tambm podem ser

    feitas nas funes de forma, uma vez que no caso um elemento com aproximao linear

    (com apenas dois ns internos), tem-se que (comprimento do elemento).

    Utilizando-se agora do conceito dos resduos ponderados (minimizando o resduo ao

    longo do domnio), temos a integrao a ser resolvida:

    Onde w a funo peso (weighted function), que, resolvida pelo mtodo de Gallerkin, ser

    a prpria funo de forma em cada ponto. Ou seja:

    aller in u

    e su stituin o em resulta em

    u

    u

    Ou seja:

    O que resulta em duas equaes:

  • E, substituindo os valores de Ni e R em (10), chegamos forma forte da declarao

    integral:

    u

    (

    u

    u ) para i

    Sabemos que:

    u

    Ou seja, na verdade temos duas equaes (i=1 e i=2), formando, assim, um sistema de

    equaes a ser resolvido para as incgnitas e .

    A questo ainda a se apontar : a integral acima (forma forte) no gerar bons resultados

    para a nossa aproximao. O motivo, neste caso, fica claro ao analisar o primeiro termo do

    resduo:

    . Ora, se a nossa funo de aproximao linear, claro que a derivada

    segunda desta funo zero. Sendo assim, este termo desaparecer na integral, para

    qualquer valor de A, o que faz com que parte da fsica do problema seja perdido.

    Conclui-se, de antemo, que a funo de aproximao linear no gerar bons resultados

    caso seja utilizada a forma forte da declarao integral.

    Para funo de aproximao u(x) quadrtica, temos as seguintes funes de forma:

    Que sero funes quadrticas, que, quando substitudas em (2), resultaro em uma

    funo aproximada do tipo u(x)=ax+bx+c (onde a, b e c sero termos de , , , ,

    e ), e, ao mesmo tempo, satisfazendo o descrito em (3). Vale ressaltar que todas estas

    condies esto sendo resolvidas para um elemento. Algumas simplificaes tambm

    podem ser feitas na equao acima (equao 4), uma vez que (comprimento

    do elemento) e

    (ponto no centro do elemento).

    Ser, agora, resolvido a integrao (resduos ponderados averaged weighted residuals):

    Onde w a funo peso (weighted function), que, resolvida pelo mtodo de Gallerkin, ser

    a prpria funo de forma em cada ponto. Ou seja:

    aller in u

    e su stituin o em resulta em

  • u

    u

    u

    Ou seja (equao 9):

    O que resulta em trs equaes:

    E, substituindo os valores de Ni e R em (7), chegamos forma forte da declarao integral:

    u

    u

    u para i

    Desta vez com a funo u(x) sendo uma funo do segundo grau:

    u

    Desta vez temos que o termo da segunda derivada no se anula, sendo assim a fsica do

    problema conservada mesmo sob a soluo com a forma forte da declarao integral.

    Porm o que ocorre que na prtica, os resultados para aproximaes com mais de um

    elemento foram insatisfatrias, pois existe uma descontinuidade nas derivadas de ordem 1

    dos elementos. Assim, a hiptese de resduo ponderado nulo ao longo do conjunto de

    elementos perde a validade quando a condio de continuidade violada nas fronteiras

    (entre elementos). A condio de continuidade implica que para um integrando de ordem

    'n' a funo e todas as derivadas de ordem (n-1) devem ser contnuas (Akin, 1986).

    Rao (1998) define elementos de continuidade C0 como elementos que possuem

    continuidade apenas da varivel de campo (u(x), por exemplo). Os elementos com

    continuidade C1 possuem continuidade da varivel de campo u(x) e de sua primeira

    derivada du/dx. Sendo assim, um elemento cbico, com quatro graus de liberdade

    necessrio para soluo do problema.

    A partir deste ponto, fica claro que vivel utilizar uma alternativa mais simples para a

    declarao integral. Ser utilizada uma estratgia para diminuir a ordem das derivadas do

    integrando (eliminar o termo de ordem 2), e assim, um elemento com continuidade C0

    (linear ou quadrtico) ser suficiente para soluo do problema.

    A estratgia consiste na integrao por partes da forma forte da declarao integral,

    reduzindo a ordem da derivada no integrando.

  • Assim, sendo:

    u

    u

    Aplica-se a integrao por partes com:

    u

    f

    f f f

    vivel, contudo, que seja utilizada a forma fraca da declarao integral. Assim, ao utilizar

    a estratgia de integrao por partes conseguimos reduzir a ordem do operador

    diferencial:

    u

    u

    u

    u

    [ u

    ]

    Sendo:

    u

    u

    A nova declarao integral que contm apenas termos de derivadas de ordem 1, ou seja,

    apenas o requisito de continuidade C0 necessrio agora. Assim, as condies de contorno

    essenciais y(i) e y(f) so aplicadas nova declarao integral (sendo i e f os pontos que

    definem inicio e final do domnio do problema).

    O termo [

    ]

    representa as condies de contorno naturais, e sero introduzidas aps

    a soluo, caso seja necessrio.

    Em resumo, temos que, a partir deste momento, ser utilizada apenas a forma fraca da

    declarao integral, que faz com que possamos usar qualquer tipo de elemento (linear ou

    quadrtico) em EDO's de ordem superiores.

    2.1.1 - Matriz de rigidez para o elemento linear

    Resolvendo agora nosso problema com elemento linear utilizando a formulao fraca,

    temos:

    E uao

    u

    u

  • E uao

    u

    u

    Com:

    u

    u

    E, para o elemento linear, conforme visto anteriormente:

    O que resulta num sistema de equaes de duas incgnitas e . Escrevendo o sistema

    de equaes na forma matricial, denominaremos a matriz dos termos dependentes de

    matriz de rigidez (k - 2x2) e o vetor dos termos independentes de vetor de foras (f -

    2x1), assim, temos que:

    ky=f

    Ou seja:

    (

    ) (

    )

    (

    ) (

    )

    Onde L representa o comprimento do elemento ( ), a coordenada global do primeiro

    n do elemento e a do segundo n. Sendo assim, escrevemos na forma matricial as

    equaes para o elemento linear:

    [

    ] [

    ] [

    ]

    Percebe-se, ento, que o mtodo resulta num sistema de equaes com uma matriz

    simtrica.

    O grande desafio agora trata-se da passagem da matriz e vetores do elemento, descritos

    acima, para termos globais, levando em conta todos os graus de liberdade do domnio em

    questo.

    A estratgia consiste em sobrepor as matrizes dos elementos, levando em conta os ns que

    tem influncia sobre um ou mais elementos. Ou seja, para um problema com 2 elementos

    (total de 3 graus de liberdade) temos:

  • [

    ]

    [

    ] [f

    f f f

    ]

    Onde os elementos da matriz de rigidez global em vermelho representam a matriz de

    rigidez do elemento 1, e os em azul a matriz de rigidez do elemento 2. Da mesma forma, no

    vetor de foras, o ndice f11 corresponde ao primeiro ndice do vetor de fora do elemento

    1, e o ndice f21 ao segundo ndice do vetor de fora do elemento 2.

    Para o caso de serem utilizados 3 elementos (4 graus de liberdade globais), temos:

    [

    ]

    [

    ] [

    f f f f f

    f

    ]

    Sendo agora as cores vermelho representando o elemento 1, azul para o elemento

    2 e verde para o elemento 3. Assim ser a superposio das matrizes de rigidez para 'n'

    elementos. Desta vez, com 3 elementos, ser necessrio resolver as duas equaes do

    ''miolo''. Observa-se, que, mesmo sendo somente necessrio resolver as duas equaes do

    meio (2 e 3), estas ainda estaro dependentes dos termos Y1 e Y4, que so as condies de

    contorno do problema. Em termos de programao, no se pode simplesmente remover as

    linhas e colunas 1 e 4 e resolver o problema, pois ao fazer isso est sendo suposto que as

    condies de contorno so necessariamente iguais a zero (para as linhas e colunas

    poderem ser eliminadas por completo). Assim, outra estratgia de programao deve ser

    implementada para incorporar as condies de contorno no termo dos vetores

    independentes. As duas equaes do ''miolo'' ficam, portanto:

    f f

    f f

    Sendo que as incgnitas do problema so Y2 e Y3 (j que Y1 e Y4 so as condies de contorno do problema), assim, tem-se que:

    f f

    f f

    Ou seja, os termos em vermelho devem ser incorporados ao sistema de equaes matricial no singular, ou seja:

    [

    ] [

    ] [f f f f

    ]

    Se Y1 e Y4 fossem 0, bastaria remover as linhas e colunas dos termos dos extremos (condies de contorno). Sendo diferente de zero (e ser necessrio generalizar na

  • programao) deve-se incorpor-los ao termos de vetores independentes reduzido (no-singular).

    Sendo o sistema acima resolvido, temos a soluo do problema para elementos lineares.

    2.1.2 - Matriz de rigidez para o elemento quadrtico

    O desenvolvimento matemtico similar ao resolvido para o elemento linear, s que agora a funo de aproximao u(x) ser quadrtica. Assim tem-se:

    u

    Resolvendo agora nosso problema com elemento linear utilizando a formulao fraca,

    temos:

    E uao

    u

    u

    E uao

    u

    u

    E uao

    u

    u

    Com:

    u

    u

    u

    E, para o elemento quadrtico, conforme visto anteriormente:

    O que resulta num sistema de equaes de trs incgnitas e e . Escrevendo o

    sistema de equaes na forma matricial, denominaremos a matriz dos termos

    dependentes de matriz de rigidez (k - 3x3) e o vetor dos termos independentes de vetor

    de foras (f - 3x1), assim, temos da mesma forma:

    ky=f

    Ou seja:

    (

    ) (

    ) (

    ) f

    (

    ) (

    ) (

    ) f

  • (

    ) (

    ) (

    ) f

    Onde L representa o comprimento do elemento ( . Os termos independentes f1, f2 e f3

    tambm ficam em funo de coordenadas locais de cada elemento ( ,

    e ), que podem

    ser convertidas em coordenadas globais. Sendo assim, escrevemos na forma matricial as

    equaes para o elemento quadrtico:

    [

    ]

    [

    ] [f f f

    ]

    Percebe-se, ento, que o mtodo sempre resulta num sistema de equaes com uma matriz

    simtrica.

    A estratgia de sobrepor as matrizes dos elementos, levando em conta os ns que tem

    influncia sobre um ou mais elementos, para transform-las em matrizes globais a

    mesma apresentada anteriormente, s que agora com matrizes dos elementos 3x3. Da

    mesma forma, tambm observa-se, que, mesmo sendo somente necessrio resolver as

    duas equaes do miolo, ainda ser necessrio implementar uma outra alterao na

    programao para incluir as condies de contorno nos termos independentes,

    exatamente da mesma forma como descrito no item anterior, para elementos lineares.

    Resolvendo o sistema de equaes global, temos a soluo do problema para elementos

    quadrticos (ou seja, cada elemento com 3 ns internos). Utilizando 2 elementos (cada um

    com 1 grau de liberdade), teramos de resolver 5 incgnitas (5 ns sendo 3 em cada

    elemento, s que um deles comum aos dois elementos), sendo das 5, duas com valores j

    conhecidos, que so as condies de contorno.

    3 Etapa elaborao do cdigo geral em Mathcad descrio e aplicao do cdigo

    3.1 Escolha do software

    Apesar do Matlab ser mais poderoso no que se refere operaes matriciais, no

    Mathcad tambm possvel realizar operaes matriciais com grande nmero de

    elementos. Assim como no Maple e no Excel, o Mathcad possui uma linguagem onde o

    cdigo muito mais compreensvel visualmente, alm de ser possvel realizar o clculo de

    expresses de forma simblica de maneira mais simples.

    Alm disso, o arquivo do cdigo simplesmente uma planilha eletrnica com

    baixssimo custo computacional. O Matlab foi abandonado por ser um programa

    extremamente pesado e com executveis no to fceis de serem gerados, e, desta

    maneira, a forma de entrada de dados (input) fica prejudicada.

    O Mathca peca apenas na soluo e ata e EDOs que resolvida apenas

    numericamente e graficamente, mas no analiticamente. Como o projeto no envolve

    soluo analtica e EDOs mas ustamente o contrrio e pelo fato o autor ter ampla

  • experincia com o software, optou-se pelo desafio de criar o cdigo MEF no Mathcad com

    dados de input dinmicos e processamento rpido (arquivo em forma de planilha). O

    resultado foi positivo.

    3.2 Algoritmos/Programao do cdigo

    3.2.1 Viso geral do programa

    Pgina 1 dados em amarelo so as entradas do cdigo.

  • Pgina 2 -Dados de sada grfico comparando as duas solues, vetor u de soluo (y1,

    y2 e y3) e o vetor xG (ns dos elementos em coordenadas globais). Percebe-se que a

    soluo em 2 elementos lineares no foi muito satisfatria.

    3.2.2 Programao/Algoritmos

    O Mathca no possui uma uanti a e e funes uilt-in ou se a pre etermina as

    no software) que fossem necessrias para a soluo deste cdigo em questo. A

    possibilidade de criar funes atravs de estruturas de repetio e condicionais so

    infinitas, porm requerem uma maior expertise no software. A seguir as funes que o

    autor considera principais para programao realizada do cdigo e soluo do problema:

    I funo REPLACE (funo criada): repl(K,K1,r,c)

    Funo que insere uma matriz menor, K1, numa matriz maior, K. Os argumentos r e c da

    funo representam a posio em que o primeiro elemento da matriz K1 ser inserida na

    matriz K. A funo necessria no procedimento para sobrepor as matrizes dos elementos

    e gerar a matriz de rigidez global.

    Programao e exemplo da funo repl

    II Funo REMOVE COLUMNS e REMOVE ROWS (funo criada): remcols(A,c) e

    remrows(A,c).

    Funo ue remove linhas ou colunas e n ice c a matriz . O ar umento c eve ser

    um vetor. uma funo necessria para remoo das linhas e colunas relativas s

    restries con ies e contorno o pro lema para erar a matriz no-singular e

    soluo do sistema de equaes reduzido (com matriz no singular).

  • Descrio da funo e exemplo

    III Funo LINEAR SOLVE (funo interna do programa): lsolve(A,B)

    A funo interna lsolve(A,B) resolve o sistema de equaes lineares do tipo Ax=B, onde A

    a matriz dos termos dependentes B vetor dos termos independentes do sistema.

    IV Funo ODESOLVE (funo interna do programa): odesolve(x,a)

    A funo interna odesolve(x,a) resolve numericamente uma EDO dadas condies de

    contornos necessrias. Para a soluo preciso criar um solve loc no Mathca a

    maneira demonstrada abaixo.

  • 3.3 Aplicao/Exemplos e validaes

    3.3.1 3 elementos lineares com EDO:

    Ou seja, A = 1, B = -1 e C=1. Condies de contorno: Y1 = 1 e Y4 = 0 e domnio: .

    Soluo no cdigo Matlab do servidor:

    Soluo cdigo Mathcad elaborado:

  • Diferena 0%.

    3.3.2 4 elementos lineares com EDO:

    Ou seja, A = 1, B = -1 e C=1. Condies de contorno: Y1 = 2 e Y5 = 2 e domnio: .

    Soluo no cdigo Matlab do servidor:

    Obs: pequena inconsistncia encontrada no cdigo Galerkin_FormaFraca_4elem_linear.m

    do servidor, nas linhas 85 a 89, onde h u3 deve ser substitudo por u4.

  • Soluo cdigo elaborado no Mathcad:

  • 3.3.3 2 elementos quadrticos com EDO:

    Ou seja, A = 1, B = -1 e C=1. Condies de contorno: Y1 = 3 e Y5 = 2 e domnio: .

    Soluo no cdigo Matlab do servidor:

  • Soluo cdigo elaborado MathCad:

    Observa-se que a aproximao com 2 elementos quadrticos j mais eficiente do que

    com 4 elementos lineares.

    3.3.4 2 elementos quadrticos / 6 elementos lineares / 8 elementos lineares, com EDO:

    Ou seja, A = 1, B = -1 e C=1. Condies de contorno: Y1 = 2 e Y5 = 2 e domnio: .

  • 2 elem quadrticos

    6 elementos lineares (pior aproximao em relao a 2 quadrticos)

  • Soluo com 8 elementos lineares.

    2. Equao diferencial proposta avanada

    A ideia agora consiste em resolver uma equao diferencial de segunda ordem mais geral,

    da seguinte forma:

    f

    Tornando a soluo do problema mais generalizada. O arquivo da proposta avanada

    EQ1DA.xmcd, e, da mesma forma do problema anterior, tambm foi possvel obter

    resultados positivos tanto para elementos quadrticos quanto para lineares.

    Observa-se, ao final, que o elemento quadrtico muito mais poderoso para a soluo do

    MEF. Em alguns casos, verificou-se que 2 elementos quadrticos geravam aproximao

    melhor do que 8 elementos lineares.