7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 1/20

ELEMENTOS FINITOS

Resumo

Será apresentado um desenvolvimento matemático para o problema da equaçãodiferencial

–  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 odesenvolvimento matemático da solução equação diferencial pelo método dos elementosfinitos é 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 entradaserã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. Posteriormente será realizada a comparação entreo método aproximado e a soluçã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 decontorno 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 Rsurgirá 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 ovalor da função exata y(xi) nos pontos (nós) do elemento, e, entre estes, será realizadauma 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 talmaneira que:

u u u u    

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 2/20

 

Assim, para que isso ocorra, é necessário que:

{ se

se    

Ou seja, , e .

Generalizando, () { se i se i  Para que a função de aproximação u(x) seja linear, temos que as funções de forma  deverão ser:

   

Garantindo, desta maneira, que:  e . E assim, afunção de aproximação u(x) será então:

u    

Onde  e  correspondem aos valores da função exata nos nós 1 e 2 do elemento, ou sejay( e y(, satisfazendo o descrito em (4). Algumas simplificações também podem serfeitas nas funções de forma, uma vez que no caso um elemento com aproximação linear(com apenas dois nós internos), tem-se que

 (comprimento do elemento).

Utilizando-se agora do conceito dos resíduos ponderados (minimizando o resíduo aolongo do domínio), temos a integração a ser resolvida:

∫  

Onde w é a função peso (weighted function), que, resolvida pelo método de Gallerkin, seráa própria função de forma em cada ponto. Ou seja:

allerin u   e sustituino em resulta em u   u  Ou seja:

∫

 

O que resulta em duas equações:

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 3/20

∫

∫

 

E, substituindo os valores de Ni e R em (10), chegamos à forma forte da declaração

integral:

∫ u u u para i  

Sabemos que: u  

Ou seja, na verdade temos duas equações (i=1 e i=2), formando, assim, um sistema deequações a ser resolvido para as incógnitas  e .

A questão ainda a se apontar é: a integral acima (forma forte) não gerará bons resultadospara a nossa aproximação. O motivo, neste caso, fica claro ao analisar o primeiro termo do

resíduo: . Ora, se a nossa função de aproximação é linear, é claro que a derivada

segunda desta função é zero. Sendo assim, este termo desaparecerá na integral, paraqualquer valor de A, o que faz com que parte da física do problema seja perdido.

Conclui-se, de antemão, que a função de aproximação linear não gerará bons resultadoscaso seja utilizada a forma forte da declaração integral.

Para função de aproximação u(x) quadrática, temos as seguintes funções de forma:

  

Que serão funções quadráticas, que, quando substituídas em (2), resultarão em umafunção aproximada do tipo u(x)=ax²+bx+c (onde a, b e c serão termos de ,  ,, ,  e ), e, ao mesmo tempo, satisfazendo o descrito em (3). Vale ressaltar que todas estascondições estão sendo resolvidas para um elemento. Algumas simplificações tambémpodem ser feitas na equação acima (equação 4), uma vez que  (comprimento

do elemento) e

⁄  (ponto  no centro do elemento).

Será, agora, resolvido a integração (resíduos ponderados – averaged weighted residuals):

∫  

Onde w é a função peso (weighted function), que, resolvida pelo método de Gallerkin, seráa própria função de forma em cada ponto. Ou seja:

allerin

u   e sustituino em resulta em

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 4/20

u   u   u    

Ou seja (equação 9):

∫

 

O que resulta em três equações:

∫

∫

∫

 

E, substituindo os valores de Ni e R em (7), chegamos à forma forte da declaração integral:

∫ u u u para i  

Desta vez com a função u(x) sendo uma função do segundo grau:

u

 

Desta vez temos que o termo da segunda derivada não se anula, sendo assim a física doproblema é conservada mesmo sob a solução com a forma forte da declaração integral.

Porém o que ocorre é que na prática, os resultados para aproximações com mais de umelemento foram insatisfatórias, pois existe uma descontinuidade nas derivadas de ordem 1dos elementos. Assim, a hipótese de resíduo ponderado nulo ao longo do conjunto deelementos perde a validade quando a condição de continuidade é violada nas fronteiras(entre elementos). A condição de continuidade implica que para um integrando de ordem'n' a função e todas as derivadas de ordem (n-1) devem ser contínuas (Akin, 1986).

Rao (1998) define elementos de continuidade C0 como elementos que possuemcontinuidade apenas da variável de campo (u(x), por exemplo). Os elementos comcontinuidade C1 possuem continuidade da variável de campo u(x) e de sua primeiraderivada du/dx. Sendo assim, um elemento cúbico, com quatro graus de liberdade énecessário para solução do problema.

A partir deste ponto, fica claro que é viável utilizar uma alternativa mais simples para adeclaração integral. Será utilizada uma estratégia para diminuir a ordem das derivadas dointegrando (eliminar o termo de ordem 2), e assim, um elemento com continuidade C0(linear ou quadrático) será suficiente para solução do problema.

A estratégia consiste na integração por partes da forma forte da declaração integral,reduzindo a ordem da derivada no integrando.

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 5/20

Assim, sendo:

∫ u

u  

Aplica-se a integração por partes com:

u   f  

∫ f f ∫ f   É viável, contudo, que seja utilizada a forma fraca da declaração integral. Assim, ao utilizara estratégia de integração por partes conseguimos reduzir a ordem do operadordiferencial:

∫ u u  

∫ u u  

u

 

Sendo:

∫

u

u  

   

A nova declaração integral que contém apenas termos de derivadas de ordem 1, ou seja,apenas o requisito de continuidade C0 é necessário agora. Assim, as condições de contornoessenciais y(i) e y(f) são aplicadas à nova declaração integral (sendo i e f os pontos quedefinem inicio e final do domínio do problema).

O termo

 representa as condições de contorno naturais, e serão introduzidas após

a solução, caso seja necessário.

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

declaração integral, que faz com que possamos usar qualquer tipo de elemento (linear ouquadrático) 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 formulação fraca,temos:

Euação ∫

u

u

 

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 6/20

Euação ∫ u u

 

Com:

u   u  

E, para o elemento linear, conforme visto anteriormente:

 

O que resulta num sistema de equações de duas incógnitas

 e

. Escrevendo o sistema

de equações na forma matricial, denominaremos a matriz dos termos dependentes dematriz de rigidez (k - 2x2) e o vetor dos termos independentes de vetor de forças (f -2x1), assim, temos que:

y= 

Ou seja:

 

 Onde L representa o comprimento do elemento (),  é a coordenada global do primeironó do elemento e   a do segundo nó. Sendo assim, escrevemos na forma matricial asequações para o elemento linear:

 

 

Percebe-se, então, que o método resulta num sistema de equações com uma matrizsimétrica.

O grande desafio agora trata-se da passagem da matriz e vetores do elemento, descritosacima, para termos globais, levando em conta todos os graus de liberdade do domínio emquestão.

A estratégia consiste em sobrepor as matrizes dos elementos, levando em conta os nós quetem influência sobre um ou mais elementos. Ou seja, para um problema com 2 elementos(total de 3 graus de liberdade) temos:

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 7/20

ff ff  

Onde os elementos da matriz de rigidez global em vermelho representam a matriz derigidez do elemento 1, e os em azul a matriz de rigidez do elemento 2. Da mesma forma, novetor de forças, o índice f11 corresponde ao primeiro índice do vetor de força do elemento1, e o índice f21 ao segundo índice do vetor de força do elemento 2.

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

ff ff ff  

Sendo agora as cores vermelho representando o elemento 1, azul para o elemento2 e verde para o elemento 3. Assim será a superposição das matrizes de rigidez para 'n'elementos. Desta vez, com 3 elementos, será necessário resolver as duas equações do''miolo''. Observa-se, que, mesmo sendo somente necessário resolver as duas equações domeio (2 e 3), estas ainda estarão dependentes dos termos Y1 e Y4, que são as condições decontorno do problema. Em termos de programação, não se pode simplesmente remover aslinhas e colunas 1 e 4 e resolver o problema, pois ao fazer isso está sendo suposto que ascondições de contorno são necessariamente iguais a zero (para as linhas e colunaspoderem ser eliminadas por completo). Assim, outra estratégia de programação deve serimplementada para incorporar as condições de contorno no termo dos vetoresindependentes. As duas equações do ''miolo'' ficam, portanto:

ff

ffSendo que as incógnitas do problema são Y2 e Y3 (já que Y1 e Y4 são as condições decontorno do problema), assim, tem-se que: ff   ff  

Ou seja, os termos em vermelho devem ser incorporados ao sistema de equações matricialnão singular, ou seja:

ff  ff

 

Se Y1 e Y4 fossem 0, bastaria remover as linhas e colunas dos termos dos extremos(condições de contorno). Sendo diferente de zero (e será necessário generalizar na

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 8/20

programação) deve-se incorporá-los ao termos de vetores independentes reduzido (não-singular).

Sendo o sistema acima resolvido, temos a solução do problema para elementos lineares.

2.1.2 - Matriz de rigidez para o elemento quadrático

O desenvolvimento matemático é similar ao resolvido para o elemento linear, só que agoraa função de aproximação u(x) será quadrática. Assim tem-se:u  

Resolvendo agora nosso problema com elemento linear utilizando a formulação fraca,temos:

Euação ∫ u u

 

Euação ∫ u u

 

Euação ∫ u u

 

Com:

u   u   u  

E, para o elemento quadrático, conforme visto anteriormente:

 

O que resulta num sistema de equações de três incógnitas

  e

  e

. Escrevendo o

sistema de equações na forma matricial, denominaremos a matriz dos termosdependentes de matriz de rigidez (k - 3x3) e o vetor dos termos independentes de vetorde forças (f - 3x1), assim, temos da mesma forma:

y= 

Ou seja:

  f  

f  

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 9/20

f  

Onde L representa o comprimento do elemento (. Os termos independentes f1, f2 e f3também ficam em função de coordenadas locais de cada elemento (

,

 e

), que podem

ser convertidas em coordenadas globais. Sendo assim, escrevemos na forma matricial asequações para o elemento quadrático:

f f f  

Percebe-se, então, que o método sempre resulta num sistema de equações com uma matriz

simétrica.

A estratégia de sobrepor as matrizes dos elementos, levando em conta os nós que teminfluência sobre um ou mais elementos, para transformá-las em matrizes globais é amesma apresentada anteriormente, só que agora com matrizes dos elementos 3x3. Damesma forma, também observa-se, que, mesmo sendo somente necessário resolver asduas equações do miolo, ainda será necessário implementar uma outra alteração naprogramação para incluir as condições de contorno nos termos independentes,exatamente da mesma forma como descrito no item anterior, para elementos lineares.

Resolvendo o sistema de equações global, temos a solução do problema para elementos

quadráticos (ou seja, cada elemento com 3 nós internos). Utilizando 2 elementos (cada umcom 1 grau de liberdade), teríamos de resolver 5 incógnitas (5 nós –  sendo 3 em cadaelemento, só que um deles comum aos dois elementos), sendo das 5, duas com valores jáconhecidos, que são as condições de contorno.

3ª Etapa – elaboração do código geral em Mathcad – descrição e aplicação do código

3.1 – Escolha do software

Apesar do Matlab ser mais poderoso no que se refere à operações matriciais, noMathcad também é possível realizar operações matriciais com grande número de

elementos. Assim como no Maple e no Excel, o Mathcad possui uma linguagem onde ocódigo é muito mais compreensível visualmente, além de ser possível realizar o cálculo deexpressões de forma simbólica de maneira mais simples.

Além disso, o arquivo do código é simplesmente uma planilha eletrônica combaixíssimo custo computacional. O Matlab foi abandonado por ser um programaextremamente pesado e com executáveis não tão fáceis de serem gerados, e, destamaneira, a forma de entrada de dados (input) fica prejudicada.

O Mathca peca apenas na solução eata e EDO’sque é resolvida apenas

numericamente e graficamente, mas não analiticamente. Como o projeto não envolvesolução analítica e EDO’s mas ustamente o contrário e pelo fato o autor ter ampla

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 10/20

experiência com o software, optou-se pelo desafio de criar o código MEF no Mathcad comdados de input dinâmicos e processamento rápido (arquivo em forma de planilha). Oresultado foi positivo.

3.2

– Algoritmos/Programação do código

3.2.1 – Visão geral do programa

Página 1

– dados em amarelo são as entradas do código.

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 11/20

Página 2 -Dados de saída – gráfico comparando as duas soluções, vetor u de solução (y1,y2 e y3) e o vetor xG (nós dos elementos em coordenadas globais). Percebe-se que asolução em 2 elementos lineares não foi muito satisfatória.

3.2.2 – Programação/Algoritmos

O Mathca não possui uma uantiae e funções ‘uilt-in’ ou sea á preeterminaasno software) que fossem necessárias para a solução deste código em questão. Apossibilidade de criar funções através de estruturas de repetição e condicionais sãoinfinitas, porém requerem uma maior expertise no software. A seguir as funções que oautor considera principais para programação realizada do código e solução do problema:

I – função REPLACE (função criada): repl(K,K1,r,c)

Função que insere uma matriz menor, K1, numa matriz maior, K. Os argumentos r e c dafunção representam a posição em que o primeiro elemento da matriz K1 será inserida namatriz K. A função é necessária no procedimento para sobrepor as matrizes dos elementose gerar a matriz de rigidez global.

Programação e exemplo da função repl

II

–  Função REMOVE COLUMNS e REMOVE ROWS (função criada): remcols(A,c) e

remrows(A,c).

Função ue remove linhas ou colunas e ínice ‘c’ a matriz . O arumento ‘c’ eve serum vetor. É uma função necessária para remoção das linhas e colunas relativas às‘’restrições’’ conições e contorno o prolema para erar a matriz não -singular esolução do sistema de equações reduzido (com matriz não singular).

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 12/20

 

Descrição da função e exemplo

III – Função LINEAR SOLVE (função interna do programa): lsolve(A,B)

A função interna lsolve(A,B) resolve o sistema de equações lineares do tipo Ax=B, onde Aé a matriz dos termos dependentes B vetor dos termos independentes do sistema.

IV – Função ODESOLVE (função interna do programa): odesolve(x,a)

A função interna odesolve(x,a) resolve numericamente uma EDO dadas condições decontornos necessárias. Para a

solução é preciso criar um ‘solve loc’ no Mathca amaneira demonstrada abaixo.

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 13/20

3.3 Aplicação/Exemplos e validações

3.3.1 – 3 elementos lineares com EDO:

 

Ou seja, A = 1, B = -1 e C=1. Condições de contorno: Y1 = 1 e Y4 = 0 e domínio: .

Solução no código Matlab do servidor:

Solução código Mathcad elaborado:

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 14/20

 

Diferença 0%.

3.3.2

–  4 elementos lineares com EDO:

 

Ou seja, A = 1, B = -1 e C=1. Condições de contorno: Y1 = 2 e Y5 = 2 e domínio: .

Solução no código Matlab do servidor:

Obs: pequena inconsistência encontrada no código Galerkin_FormaFraca_4elem_linear.mdo servidor, nas linhas 85 a 89, onde há u3 deve ser substituído por u4.

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 15/20

 

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 16/20

 

Solução código elaborado no Mathcad:

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 17/20

3.3.3 –  2 elementos quadráticos com EDO:  

Ou seja, A = 1, B = -1 e C=1. Condições de contorno: Y1 = 3 e Y5 = 2 e domínio:

.

Solução no código Matlab do servidor:

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 18/20

Solução código elaborado MathCad:

Observa-se que a aproximação com 2 elementos quadráticos já é mais eficiente do quecom 4 elementos lineares.

3.3.4 – 2 elementos quadráticos / 6 elementos lineares / 8 elementos lineares, com EDO:  

Ou seja, A = 1, B = -1 e C=1. Condições de contorno: Y1 = 2 e Y5 = 2 e domínio: .

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 19/20

 

2 elem quadráticos

6 elementos lineares (pior aproximação em relação a 2 quadráticos)

7/17/2019 Elementos finitos: estudo de caso

http://slidepdf.com/reader/full/elementos-finitos-estudo-de-caso 20/20

 

Solução com 8 elementos lineares.

2. Equação diferencial

– proposta avançada

A ideia agora consiste em resolver uma equação diferencial de segunda ordem mais geral,da seguinte forma:

f Tornando a solução do problema mais generalizada. O arquivo da proposta avançada éEQ1DA.xmcd, e, da mesma forma do problema anterior, também foi possível obterresultados positivos tanto para elementos quadráticos quanto para lineares.

Observa-se, ao final, que o elemento quadrático é muito mais poderoso para a solução doMEF. Em alguns casos, verificou-se que 2 elementos quadráticos geravam aproximaçãomelhor do que 8 elementos lineares.