Apostila metodos numericos_formatada_vf

INTRODUÇÃO AOS MÉTODOS

NUMÉRICOS

Celme Torres F. da Costa UNIVERSIDADE FEDERAL DO CEARÁ

Campus Cariri 2010

Introdução aos Métodos Numéricos

2

Sumário

1. Introdução aos Métodos Numéricos ................................................................... 4

1.1. Soluções Analíticas .................................................................................... 5 1.2. Método das diferenças Finitas ................................................................... 6

1.3. Método dos elementos finitos .................................................................... 6 1.4. Método dos momentos (MOM) .................................................................. 7

1.5. Método dos elementos analíticos ............................................................... 7 1.6. Método das camadas finitas ...................................................................... 7

1.7. Método de elementos de contorno (MEC) ................................................. 8 1.8. Método dos volumes finitos........................................................................ 8

2. Condições de Contorno .................................................................................... 10 2.1. Condições de contorno naturais .............................................................. 11

2.2. Condições de contorno geométricas ........................................................ 11

3. Sistemas de Equações Diferenciais ................................................................. 12

3.1. Sistemas de Equações Diferenciais ......................................................... 12 3.1.1. Sistemas homogêneos ...................................................................... 12

3.1.2. Sistemas não homogêneos ............................................................... 14 3.2. Conceito de autovalor e autovetor ........................................................... 15

3.3. Diagonalização de matrizes ..................................................................... 18 3.4. Solução geral de sistemas homogêneos ................................................. 19

3.5. Solução geral de sistemas não homogêneos .......................................... 20 3.5.1. Exercícios e Aplicações de sistemas de equações diferenciais ........ 21

3.6. Solução de sistemas lineares .................................................................. 22 3.6.1. Exercícios .......................................................................................... 27

4. Métodos dos Resíduos Ponderados ................................................................ 28 4.1. Método de colocação pontual .................................................................. 29

4.2. Método dos subdomínios ......................................................................... 32 4.3. Método dos mínimos quadrados .............................................................. 34

4.4. Método de Galerkin .................................................................................. 36 4.5. Exercícios ................................................................................................ 38

5. Introdução aos Métodos Variacionais............................................................... 39 5.1. Cálculo variacional ................................................................................... 39

5.2. Máximo e mínimo de funções .................................................................. 39 5.3. Funcional ................................................................................................. 43

5.3.1. Funcional de uma função .................................................................. 44 5.3.2. Determinação do Funcional - Equação de Euler ............................... 47

5.3.3. Exercícios .......................................................................................... 51

6. Métodos Variacionais ....................................................................................... 53

6.1. Método de Ritz ......................................................................................... 53 6.1.1. Exercícios .......................................................................................... 58

6.2. Método de Rayleigh-Ritz .......................................................................... 59

7. Métodos das Diferenças Finitas ....................................................................... 66

7.1. Derivadas de ordem superior ................................................................... 68 7.2. Exercícios ................................................................................................ 73

8. Métodos dos Elementos Finitos ....................................................................... 75

Profª Celme Torres

3

8.1. Sistemas discretos e funções de base ..................................................... 78

8.1.1. Funções de base ................................................................................ 78 8.1.2. Exemplos de sistemas discretos ........................................................ 81

8.1.2.1. Mola elástica ................................................................................ 81 8.1.2.2. Torção em barra circular .............................................................. 82

8.1.2.3. Fluxo de um fluido em um tubo .................................................... 83

9. Métodos dos Elementos Finitos Bidimensional ................................................. 94

9.1. Modelos de elementos finitos ................................................................... 95 9.1.1. Funções de forma .............................................................................. 95

9.1.2. Elemento triangular linear................................................................... 96 9.1.3. Teorema de Green-Gaus ................................................................... 98

9.1.4. Imposição das condições de contorno ............................................. 100 9.1.5. Determinação da matriz de rigidez ................................................... 101

9.1.6. Determinação do vetor de carga nodal............................................. 102


4

1. Introdução aos Métodos Numéricos

Os modelos são ferramentas projetadas para representar uma versão simplificada de um

problema real. É uma tentativa para compreensão dos processos físicos, químicos e bioló-

gicos traduzidos em termos matemáticos. O objetivo da modelagem é o de prever ou predi-

zer cenários onde estão envolvidos variáveis desconhecidas, como, por exemplo, a varia-

ção da carga hidráulica ou distribuição de concentrações de espécies químicas em um sis-

tema aqüífero no tempo e no espaço (BEDIENT et al., 1994), variação de temperaturas,

análise de estruturas, entre outras.

No desenvolvimento de um modelo, o primeiro passo é o desenvolvimento de um modelo

teórico ou conceitual, que consiste na descrição dos processos envolvidos e no conheci-

mento do comportamento do sistema analisado. O próximo passo é traduzir esse modelo

teórico em termos matemáticos, o chamado modelo matemático, que é o conjunto de equa-

ções diferenciais associada ao conjunto de equações de condições iniciais e de condições

de contorno. Finalmente, a solução das equações sujeitas às condições de contorno impos-

tas é obtida utilizando métodos analíticos ou numéricos.

Quando os modelos consideram as variáveis de entrada e/ou os parâmetros dependentes de

funções aleatórias (probabilidade) esses modelos são definidos como modelos ESTO-

CÁSTICOS ou PROBABILÍSTICOS, do contrário, são DETERMINÍSTICOS. Os mo-

delos baseados em resultados experimentais, são ditos modelos EMPÍRICOS (COSTA,

2000).

Uma vez que o modelo teórico é formulado e os parâmetros apropriados tenham sido de-

terminados, o próximo passo é achar uma solução para as equações que governam o pro-

blema, considerando as condições iniciais e de contorno. As técnicas mais utilizadas na

solução destas equações são: solução analítica, método dos momentos, método das cama-

das finitas, elementos de contorno, diferenças finitas e elementos finitos, descritas a seguir.

Profª Celme Torres

5

1.1.1. Soluções Analíticas

Qualquer função f definida em algum intervalo I, que, quando substituída na equação dife-

rencial reduz a equação a uma identidade, é chamada de solução analítica da equação di-

ferencial no intervalo. Porém, nem toda equação diferencial possui necessariamente uma

solução analítica.

Quando resolvemos uma equação diferencial, normalmente obtemos uma família de curvas

ou funções contendo um parâmetro arbitrário tal que, cada membro da família é uma solu-

ção da equação diferencial.

As soluções das equações diferenciais são divididas em Explícitas ou Implícitas. Uma so-

lução explícita é qualquer função y(x) que verifique a equação num intervalo a x b .

Ou seja, uma solução na qual a variável dependente é expressa somente em termos da vari-

ável independente e das constantes. Para nossos propósitos, vamos pensar que uma solução

explicita seja da forma y f (x) , a qual pode ser manipulada, calculada e diferenciável por

meio das regras padrão. Isto é, a variável dependente pode ser isolada e igualada a uma ex-

pressão, a qual é função apenas da variável independente, não ambígua.

Uma solução implícita é uma relação g(x, y) 0 que verifique a equação. Dizemos que

uma relação g(x, y) 0 é uma solução implícita de uma EDO em um intervalo I, se ela

define uma ou mais soluções explícitas em I. Isto é, a variável dependente (função) não

pode ser isolada e igualada a uma expressão que dependa apenas da variável independente,

ou quando isto for possível então a expressão será ambígua.

As soluções analíticas são ideais para um rápido cálculo preliminar, são válidas para estu-

dos de sensibilidade e úteis para checar resultados de análises mais complexas. Toda base

teórica para determinação das soluções analíticas de equações diferenciais foram descritas

e assimiladas em disciplinas anteriores envolvendo equações diferencias ou matemática

aplicada.


6

1.1.2. Método das diferenças Finitas

O método das diferenças finitas substitui a equação diferencial por um conjunto de equa-

ções lineares algébricas que são resolvidas de forma aproximada (McDONALD; HAR-

BAUGH, 1988). As idéias fundamentais foram estabelecidas e usadas por matemáticos do

século XVIII, tais como Taylor e Lagrange (GOMES, 2000). No método das Diferenças

Finitas, uma equação diferencial (de natureza contínua) é substituída por uma série de e-

quações algébricas, chamadas de diferenças finitas, em pontos discretos. Esse método é o

mais antigo e o mais divulgado e foi o primeiro a ser usado para a solução sistemática de

problemas de água subterrânea.

Figura 1. Malha de diferenças finitas.

1.1.3. Método dos elementos finitos

O método dos elementos finitos substitui a equação diferencial por uma formulação vari-

acional do problema, que leva a resolução aproximada de um sistema de equações lineares

algébricas (ZIENKIEWICZ, 1977). Baseia-se na divisão do domínio do problema em sub-

domínios, o que permite a representação de domínios complexos como um conjunto de

subdomínios mais simples – chamados Elementos Finitos. Cada elemento finito é conecta-

do aos elementos vizinhos através de nós. O conjunto dos elementos e dos nós é chamado

de malha. A Figura 1.0 mostra uma malha de elementos finitos.

i, ji-1, j i+1, j

i, j+1

i, j-1i-1, j-1 i+1, j-1

i-1, j+1 i+1, j+1 i+2, j+2

i+2, j-2

x

y

x

y

Profª Celme Torres

7

Figura 2 – Malha de elementos finitos.

1.1.4. Método dos momentos (MOM)

O método dos momentos (MOM) tem a sua origem – como muitas outras idéias em esta-

tística - no trabalho precursor de Lord Karl PEASON (1857). Pearson havia construído o

seu sistema de famílias de distribuições de probabilidade e necessitava um procedimento

automático para, uma vez escolhida a família apropriada, estimar os parâmetros da particu-

lar curva (ou função de densidade) que melhor se ajustaria aos dados. Como tais parâme-

tros eram, geralmente, relacionados aos momentos da distribuição através de funções algé-

bricas não muito complicadas, ele propôs que se igualasse os valores dos momentos amos-

trais aos momentos teóricos, conforme descritos por tais funções. Desse modo, indo-se até

o momento de ordem igual ao número de parâmetros a estimar, se obtinha um sistema de

tantas equações quanto incógnitas – os parâmetros – a resolver.

1.1.5. Método dos elementos analíticos

O método dos elementos analíticos é baseado no princípio da superposição de várias fun-

ções analíticas. Este princípio consiste da adição de soluções individuais para cada elemen-

to. A solução do problema é obtida através da adição de todas as influências dos elementos

analíticos individuais, representados pelas funções analíticas, que correspondem às caracte-

rísticas do sistema em estudo.

1.1.6. Método das camadas finitas

O método das camadas finitas, descrito por ROWE e BOOKER (1985 a, b; 1987), se a-

plica a situações onde a hidroestratigrafia pode ser idealizada como composta de camadas


8

horizontais, com as propriedades do solo sendo as mesmas em qualquer local na camada. O

método das camadas finitas é muito utilizado em problemas de dimensionamento de pavi-

mentos.

1.1.7. Método de elementos de contorno (MEC)

O Método dos elementos de contorno ou Boundary Element Method – BEM, é um méto-

do computacional para a solução de sistemas de equações diferenciais, formuladas em

forma integral. É aplicado em diversas áreas da engenharia, como em mecânica dos flui-

dos, acústica, eletromagnetismo e estudo de fraturas.

O método possui melhor desempenho que o método dos elementos finitos em certas cir-

cunstâncias, como por exemplo, quando o domínio de estudo for infinito ou semi-infinito.

No método de elementos de contorno (MEC), primeiramente ocorre a transformação da

equação diferencial parcial que descreve o comportamento da incógnita no interior e no

contorno do domínio, em uma equação integral que envolve somente incógnitas em seu

contorno. A seguir, faz-se a discretização do contorno em elementos de superfície e, por

fim, encontra-se a solução do sistema de equações algébricas resultantes. As principais

vantagens do MEC estão no fato de que somente o contorno (ou contornos) do domínio

deve ser dividido em sub-regiões (discretizado). Nos outros métodos (Elementos Finitos e

Diferenças Finitas), todo o domínio deve ser discretizado. Desta forma, a dimensão do

problema é efetivamente reduzida em uma dimensão (BREBBIA e SKERGET, 1984).

1.1.8. Método dos volumes finitos

O método dos volumes finitos é um método de resolução de equações às derivadas parci-

ais baseado na resolução de balanços de massa, energia e quantidade de movimento a um

determinado volume de meio contínuo.

Este método evoluiu das diferenças finitas e não apresenta problemas de instabilidade ou

convergência, por garantir que, em cada volume discretizado, a propriedade em questão

(por exemplo, a massa) obedece à lei da conservação de massa. É um método largamente

Profª Celme Torres

9

utilizado na resolução de problemas envolvendo transferência de calor ou massa e em me-

cânica dos fluidos.

Figura 1.1 - Diagrama de métodos de solução de um problema físico.


10

2. Condições de Contorno

O objetivo da maioria das análises numéricas é determinar funções desconhecidas, chama-

das variáveis dependentes, que satisfaçam um dado conjunto de equações diferenciais em

um domínio ou região e algumas condições de contorno na borda do domínio. Um domínio

é uma coleção de pontos no espaço com propriedades que se P é um ponto no domínio, en-

tão todos os pontos próximos a P pertencem ao mesmo domínio. Esta definição implica

que um domínio consiste somente de pontos internos. Os símbolos Ω são utilizados para

denotar um domínio arbitrário e Γ para determinar seu contorno.

Γ (contorno)

Ω (domínio)

Figura 1.2 – Esquema de apresentação do domínio e seu contorno.

Em matemática, no ramo de equações diferenciais, um problema de valor de contorno ou

de fronteira é uma equação diferencial provida de um conjunto de restrições adicionais,

as chamadas condições de contorno. Uma solução para um problema de valor de contorno

é aquela que seja solução da equação diferencial e satisfaça as condições de contorno.

Problemas de valor de contorno surgem em diversos ramos, assim como toda equação dife-

rencial também o terá. Problemas envolvendo a equação de onda, bem como a determina-

ção dos modos normais, são frequentementes classificados como problemas de valor de

contorno. Um problema de valor de contorno deve ser bem determinado. Isto é, dado certas

condições para o problema, haverá então solução única, que depende continuamente das

condições citadas.

A solução de qualquer equação diferencial dependente do tempo (transiente) requer a espe-

cificação de condições na fronteira do sistema estudado (condições de contorno) e as con-

Profª Celme Torres

11

dições no início do processo físico (condição inicial). A definição precisa das condições de

contorno e inicial é a parte mais importante para os processos de modelagem.

Quando as variáveis dependentes são funções de uma variável independente (x), o domínio

é um segmento de reta, dito, nesse caso, unidimensional, e os pontos limites são chamados

pontos de contorno. Quando as variáveis dependentes são funções de duas variáveis inde-

pendentes (x, y), o domínio é bidimensional, sendo dado por uma superfície ou área, ge-

ralmente um plano.

2.1.1. Condições de contorno naturais

As condições de contorno naturais se referem à derivada da função que está sendo analisa-

da.

No caso de escoamento corresponde aos valores de fluxo no contorno do domínio, assim

como os pontos de perda ou ganho (poços, vazões, carga hidráulica). No caso de análise de

estruturas as condições de contorno naturais representam as forças concentradas e distribu-

ídas existentes no sistema físico analisado.

2.1.2. Condições de contorno geométricas

As condições de contorno geométricas são os valores da função ao longo do contorno ou

domínio. No caso de análise de estruturas correspondem aos apoios e as restrições do sis-

tema.


12

3. Sistemas de Equações Diferenciais

3.1.1. Sistemas de Equações Diferenciais

As equações diferenciais ordinárias simultâneas envolvem duas ou mais equações que con-

tém derivadas de duas ou mais funções incógnitas de uma única variável independente. Se

x, y e z são funções da variável t, então:

2

2

2

2

d x4 5x y

dt

d y2 3x y

dt

e, 2

x ' 3x y ' z ' 5

x ' y ' 2z ' t

x y ' 6z ' t 1

São exemplos de sistemas de equações diferenciais.

Soluções de um Sistema

Uma solução de um sistema de equações diferenciais é um conjunto de funções diferenciá-

veis x f (t) , y g(t) , z h(t) , etc., que satisfaz cada equação do sistema em algum in-

tervalo I.

Os sistemas de equações diferenciais podem ser resolvidos por:

(i) Método dos operadores, através da eliminação algébrica sistemática;

(ii) Método da transformada de Laplace

(iii) Método de autovetores e autovalores de matrizes

Sistemas homogêneos

Um sistema de equações diferenciais é um conjunto de n ED, com uma variável indepen-

dente e n variáveis dependentes, que podem ser escritos da seguinte forma:

' ' '11 1 2 n 1 2 n

' ' '22 1 2 n 1 2 n

' ' 'nn 1 2 n 1 2 n

dyF y , y ,..., y , y , y ,..., y , x

dx

dyF y , y ,..., y , y , y ,..., y , x

dx

dyF y , y ,..., y , y , y ,..., y , x

dx

Profª Celme Torres

13

onde F1, F2, ..., Fn são quaisquer funções de (2n+1) variáveis reais, que definem o sistema.

O sistema de EDO lineares pode ser escrito numa forma simplificada, usando a notação

vetorial,

onde y é um vetor com n componentes, cada uma delas função de x, e F é um vetor com n

componentes, funções de y, x e y’.

Um caso especial na teoria de ED é quando o vetor de funções F, no sistema de equações

tem a forma:

onde A é uma matriz quadrada n x n de funções que dependem unicamente de x, f é um

vetor

com n componentes dependentes de x.

Exemplo 3.1

1. Escreva em forma matricial os sistemas

(a) 1 1 2

2 1 2

y ' x x

y ' 4x x

Na forma vetorial, o sistema pode ser escrito como:

O vetor f é nulo e a matriz A é igual a 1 1

A4 1

dy dyF y, x,

dx dt

F Ax f

1 1

2 2

y ' x1 1

y ' 4 1 x

Princípio da Superposição

Seja 1 2 nx , x ,..., x um conjunto de vetores solução de um sistema de EDO homogêneo

em um intervalo I. Então a combinação linear

1 1 2 2 n nx c x c x ... c x

Onde, c1, c2 e cn são constantes arbitrárias, é também solução do sistema no intervalo.


14

(b)

t

t

dx3x 4y e sin 2t

dt

dy5x 9y 4e cos 2t

dt

A forma matricial temos:

t

t

x ' 3 4 x e sin 2t

y ' 5 9 y 4e cos 2t

t

t

3 4 e sin 2tdXX

5 9dt 4e cos 2t

Onde, x

Xy

2. Verifique que

2t

1 2t

ex

e

e

6t

2 6t

3ex

5e

são soluções do sistema de EDO’s

1 3X ' X

5 3

, no intervalo ( , ) , onde x

Xy

.

Temos que,

2t

1 2t

2ex '

2e

2t 2t 2t 2t

1 12t 2t 2t 2t

1 3 2e e 3e 2eAx x '

5 3 2e 5e 3e 2e

6t

2 6t

18ex '

30e

6t 6t 6t 6t

2 26t 6t 6t 6t

1 3 3e 3e 15e 18eAx x '

5 3 5e 15e 15e 30e

Logo, x1 e x2 são soluções da EDO.

Sistemas não homogêneos

Para sistemas não homogêneos, uma solução particular xp em um intervalo I é qualquer

vetor, sem parâmetros arbitrários, cujos elementos são funções que satisfazem o sistema,

dXA(t)X F(t)

dt

Profª Celme Torres

15

Exemplo 3.2

Verifique que o vetor p

3t 4X

5t 6

é uma solução particular do sistema não-homogêneo

1 3 12t 11X ' X

5 3 3

, no internalo ( , ) .

Solução

Temos que, p

3X '

5

p

1 3 12t 11 1 3 3t 4 12t 11X

5 3 3 5 3 5t 6 3

12t 14 12t 11 3

2 3 5

Logo, o vetor p

3t 4X

5t 6

é solução particular do sistema homogêneo.

3.1.2. Conceito de autovalor e autovetor

Dado o sistema:

8 2 1 14

2 11 3 35

Teorema

Seja 1 2 nx , x ,..., x um conjunto de vetores solução de um sistema de EDO homogêneo

em um intervalo I e seja Xp um vetor arbitrário solução do sistema não homogêneo no

mesmo intervalo. Então a combinação linear

1 1 2 2 n n px c x c x ... c x X

É também uma solução do sistema não-homogêneo no intervalo, para quaisquer valores

das constantes c1, c2 e cn .

vetor de entrada

vetor de saída


16

Tendo,

O vetor de saída está na mesma direção e sentido que o vetor de entrada. Neste caso:

Sendo um escalar e v um vetor.

Logo, v é chamado de AUTOVETOR de A e um AUTOVALOR de A associado ao ve-

tor v.

Considere a matriz quadrada . Se for um número complexo e v um vetor com-

plexo não nulo de dimensões n, satisfazendo a identidade:

8 2 1 12 112

2 11 2 24 2

Av v

n nA R

Av v

vetor de saída vetor de entrada

vetor de saída

vetor de entrada

Teorema de Cayley-Hamilton

Qualquer matriz quadrada satisfaz sua própria equação característica. Isto é, se e somen-

te se,

n n 1 2

n n 1 2 1 odet(A I) b b ... b b b

Então, n n 1 2

n n 1 2 1 ob A b A ... b A b A b I 0

Profª Celme Torres

17

Então é um autovalor de A e v é um autovetor de A associado ao autovalor de .

Se A é uma matriz n x n sobre R e I é a matriz identidade de mesma ordem de A, defini-

mos o polinômio característico de A como:

f ( ) det( I A)

Os autovalores de A podem ser obtidos a partir da existência de escalares e vetores não-

nulos v para os quais Av v .

Este sistema pode ser escrito como: Av Iv

(A I)v 0

Sendo I é a matriz identidade. Para que a solução não seja nula, o determinante deve ser

igual a zero (equação característica), logo:

det( I A) 0

OBS: O número de autovalores é sempre igual a ordem da matriz.

Exemplo 3.3

(1) Calcule o autovalor da EDO: 2y'' 6y ' 4y 0

Pelo teorema de Cayley-Hamilton

2

1

2

2 6 4 0

logo,

2

1

(2) Calcule os autovalores e autovetores da matriz 3 4

A1 7

Pelo teorema de Cayley-Hamilton

2 2

1 2

3 4det(A I) ( 5) 0

1 7

5

Para encontrar os autovetores correspondentes ao autovalor, fazemos:


18

1

2

1 2

1 2

(A I)k 0

(A 5I)k 0

k3 4 1 05 0

1 7 0 1 k

2k 4k 0

k 2k 0

Por esse sistema temos que 1 2k 2k . Escolhendo a solução trivial 2k 1 , temos o úni-

co autovetor do sistema, dada por:

1

2K

1

3.1.3. Diagonalização de matrizes

Seja C uma matriz qualquer

a b

Cc d

Se multiplicarmos C por 1

0

a b 1 a

Cc d 0 c

1ª coluna de C

Multiplicando C por 0

1

a b 0 b

Cc d 1 d

2ª coluna de C

Suponha uma matriz A2 x 2 cujos autovetores e autovalores são:

1 1

2 2

v

v

Se uma matriz E 2 x 2 é formada pelos autovetores de A e se,

1

av

c

e 2

bv

d

Profª Celme Torres

19

a b

Ec d

Considere a matriz resultante de 1E AE , multiplicada por 1

0

1

1

av

c

1E A E

0

.

Temos, 1

1E Av

Onde, 1 1 1Av v

1 1

1 1 1 1E v E v

Se, 1

1E v

0

1

1

1E v

0

Assim, 11

1 1E v0

1ª coluna de 1

1E Av

Do mesmo modo se multiplicarmos 1E AE

por 0

1

, obtemos a 2ª coluna de 1

1E Av

Portanto, 11

2

0E AE

0

, diagonaliza a matriz A.

3.1.4. Solução geral de sistemas homogêneos

Seja 1 2 n, ,..., n autovalores reais distintos de uma matriz de coeficientes A do sistema

X' AX , onde x

Xy

, e sejam 1 2 nk ,k ,...,k os autovetores correspondentes. Então a

solução geral do sistema no intervalo ( , ) é dada por:

1 2 nt t t

1 1 2 2 n nX c k e c k e ... c k e

Exemplo 3.4

Resolva

dx2x 3y

dt

dy2x y

dt

A equação característica é


20

22 3det(A I) 3 4 0

2 1

Ao autovalores do sistema são: 1 1 e 2 4

Para 1 1 , temos para (A I)K 0

1 2

1 2

3k 3k 0

2k 2k 0

Resolvendo o sistema: 1 2k k

Quando 2k 1 , o autovetor correspondente é 1

1k

1

Para 2 4 , temos para (A I)K 0

1 2

1 2

2k 3k 0

2k 3k 0

Resolvendo o sistema: 21

3kk

2

Quando 2k 2 , o autovetor correspondente é 2

3k

2

Como a matriz A de coeficientes é uma matriz 2 x 2, achamos duas soluções linearmente

independentes na forma:

t

1

1X e

1

e 4t

2

3X e

2

A solução geral do sistema de EDO homogêneo é dada por:

1 2t t t 4t

1 1 2 2 1 2

1 3X c k e c k e c e c e

1 2

t 4t

1 2

t 4t

1 2

x(t) c e 3c e

y(t) c e 2c e

3.1.5. Solução geral de sistemas não homogêneos

Seja Xp uma solução do sistema não-homogêneo X' A(t)X F(t) em um intervalo I, e

denotemos por

p 1 1 2 2 n nX c X c X ... c X

A solução geral no mesmo intervalo, do sistema homogêneo X' A(t)X correspondente.

Define-se solução geral do sistema não-homogêneo no intervalo como c pX X X .

Profª Celme Torres

21

Exemplo 3.5

Verificamos que o vetor p

3t 4X

5t 6

é uma solução particular do sistema não- homo-

gêneo 1 3 12t 11

X ' X5 3 3

, no internalo ( , ) .

Para a parte homogênea da EDO temos

dx1x 3y

dt

dy5x 3y

dt

Resolvendo aplicando autovetor e autovalores de matriz, temos:

2t

1

1X e

1

e 6t

2

3X e

5

A solução não homogênea Xn é dada por 2t 6t

n

1 3X e e

1 5

A solução geral do sistema de EDO’s pela definição é n pX X X , logo:

2t 6t1 3 3t 4

X e e1 5 5t 6

2t 6t

2t 6t

x(t) e 3e 3t 4

y(t) e 5e 5t 6

Exercícios e Aplicações de sistemas de equações diferenciais

(1) Encontre a solução geral dos sistemas de equações diferenciais

a)

dxx 2y

dt

dy4x 3y

dt

b) 0,1 0,075 1,5

Q' Q0,1 0,2 3

c)

dx2x y

dt

dy5x 4y

dt


22

d) 1 5 sin t

X ' X1 1 2cos t

(2) Sejam T1 = T1(t) e T2 = T2(t) as temperaturas, no instante t, nos ambientes 1 e 2, res-

pectivamente. Admita que as temperaturas estejam variando nas seguintes taxas

dT1T1 4T2

dt

dT22T1 3T2

dt

Suponha que T1(0) = 20 e T2(0) = 10.

Determine as temperaturas, no instante t, nos ambientes 1 e 2.

(3) O tanque A contém 50 litros de água que foram dissolvidos 25 gramas de sal. Um se-

gundo tanque B, contém 50 litros de água pura. Bombeia-se o líquido para dentro e para

fora dos tanques as taxas indicadas na figura. O sistema de equações diferenciais que des-

creve o problema de tanques de mistura é

11 2

21 2

dx 2 1x x

dt 25 50

dx 2 2x x

dt 25 25

onde as condições iniciais são dadas por 1x (0) 25 e 2x (0) 0 . Aplicando autovetor e

autovalor de matriz encontre a solução do sistema de equações diferenciais

3.1.6. Solução de sistemas lineares

Vamos considerar a solução do sistema linear: usando o fato de que qualquer

matriz quadrada pode ser expressa como um produto de uma matriz triangular infe-

rior [L] e uma matriz triangular superior [U], temos:

(3.1)

Considerando n = 3,

(3.2)

[A]x B

n n[A]

[A] [L][U]

11 12 13 11 12 13

21 22 23 21 22 23

31 32 33 31 32 33

a a a l 0 0 1 u u

a a a l l 0 0 1 u

a a a l l l 0 0 1

Profª Celme Torres

23

Resolvendo o sistema para as variáveis e teremos:

(3.3)

(3.4)

(3.5)

(3.6)

Uma vez que decomposição da matriz A é feita, a resolução do sistema é feita da seguinte

forma:

(3.7)

(3.8)

Definindo um vetor [Z], como:

(3.9)

Portanto,

(3.10)

O sistema pode ser resolvido por substituição:

onde, i = 1,..., n (3.11)

Encontrados os valores de Z podemos calcular X usando:

(3.12)

Ou,

(3.13)

A solução é dada então por:

ijl iju

11 11l a 21 21l a 31 31l a

1212

11

au

l 22 22 21 12l a l u

32 32 31 12l a l u

j 1

ij ij ik kj

k 1

i 1

ij ij ik kj

k 1ii

l a l u se i j

1u a l u se i j

l

[A]X B

[L][U]X B

[U]X Z

[L]Z B

i 1

i ik k

k 1i

ii

b l z

zl

[U]X Z

1 12 2 13 3 1n n 1

2 23 3 2n n 2

3 3n n 3

n n

x u x u x ... u x z

x u x ... u x z

x ... u x z

x z


24

(3.14)

A solução de sistemas lineares pode ser resolvida também pelo Método de Gauss.

Exemplo 3.6

Seja o sistema:

Para eliminar x1 da segunda linha, basta multiplicar a primeira linha por 5/10 , isto é, por

a21/a11 , o que torna o coeficiente de x1 na primeira linha igual ao de x1 na segunda linha.

Em seguida basta subtrair a primeira linha modificada da segunda, o que zerará o coefi-

ciente de x1, eliminando-o da segunda equação.

Repetir essa operação para a terceira linha, multiplicando a primeira por a31/a11 , antes

de subtraí-la da terceira. O mesmo para a quarta linha, quinta linha etc...

Na enésima linha multiplica-se a primeira por an1/a11 e subtrai-se da enésima, eliminando-

se, dessa linha, a variável x1.

O mesmo raciocínio é repetido para a variável x2 da terceira equação em diante, para a

variável x3 da quarta em diante, até que a penúltima variável, xn-1 seja eliminada da última

equação.

Nesse momento o sistema estará triangularizado.

A partir desse ponto começa a segunda fase, “backward”, quando se calcula, na última

equação, a última variável, xn ; leva-se xn à penúltima equação e se calcula a penúltima

variável, xn-1 ; leva-se xn e xn-1 à antepenúltima equação e se calcula a antepenúltima vari-

ável, xn-2 , etc... até ser calculada a primeira variável na primeira equação.

Sendo Li a linha i , a eliminação da variável xj dessa linha se dará pela operação:

Nova Li = Li – aij/ajj . Lj , variando-se j de 1 a n-1 e i de j+1 até n.

No sistema apresentado:

n 1

i i i,k 1 k 1

k 1

x z u x

1 2 3

1 2 3

1 2 3

10x 5x 3x 17

5x 8x 2x 19

2x 3x 8x 0

1 2 3

1 2 3

1 2 3

10x 5x 3x 17

5x 8x 2x 19

2x 3x 8x 0

Profª Celme Torres

25

Observe que após eliminarmos x1 da segunda à última equação, esquecendo-se de L1, fica-

se com um sistema de 2 equações a 2 incógnitas, isto é, diminuiu-se de 1 a ordem do sis-

tema anterior; no caso passamos a ter um sistema de duas equações a duas incógnitas.

Para se eliminar x2 da terceira equação, trabalha-se da mesma maneira:

Tem-se agora um sistema de 1 equação a 1 incógnita.

O sistema completo fica sendo:

O sistema está triangularizado e começa, nesse momento, a segunda parte, a solução

propriamente dita do sistema.

2 2 1L L 5/10L

2 1 2 3L 0x (8 2,5)x (2 1,5)x 19 8,5

2 1 2 3L 0x 5,5x 0,5x 10,5

3 1 2 3L 0x 2,0x 7,4x 3,4

2 3

2 3

5,5x 0,5x 10,5

2,0x 7,4x 3,4

3 3 2L L 2,0 / 5,5L

3 2 3L 0x (7,4 0,182)x 3,4 3,818

37,218x 7,218

1 2 3

1 2 3

1 2 3

10x 5x 3x 17

0x 5,5x 0,5x 10,5

0x 0x 7,218x 7,218


26

Solução:

Algumas Considerações:

Numa dada matriz, quando se subtrai de uma linha outra linha multiplicada por uma

constante, o determinante da matriz não se altera. Como no processo de eliminação,

o que se fez foi sempre subtrair de uma linha outra linha multiplicada por uma cons-

tante, conclui-se que ao ser triangularizada, a matriz original manteve constante o

valor do determinante.

O determinante de matriz triangular é dado pelo produto dos elementos da diagonal

principal. Dessa maneira, o determinante da matriz original, antes de ser triangulari-

zada pelo método de Eliminação de Gauss, será dado pelo produto dos elementos da

diagonal principal da matriz triangular resultante.

No processo de eliminação da variável xj da linha i, faz-se a transformação

ij

i j

jj

aL L

a . No caso de ajj ser zero, deve-se renumerar as linhas para se conseguir ajj

diferente de zero.

Na verdade, o ideal seria conseguir uma linha j tal que, além de ter ajj diferente de

zero, tivesse esse elemento com o maior módulo possível e os demais elementos da

linha j o menor possível para que alterasse o mínimo a linha i, da expressão

ij

i j

jj

aL L

a . Dessa maneira diminuem os erros propagados de uma linha para outra no

processo de eliminação. Levar em conta a escolha adequada da linha que será a linha

j, para minimizar a propagação dos erros aumentando a precisão do método, é a

chamada Condensação Pivotal.

3

2

1

x 1,0

x (10,5 0,5( 1,0)) / 5,5 2,0

x (17 5.2,0 3.( 1,0)) /10 1,0

1x 1,0 2x 2,0 3x 1,0

Profª Celme Torres

27

Exercícios

Resolva os sistemas abaixo aplicando o método de Gauss

a) b)

1 2 3

1 2 3

1 2 3

2x 6x x 7

x 2x x 1

5x 7x 4x 9

x 3y 2z 7

4x y 3z 5

2x 5y 7z 19


28

4. Métodos dos Resíduos Ponderados

São métodos de aproximação utilizados para resolver equações diferenciais. Estes métodos

têm diversos procedimentos, dentre os quais se destacam o Método de Galerkin, o Método

da Colocação, o Método do Subdomínio e o Método dos Mínimos Quadrados.

A solução aproximada, para o caso de problemas unidimensionais, é do tipo:

(4.1)

Onde são as funções que satisfazem as condições de contorno do problema e são

coeficientes que serão determinados a partir da imposição de que o erro resultante do uso

desta solução aproximada seja mínimo. A seguir essa solução é substituída diretamente na

equação diferencial.

Como a solução aproximada não satisfaz a equação, haverá como resultado da

substituição um erro ou resíduo .

Esse erro é então multiplicado por uma função ponderadora e a integral desse pro-

duto é igualada a zero:

(4.2)

Onde D é o domínio do problema e o número de funções ponderadoras é igual ao

número de coeficientes desconhecidos .

Os métodos dos resíduos ponderados mais conhecidos são o método de colocação pontual,

método dos subdomínios, método dos mínimos quadrados e método de Galerkin. Veremos

a seguir cada um deles.

n

n i i

i 1

h (x) a N (x)

iN (x) ia

nh (x)

R(x)

iW (x)

i

D

W (x)R(x)dx 0

iW (x)

ia

Profª Celme Torres

29

4.1.1. Método de colocação pontual

Para este método, as funções ponderadoras são funções ponderadoras do tipo de

Delta-Dirac .

Uma das propriedades da função Delta-Dirac é:

(4.3)

Sendo assim, se a função aproximada tem n coeficientes , então:

(4.4)

Portanto,

(4.5)

Ou ainda,

(4.6)

Expandindo a expressão acima teremos:

(4.7)

iW (x)

ix x

i if (x) x x dx f (x )

ia

1 2 nR R x,a ,a ,...,a

i 1 2 n

D

W (x)R x,a ,a ,...,a dx 0

i 1 2 n

D

x x R x,a ,a ,...,a dx 0

1 1 2 n

2 1 2 n

3 1 2 n

n 1 2 n

R x ,a ,a ,..., a 0

R x ,a ,a ,..., a 0

R x ,a ,a ,..., a 0

R x ,a ,a ,..., a 0


30

Em outras palavras, este método impõe que o resíduo seja zero em determinados pontos do

domínio . O número de pontos n é o mesmo número de coeficientes da solu-

ção aproximada.

Exemplo 4.1

Vamos considerar nossa já conhecida equação diferencial:

Sujeita as condições de contorno

e

Uma solução aproximada para esta equação pode ser:

Observe que satisfaz as condições de contorno no problema.

Usando o método de colocação pontual, nos vamos resolver o problema de forma que o

resíduo R(x) seja zero para , assim:

Logo,

Impondo a condição de R(x) = 0

Deste modo,

1 2 nx , x ,..., x

22

2

d h10x 0

dx 0 x 1

h(0) 0 h(1) 0

1 1h (x) a x 1 x

1h (x)

1x 0,5

22

2

21 1 1

1 1

2

12

d h10x R(x)

dx

h a x(1 x) a x a x

dha 2a x

dx

d h2a

dx

21 1R(x,a ) 2a 10x

1

21

1

R(x,a ) 0

2a 10x 0

a 1,25

1h (x) 1,25x 1 x

Profª Celme Torres

31

Exemplo 4.2

Considerando a mesma equação diferencial, com o objetivo de verificar o acréscimo de

precisão da solução ao se aumentar o número de termos da função de aproximação, va-

mos agora considerar uma solução aproximada do tipo:

Assim,

Logo,

Simplificando,

Como agora existem temos dois coeficientes a determinar, vamos impor que

em dois pontos do domínio do problema. Podemos escolher esses pontos

aleatoriamente.

Escolhendo e , temos:

Resolvendo o sistema acima para a1 e a2,

Assim ,

22 1 2h (x) a x 1 x a x 1 x

22

2

2 2 2 31 2 1 1 2 2

21 1 2 2

2

1 2 22

d h10x R(x)

dx

h a x 1 x a x 1 x a x a x a x a x

dha 2a x 2a x 3a x

dx

d h2a 2a 6a x

dx

21 2 1 2 2R(x,a ,a ) 2a 2a 6a x 10x

21 2 1 2 2R(x,a ,a ) a a 3a x 5x

1 2R(x,a ,a ) 0

1x 0,3 2x 0,6

21 2 2

21 2 2

a a 3a (0,3) 5(0,3) 0

a a 3a (0,6) 5(0,6) 0

1

2

a 0,6

a 1,5

22h (x) 0,6x 1 x 1,5x 1 x


32

Comparação gráfica entre e e a solução exata.

Figura 5.2 – Comparação gráfica entre as soluções encontradas utilizando o Método de Colocação Pontual e a solução exata.

4.1.2. Método dos subdomínios

Neste caso, todas as funções ponderadoras são a unidade . Isto equivale a reque-

rer que a integral do resíduo desapareça em determinados intervalos do domínio.

O número desses intervalos é igual ao número de coeficientes indeterminados da solução

aproximada.

Vamos aplicar esse método à equação diferencial dos exemplos anteriores.

Exemplo 4.3

Vamos considerar nossa já conhecida equação diferencial:


1h (x) 2h (x)

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0,35

0,40

0,45

0 0,2 0,4 0,6 0,8 1 1,2

Solução h1

Solução h2

Solução Exata

iW (x) 1

22

2

d h10x 0

dx 0 x 1

Profª Celme Torres

33

e

Considerando a solução aproximada,

Como há apenas um coeficiente a ser determinado, o intervalo de integração será o pró-

prio domínio do problema .

Vimos que,

Assim,

Portanto,

No caso de escolhermos uma função do tipo:

Vimos que,

Como agora existem dois coeficientes a serem determinados, usaremos dois intervalos

e .

h(0) 0 h(1) 0

1 1h (x) a x 1 x

0 x 1

21 1R(x,a ) 2a 10x

1

1

0

R(x,a )dx 0

13

1

0

102a x x 0

3

1

5a

3

1

5h (x) x 1 x

3

22 1 2h (x) a x 1 x a x 1 x

21 2 1 2 2R(x,a ,a ) 2a 2a 6a x 10x

0 x 0,5 0,5 x 1,0

0,5

1 2

0

1,0

1 2

0,5

R(x,a ,a )dx 0

R(x,a ,a )dx 0


34

O que dá origem ao seguinte sistema de equações:

Resolvendo o sistema, temos:

e

Assim,

Comparação gráfica entre e e a solução exata utilizando o método dos

subdomínios.

Figura 5.3 – Comparação gráfica entre as soluções encontradas utilizando o Método dos Subdomínios e a solução exata.

4.1.3. Método dos mínimos quadrados

Este método requer que a integral do quadrado do resíduo seja minimizado. A integral do

quadrado do resíduo é dada por:

(4.8)

1 2

1 2

a 0,25a 0,4167

a 1,25a 2,917

1

5a

6 2

5a

3

22

5 5h (x) x 1 x x 1 x

6 3

1h (x) 2h (x)

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0,35

0,40

0,45

0,50

0 0,2 0,4 0,6 0,8 1 1,2

Solução h1

Solução h2

Solução Exata

2

1 2 n

D

Q R(x,a ,a ,...,a ) dx

Profª Celme Torres

35

Minimizando essa integral, temos:

(4.9)

Ou,

(4.10)

Ou seja, por esse método, as funções ponderadoras são:

(4.11)

2

1 2 n

D

1 1

R(x,a ,a ,...,a ) dxQ

0a a

2

1 2 n

D

2 2

2

1 2 n

D

n n

R(x,a ,a ,..., a ) dxQ

0a a

R(x,a ,a ,..., a ) dxQ

0a a

2

1D

R dx 0a

1D

RRdx 0

a

2

2D

R dx 0a

2D

RRdx 0

a

2

nD

R dx 0a

nD

RRdx 0

a

1

1

2

2

n

n

RW (x)

a

RW (x)

a

RW (x)

a


36

Exemplo 4.4

Aplicando o método dos mínimos quadrados a nossa já conhecida equação diferencial:


e

Usado a função aproximada,

Vimos que,

Logo,

Portanto,

Que é o mesmo resultado obtido usando o método dos subdomínios. Desta forma, se usar-

mos uma função de aproximação com dois coeficientes, obteremos o mesmo resultado do

método dos subdomínios.

É importante ressaltar que essa igualdade dos resultados entre os métodos acontece para a

equação diferencial usada nos problemas e não se estende genericamente para qualquer

tipo de equação diferencial.

4.1.4. Método de Galerkin

Neste método, as funções ponderadoras são iguais as funções de aproximação

. Assim,

22

2

d h10x 0

dx 0 x 1

h(0) 0 h(1) 0

1 1h (x) a x 1 x

21 1R(x,a ) 2a 10x

1

1

RW (x) 2

a

1 12

1

10 0

RRdx a 5x dx 0

a

1

5a

3

iW (x)

iN (x)

Profª Celme Torres

37

(4.12)

Sendo,

o conjunto , com . Podendo

o mesmo ser escrito ma forma:

(4.13)

Exemplo 4.5

Vamos aplicar o método de Galerkin para resolver a equação diferencial:


e

Usado como função de aproximação,

Vimos que,

e,

Portanto,

Logo,

Usando agora uma função de aproximação com dois coeficientes, temos:

E repetindo o mesmo procedimento, teremos: e,

i i

D D

W (x)R(x)dx N (x)R(x)dx 0

1 2 n[N] N , N ,..., N i

D

N (x)R(x)dx 0 i 1,2,...,n

T

D

N R(x)dx 0

22

2

d h10x 0

dx 0 x 1

h(0) 0 h(1) 0

1 1h (x) a x 1 x

1N (x) x(1 x) 21R(x) 2a 10x

12

1

0

x(1 x)( 2a 10x )dx 0

1

3a

2

22 1 2h (x) a x 1 x a x 1 x

1

3a

2 2

5a

3


38

4.1.5. Exercícios

(1) Aplique o método de colocação pontual e o método dos mínimos quadrados para

encontrar a solução aproximada da equação diferencial no intervalo de 0 x 20 .

Utilize um programa gráfico (Maple ou Excel) para representar graficamente as

soluções indicando qual mais se aproxima da solução exata da EDO.

y'' 0,25x 0 y(0) 2 e y'(0) 2

(2) Aplique o método de Galerkin e o método dos subdomínios utilizando dois

coeficientes para encontrar a solução aproximada das equações diferenciais no

intervalo de 0 x 1 . Utilize um programa gráfico (Maple ou Excel) para representar

graficamente as soluções indicando qual mais se aproxima da solução exata da EDO.

2

2

d y dy2x

dx dx y(0) 1 e y'(0) 3

(3) Encontre a solução aproximada da equação diferencial no intervalo indicando

utilizando os métodos dos resíduos ponderados, colocação pontual, mínimos

quadrados, subdomínios e o método de Galerkin. Represente graficamente as soluções

aproximadas comparando com a solução exata da EDO.

(a) 216y'' 8y ' 0,150x 0 x 10 y(0) 2 e y'(0) 1

(b) ty'' 3y ' 4y 2e 0 x 1 y(0) 1 e y'(0) 0

(c) 2

2

d y4y 3sin 2t

dt 0 x 2 y(0) 2 e y'(0) 1

(4) Considere a equação diferencial que rege a distribuição da temperatura T(x), em uma

barra trocando calor com o ambiente, cuja temperatura é de 25ºC.

Sujeita as condições de contorno:

T(0) 150

dT(2) 0

dx

Proponha uma solução com um termo e outra com dois termos e ache os valores dos coefi-

cientes dessas soluções usando o método de colocação pontual, o método dos subdomínios,

dos mínimos quadrados e o método de Galerkin.

Compare a graficamente a solução com dois termos de cada método com a solução analíti-

ca e determine qual método mais preciso, neste caso.

2

2

d T3(T 25) 0

dx 0 x 2

Profª Celme Torres

39

5. Introdução aos Métodos Variacionais

A análise de problemas de engenharia envolve o isolamento de um elemento diferencial da

região ou domínio onde o fenômeno está ocorrendo.

MODELO MATEMÁTICO RESULTANTE EQUAÇÃO DIFERENCIAL

Tipos de solução de uma Equação Diferencial

1. Solução analítica ou Exata

2. Solução aproximada ou solução numérica.

A solução analítica é a ideal, pois ela se constitui de uma expressão (equação) que satisfaz

a equação diferencial e suas condições de contorno, portanto para obtermos o valor exato

da variável dependente da equação diferencial em qualquer ponto do domínio do problema,

basta substituirmos as coordenadas desse ponto na equação solução.

5.1.1. Cálculo variacional

O cálculo variacional é uma parte da análise matemática que trata de problemas de máxi-

mos e mínimos de tipos especiais de funções, chamadas FUNCIONAIS.

Vários são os motivos que tornam convenientes o estudo do cálculo variacional. O princi-

pal deles é que este se constitui em uma alternativa que permite definir com mais clareza

os conceitos básicos que correspondem a problemas de engenharia governados por siste-

mas de equações diferenciais.

O cálculo variacional é especialmente aplicável a problemas de mecânica dos meios contí-

nuos. Através do uso de princípios de energia, é possível formular os problemas na forma

variacional.

5.1.2. Máximo e mínimo de funções

Suponha a seguinte função apresentada no gráfico.


40

O domínio de F(x) é a bx x x . A relação funcional implica que, para cada valor da va-

riável x, um único valor será obtido para a função F(x), dentro do domínio.

O menor valor de F(x) dentro do domínio é o mínimo absoluto da função F(xa) e o máxi-

mo absoluto é dado por F(xb). Além disso, podem existir valores mínimos e máximos rela-

tivos.

Por exemplo; a função F(x) terá mínimo relativo ou um máximo relativo, para um certo

valor da variável independente x, quando todos os valores de F(x) na vizinhança infinite-

simal de x: x x , com x tão pequeno quanto se queira, seja maior ou menor, respecti-

vamente, que F(x). Para estudar as condições que determinam os valores extremos relat i-

vos, é conveniente expandir F(x) em uma Série de Taylor em torno de x.

2 32 3

2 3x

dF 1 d F 1 d FF(x x) F(x) x x x ...

dx 2! 3!dx dx (4.1)

O incremento F será dado por:

2 32 3

2 3x

dF 1 d F 1 d FF(x x) F(x) x x x ...

dx 2! 3!dx dx (4.2)

Observe que haverá um mínimo relativo quando:

F(x) é sempre menor que F(x x) , não interessando o valor de x , para F 0

F positivo;

Para F 0 F negativo máximo relativo.

F(x)

x

xb x2 x1 xa

Profª Celme Torres

41

Como x é infinitamente pequeno, x é sempre maior que 2x e que

3x e assim su-

cessivamente.

Portanto,

Se x

dF0

dx , o sinal de F será indefinido, pois depende do sinal de x ;

Se x

dF0

dx , não teremos ponto de mínimo e de máximo.

Então, para se garantir que a função terá pontos de mínimo e de máximo relativo a primei-

ra derivada obrigatoriamente será zero x

dF0

dx .

Os pontos que satisfazem essa condição são chamados de pontos estacionários. Neste ca-

so, como 2x é sempre positivo, a característica de um ponto estacionário dependerá do

sinal da 2ª derivada de F(x): 2

2

x

d F

dx.

Se 2

2

x

d F0

dx mínimo relativo

Se 2

2

x

d F0

dx máximo relativo

Se a derivada segunda for nula, haverá um ponto neutro. Neste caso, devem ser analisados

os sinais das derivadas de ordem superior, para se determinar a natureza do ponto estacio-

nário. Se todas as derivadas até a ordem n 1 forem nulas, a natureza do ponto estacioná-

rio dependerá do sinal da enésima derivada.

Se n for par, será aplicado para o último termo da Série de Taylor o mesmo critério que se

usa para a 2ª derivada.

Se n for impar, o termo n

n

n

x

1 d Fx

n! dx , vaiará com o sinal de x e não se terá um máximo

ou mínimo relativo, e sim um ponto indefinido.


42

No caso da 2ª derivada ser nula, tais pontos são denominados de pontos de sela ou pontos

de inflexão.

Uma análise prévia pode ser estendida a funções de mais de uma variável independente.

Por exemplo, funções do tipo: F(x, y) , com duas variáveis independentes. Ao se expandir

esta função em uma Série de Taylor, temos:

2 2 22 2

2 2x,y x,y x,y x,y x,y

F F(x x, y y) F(x, y)

F F 1 F F FF x y x 2 x y y ...

x x 2! x yx y

(4.4)

Como visto anteriormente, haverá um ponto estacionário quando: x,y x,y

F F0

x y

E a característica do ponto estacionário será determinada, na forma:

2 2

22 2 2

x,y x,y2 2

2 2 2 22x,y x,y x,y

2

x,y x,y

F F

x yx xF F Fx 2 x y y x y

yx yx y F Fy

x y y

(4.5)

A fórmula acima é quadrática e o sinal depende da natureza da matriz das derivadas de 2ª

ordem, assim:

Quando a matriz é positiva definida, com autovalores positivos, tem-se um ponto

de mínimo.

Quando a matriz é negativa definida, com autovalores negativos, tem-se um ponto

de máximo.

Exemplo 5.1

Considere a função

3 23f (x) x x 6x

2

Os pontos estacionários dessa função são dados por:

Profª Celme Torres

43

2df3x 3x 6 0

dx

A qual será definida para 1

2

x 1

x 2

Ao se verificar a segunda derivada : 2

2

d f6x 3 0

dx

Pode verificar que f(1) é um ponto de mínimo relativo e f(-2) é um ponto de máximo relati-

vo.

Exemplo 5.2

Considere a função

2 2f (x, y) x 2xy 4y 2x 8y

Haverá um ponto estacionário quando:

f2x 2y 2 0

x

f2x 8y 8 0

y

Reescrevendo temos:

2x 2 2

2x 8y 8

A matriz da 2ª derivada é dada por:

2 2

2

2 2

2

F F

x y 2 2x

2 0F F

x y y

Cujos autovalores são positivos ponto estacionário é um ponto de mínimo.

5.1.3. Funcional

A seção anterior tratava de determinar pontos estacionários de funções. O problema do cál-

culo variacional também se refere à determinação de pontos estacionários, mas definidos

para funções especiais, chamadas funcionais.


44

Por exemplo, seja a seguinte função: f (x, y,z,...)

Ao se dar valores às variáveis independentes, um valor numérico é obtido para a função f.

Determina-se assim, um grupo de valores de x, y, z, ..., para o qual define-se um ponto es-

tacionário para f.

O FUNCIONAL é um tipo especial de função, no qual as variáveis independentes não são

simples números, mas funções desconhecidas. Ou seja, um funcional é um grupo (ou uma

função) que depende de outras funções. Funcionais são chamados de F ou .

1 2

20

1

x xx

0

3df 8d fF (2f )dx

dx dx

F G(g,g ,g , x)dx

Onde, x

dff

dx e

2

xx 2

d fg

dx

Tanto F quanto G são considerados FUNCIONAIS. Quando se adota diferentes valores

para as funções desconhecidas f, se obtém diferentes valores para os funcionais.

Em geral um funcional pode depender de várias funções, as quais podem depender de uma

ou mais variáveis independentes.

5.1.4. Funcional de uma função

Considere o funcional no domínio

a bx x x

b

a

x

x

x

F G(g,g , x)dx

Condições de contorno para g: a a

b b

g(x ) g

g(x ) g

Profª Celme Torres

45

Considere agora que g é uma função que faz o funcional F ser estacionário.

Definindo um grupo de funções:

g( ,x) g(x) (x)

onde, é um escalar infinitesimal e (x) é uma função conhecida com as seguintes condi-

ções:

a

b

(x ) 0

(x ) 0

Podendo escrever então:

f g g

onde, g (x)

Sendo o símbolo variacional e significa a 1ª variação de g.

No cálculo variacional α corresponde a x da Série de Taylor. Dessa forma, o funcional é

dado por:

b

a

x

xxx

G GF g g dx

g g

(4.6)

5.1.5. Propriedades de um Funcional:

d( f ) (df )

(f g) f g

xa xb

𝑛(𝑥)

𝑓(𝑥)

𝑔 𝑥 = 𝑓 𝑥 + 𝑛(𝑥)

𝛼𝑛(𝑥)


46

(f.g) f g fg

2

f f g fg

g g

n n 1(f ) nf f

Pretende-se agora determinar as condições para que o funcional seja estacionário através

do estudo da variação de f f . Isto é semelhante ao caso de funções onde a existência de

um ponto estacionário é determinada para um certo valor de x através do estudo variacio-

nal de f(x) no ambiente x x .

Se considerarmos,

b

a

x

x x

x

(g) F(f , f , x)dx (4.8)

Onde f e são conhecidos, portanto o funcional só depende de α.

Expandindo em torno de f através de uma Série de Taylor, temos:

2 3

2 3

2 30 0 0

d 1 d 1 d(g) (f ) ...

d 2! 3!d d

(4.9)

2 3

2 3

2 30 0 0

d 1 d 1 d(g) (f ) ...

d 2! 3!d d

(4.10)

Dessa forma é possível definir a variação do funcional , de acordo com a ordem da equa-

ção diferencial:

1ª ordem 0

d

d

2ª ordem 2

2 2

2

0

d

d

Profª Celme Torres

47

3ª ordem 3

3 3

3

0

d

d

A equação para variacional total pode ser escrita na forma:

2 31 1...

2! 3! (4.11)

De maneira semelhante aos caso de funções, a natureza de ponto estacionário, isto é, se é

um ponto de máximo, de mínimo ou de inflexão, dependerá do sinal da derivada de 2ª

ordem.

5.1.6. Determinação do Funcional - Equação de Euler

Considere uma integral da forma:

(4.12)

Com,

(4.13)

Onde o integrando é uma dada função de x, h e hx. Tanto I(h) como F são

chamados de funcional. A palavra funcional significa função de funções.

Considere o funcional . Uma mudança de em h, onde é uma constante e

w é uma função, é chamada de variação de h e é denominada:

(4.14)

Onde é chamado de símbolo variacional.

A primeira variação de F em h é definida por:

(4.15)

b

x

a

I(h) F(x,h,h )dx

h h(x) x

dhh

dx

xF(x,h,h )

xF(x,h,h ) w

h w

xx

F FF h h

h h


48

O símbolo age como um operador diferencial com respeito às variáveis dependentes,

com a seguinte propriedade:

(4.16)

Considere agora o problema de se determinar o mínimo do funcional:

(4.17)

A condição necessária para minimizar esse funcional é que a primeira variação de I deve

ser zero, ou:

(4.18)

Usando a condição para o funcional ser estacionário, temos:

(4.19)

Podemos ainda escrever:

(4.20)

Assim,

(4.21)

Usando integral por partes ( ) com, e

(4.22)

b b

a a

dh d( h)

dx dx

h(x)dx h(x)dx

b

x

a

I(h) F(x,h,h )dx

b

x

a

I F(x,h,h )dx 0

b

xx

a

F FI h h dx 0

h h

x

dh d( h)h

dx dx

b

xa

F F d( h)I h dx 0

h h dx

udv uv vdu v h x

Fu

h

b b

x x aa

F d F FI h h dx h

h dx h h

Profª Celme Torres

49

(4.23)

O que só pode ser verdade para qualquer variação de se:

(4.24)

com, Essa expressão é gerada pela condição de contorno do problema.

Este resultado é uma equação diferencial com h e hx como pseudovariáveis independentes.

A equação 4.24 é chamada de Equação de Euler ou Euler-Lagrange.

É importante salientar que, a função h(x) que minimiza o funcional I é a mesma função que

satisfaz a equação diferencial de Euler.

Exemplo 5.3

Determine o funcional da equação diferencial:


e

Como h é uma constante nos extremos do intervalo , a variação de é zero

nestes pontos. Este tipo de condição de contorno é denominada essencial ou geométrica.

Neste caso,

b b

x x aa

F d F FI h dx h 0

h dx h h

h

x

F d F0

h dx h

b

x a

F0

h

22

2

d h10x 0

dx

h(0) 0 h(1) 0

0 x 1 h

1

x 0

Fh 0

h


50

Determinação do funcional F usando a equação de Euler:

Assim, usando a equação (10), temos:

Integrando a 2ª integral por partes e considerando e , temos:

(E.1)

Analisando novamente a 2ª integral da expressão acima, usando d( h) (dh) , a seguinte

identidade é considerada:

Se, na equação (E.1) assumirmos que,

, teremos:

Como a variação de x2 é zero, e usando , temos:

Portanto, o funcional F associado à equação diferencial do problema é:

22 2

2

F d F d h d dh10x 10x

h dx x dx dxdx

1

2

0

1 1

2

0 0

d dhI 10x hdx 0

dx dx

d dhI 10x hdx hdx 0

dx dx

u h dh

vdx

1 11

2

00 0

d hdh dhI 10x hdx h dx 0

dx dx dx

1 1 1 2

0 0 0

d hdh dh dh 1 dhI dx dx dx

dx dx dx dx 2 dx

1

0

dhh 0

dx

1 22

0

1 dhI 10x h dx

2 dx

b b

a a

Fdx Fdx

1 22

0

1 dhI 10x h dx

2 dx

Profª Celme Torres

51

É bom enfatizar que a função h(x) que minimiza o funcional acima é também a solução da

ED que rege o problema. O funcional acima é chamado de forma variacional ou forma

fraca da ED. Lembrando que I não é uma solução, apenas o funcional da ED.

5.1.7. Exercícios

(1) Determine o funcional associado às equações diferenciais sujeita as condições de con-

torno apresentadas.

a) sendo,

x x

h(0) 0 h(5) 100

h (0) 1 h (5) 0

b) 2

2 2

2

d f9f 5x

dx onde, 0 x 5 sendo,

f (0) 9 e f (5) 0

c) 3

3

d f5f 0

dx onde, 0 x 4 sendo,

df (0)0

dx e

f (4) 0

(2) Determine a equação de Euler associada ao funcional

Sabendo que:

1 22

0

1 dhI 10x h dx

2 dx

4 2

4 2

d h d h8 4h 10

dx dx 0 x 5

b

x xx

a

I F(x,h,h ,h )dx 0

b

x xx a

F d Fh 0

h dx h


52

b

xx a

Fh 0

h

Profª Celme Torres

53

6. Métodos Variacionais

6.1.1. Método de Ritz

O primeiro passo do método de Ritz é assumir uma solução aproximada para a equação

diferencial na forma de uma série finita:

(6.1)

Onde as constantes ai são chamadas coeficientes de Ritz e as funções Ni(x) são conhecidas

como funções de aproximação.

Como hn é uma solução da equação diferencial, hn deve satisfazer as condições de contorno

essenciais. Qualquer condição de contorno natural já está incluída na formulação variacio-

nal.

O passo seguinte no Método de Ritz usa o princípio de que a função que minimiza o

funcional também é a solução da equação diferencial associada ao funcional, assim

como:

(6.2)

Após a substituição de em , este funcional se torna função dos parâmetros

ou e a condição necessária para minimizar I será:

(6.3)

n

n i i

i 1

h a N (x)

h(x)

I(x)

n

n i i

i 1

h a N (x)

nh I(h)

1 2 na ,a ,...,a 1 2 nI(a ,a ,...,a )

1

2

n

I0

a

I0

a

I0

a


54

O que resultará em um sistema de n equações com n coeficientes a determinar

.

Uma vez conhecidos os valores desses coeficientes (conhecidos como coeficientes de

Ritz), então a solução aproximada será conhecida.

Exemplo 6.1

Considere a equação diferencial:


e

Como vimos o funcional associado a essa equação diferencial é:

Usando o método de Ritz, vamos assumir a seguinte solução aproximada:

Observe que h1 satisfaz as condições essenciais e de contorno, portanto:

i(a ,i 1,2,...,n)

nh

22

2

d h10x 0

dx 0 x 1

h(0) 0 h(1) 0

1 22

0

1 dhI 10x h dx

2 dx

1 1h a x(1 x)

1 1

21 1

11

h a x(1 x)

h a (1x x )

dha (1 2x)

dx

Profª Celme Torres

55

Logo,

Dessa forma:

Minimizando então I:

Portanto,

Vamos agora usar a seguinte a solução aproximada:

Usando os dois resultados acima, obtemos:

Minimizando agora I com respeito a e obteremos o seguinte sistema de equações:

2 221 1

1

dh a1a (1 2x)

2 dx 2

12

23 4 11

0

aI 10a x x 1 2x dx

2

1 14 52 2 3

1 1

0 0

21 1

x x 4I 10a a x 2x x

4 5 3

a aI

2 6

1

I0

a

1a10

2 3 1

3a

2

1

3h x 1 x

2

22 1 2

221 2

h a x(1 x) a x (1 x)

dha (1 2x) a (2x 3x )

dx

2 21 1 2 1 2 2a a a a a a

I6 2 3 6 15

1a 2a


56

Resolvendo o sistema acima para e , obtemos:

e,

Dessa forma,

A seguir faremos uma comparação gráfica, usando o Maple, entre e e a so-

lução exata:

Figura 6.1 – Comparação gráfica entre as soluções encontradas utilizando o Método de Ritz e a solução exata.

Verificamos que a função aproximada encontrada ao se usar o Método de Ritz através do

método variacional coincidiu com a solução obtida aplicando o Método de Galerkin. Tal

fato não foi uma coincidência. Considerando a equação resultante da minimização do fun-

cional I,

1 2

1

1 2

2

a aI 10

a 2 3 6

a 2aI 10

a 3 6 15

1a 2a

1

3a

2 2

5a

3

22

2 5h x 1 x x 1 x

3 3

1h (x) 2h (x)

410h(x) x x

12

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0,35

0,40

0,45

0 0,2 0,4 0,6 0,8 1 1,2

Solução h1

Solução h2

Solução Exata

Profª Celme Torres

57

Usando a solução aproximada,

Onde,

e,

Neste caso, como as funções de aproximação são invariáveis,

Substituindo a expressão acima na equação do funcional, temos:

Como, não depende de x,

Considerando que, é a equação diferencial associada ao funcional I, e que

a solução aproximada foi usada,

Teremos,

Para a expressão acima ser verdadeira para qualquer valor de ,

b b

x x aa

F d F FI h dx h 0

h dx h h

nT

m i i

i 1

h h (x) a N (x) [a][N]

1 2 n[a] [a ,a ,...,a ]

1 2 n[N] [N , N ,..., N ]

T[N]

Th [a][N]

b bT T

x x aa

F d F FI [a][N] dx [a][N] 0

h dx h h

[a]

b bT T

x x aa

F d F FI [a] [N] dx [a][N]

h dx h h

x

F d F

h dx h

mh (x)

m x

F d Fdx R(x)

h dx h

b bT T

x aa

FI [a] [N] R(x)dx [N]

h

[a]


58

O primeiro termo da expressão acima é o princípio básico do Método de Galerkin. O se-

gundo termo refere-se às condições de contorno naturais do problema.

Podemos então concluir que havendo um princípio variacional para o problema, a solução

através do uso do método de Galerkin coincide com a solução do método variacional. Essa

é a razão pelo qual o Método de Galerkin é o escolhido como fundamento para o Método

dos Elementos Finitos (MEF).

6.1.2. Exercícios

(1) Considere a seguinte equação diferencial

Sujeitas as condições de contorno:

a) Ache o funcional associado a essa equação diferencial

b) Pode a seguinte equação 1

xT (x) sin

10

ser usada na solução aproximada? Se puder,

usando o Método de Ritz, ache a solução aproximada usando a função acima e compare

graficamente com a solução exata.

c) Usando a função aproximada , ache os coeficientes e

.

d) Faça um gráfico comparando as soluções (b) e (c) com a solução exata.

(2) Dada a seguinte equação diferencial


b bT T

x aa

F[N] R(x)dx [N] 0

h

2

2

d T100 0

dx 0 x 10

T(0) 0

T(10) 0

2 1 2

x 3 xT (x) a sin a sin

10 10

1a

2a

2

2

d y6y 10x

dx 0 x 2

Profª Celme Torres

59

y(0) 1

y(2) 0

Encontre o funcional associado à equação diferencial e encontre uma solução aproxi-

mada usando o Método de Ritz.

6.1.3. Método de Rayleigh-Ritz

No método de Ritz apresentado anteriormente, as funções de Ni(x) eram válidas em todo o

domínio do problema. No método de Rayleigh-Ritz, essas funções de aproximação são vá-

lidas apenas em trechos do problema. Portanto, o domínio do problema é dividido em ele-

mentos e as funções de aproximação são específicas para esses elementos.

Vamos analisar agora como encontrar essas funções de aproximação Ni(x), também co-

nhecidas como funções de forma.

Considere um problema cujo domínio é o intervalo que será dividido em M in-

tervalos ou elementos finitos. Nos limites de cada elemento, há dois nós – denominados

genericamente de i e j. Como, nesse caso, cada elemento e tem dois nós, uma interpolação

polinomial linear deve ser usada como função de interpolação no elemento, portanto:

De acordo com a Figura 5.4, podemos constatar que:

(6.4)

Portanto, (6.5)

a x b

eh (x) cx d

ei i

ej j

x x h h

x x h h

i i

j j

h cx d

h cx d


60

Figura 6.2 – Divisão de um domínio qualquer em elementos lineares.

Resolvendo o sistema para achar c e d, temos:

e (6.6)

Assim,

(6.7)

Onde e são as funções de forma.

Logo,

(6.8)

h (x)e

hi

hj

xjxi

h(x)

x

x = a x = b

1 2 3 i j ne1 e

2 em

e

j i

j i

h hc

x x

i j j i

j i

h x h xd

x x

je ii j i i j j

j i j i

x x x xh (x) h h N (x)h N (x)h

x x x x

iN (x) jN (x)

j i

j i j i

x x x x[N]

x x x x

Profª Celme Torres

61

Figura 6.3 – Interpolação linear de um elemento.

Exemplo 6.2

Considere a já conhecida equação diferencial:


e

Como vimos o funcional associado a esse equação diferencial é:

Se dividirmos o domínio em vários trechos (elementos), temos:

Sendo,

e,

Minimizando o funcional do elemento, temos:

h =N (x)h + e

i i N (x)h j j

hi

hj

x

i je

1,0 1,0N (x)i N (x)j

22

2

d h10x 0

dx 0 x 1

h(0) 0 h(1) 0

1 22

0

1 dhI 10x h dx

2 dx

jj

ii

x 2x ee e 2 e

xx

dh 1 dhI h 10x h dx

dx 2 dx

ne

e 1

I I

eeh N a

ie

j

aa

a


62

Mas como, podemos reescrever o termo da seguinte forma:

Sendo a matriz transposta de .

Da mesma forma podemos simplificar o segundo termo da integral, onde:

Portanto,

Ou ainda,

O que pode ser escrito de forma simplificada:

Sendo, [K]e a matriz de rigidez do elemento; f

e é conhecido como vetor de carregamento

nodal do elemento e ae é o vetor dos valores aproximados a função no nó do elemento.

Depois que as matrizes [K]e e os vetores f

e e a

e de todos os elementos serem conheci-

dos, encontra-se o sistema global do problema, na forma:

jj

ii

x 2xe e e e2

e e e ex

x

dI dh dh dh 1 d dh10x dx 0

dx 2 dxd a d a d a d a

eeh N a

eT

e

dhN

d a

T

N N

2e e e e

e e e

e e

e

e e

e

dh 1 d dh 1 dh d dh(2)

2 dx 2 dx dxd a d a d a

d dh dh

dx dxd a

d N d Nd a a

dx dxd a

Ted N d N

adx dx

j jj

ii i

x xx TeT T e2

ex

x x

d N d NdI dhN 10x N dx a dx 0

dx dx dxd a

j jj

ii i

x xxTe T T2

xx x

d N d N dhdx a N 10x N dx

dx dx dx

e e e

K a f

K a f

Profª Celme Torres

63

Agora as condições de contorno podem ser impostas. Voltando ao nosso exemplo. Se divi-

dirmos o domínio em dois elementos, temos dois elementos e três nós, conforme a Figura

5.6.

Figura 6.4 – Elementos lineares

Para cada elemento, temos:

Para os dois elementos do problema:

Analisando o termo,

Analisando agora a integral,

j

i

x Te

j ix

1 1d N d N 1K dx

1 1dx dx x x

1 2 2 2K K

2 2

j

j

ii

xxT Te 2

xx

dhf N 10x N dx

dx

j

i

x e ei jT i ij i

j ij jx

e ej i

ei

ej

N (x ) N (x )dh (x ) dh (x )dhN

N (x )N (x )dx dx dx

dh (x )0 1 dh (x )

1 0dx dx

dh (x )

dx

dh (x )

dx

e1 e

2

1 2 3

0,5 0,5


64

Assim, para o elemento 1, temos:

Para o elemento 2,

Fazendo a montagem do sistema global:

Levando em consideração a equação da continuidade, podemos impor, que:

Dessa forma o sistema global se reduz a forma:

j

i

x 4 3 4j j i iT2

4 4j i j i j ix

x 4x x 3x1010x N dx

10 x x 3x 4x x x

11

1

12

dh (x )0,104

dxf

dh (x )0,313

dx

22

2

23

dh (x )1,146

dxf

dh (x )1,771

dx

11

1 1 22 2

2

3 23

dh (x )0,104

dxh2 2 0

dh (x ) dh (x )2 2 2 2 h 1,459

dx dx0 2 2 h

dh (x )1,771

dx

1 21 2dh (x ) dh (x )

dx dx

Profª Celme Torres

65

Usando a 2ª linha do sistema encontramos o valor para h2

A solução gráfica apresenta-se da seguinte forma:

Figura 6.5 – Comparação gráfica entre as soluções encontradas utilizando o Método de Rayleigh-Ritz e a solução exata.

11

1

2

23 3

dh (x )0,104

h2 2 0 dx

2 4 2 h 1,459

0 2 2 h dh (x )1,771

dx

2h 0,365

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0,35

0,40

0,45

0 0,25 0,5 0,75 1 1,25

Solução Exata

Rayleigh-Ritz


66

7. Métodos das Diferenças Finitas

Como mencionado anteriormente, a idéia básica do método das diferenças finitas (MDF)

consiste na aproximação das derivadas parciais de uma equação diferencial através de fór-

mulas de diferenças, e satisfação dessa equação aproximada de diferenças em determina-

dos pontos do domínio (região). Isto conduz a um sistema de equações algébricas cuja so-

lução fornece o valor da função desejada nos pontos pré-determinados.

Por definição, a derivada de uma função u(x) num ponto é dada por:

(7.1)

De forma aproximada, pode-se tomar:

(7.2)

Onde h é um incremento pequeno, porém finito. Esta aproximação é chamada de diferença

progressiva, já que no seu cálculo aparecem os valores de u nos pontos e

. De maneira análoga, pode-se definir uma aproximação de diferenças regressi-

vas na forma:

(7.3)

E um esquema de diferença central:

(7.4)

As aproximações acima podem ser visualizadas na Figura 6.1.

1x x

i

i i

h 0x x

u(x h) u(x )dulim

dx h

i

i i

x x

u(x h) u(x )du

dx h

1x x

1x x h

i

i i

x x

u(x ) u(x h)du

dx h

i

i i

x x

u(x h) u(x h)du

dx 2h

Profª Celme Torres

67

Uma forma alternativa de se obter fórmulas aproximadas de diferenças finitas é através de

séries de Taylor. Desta maneira, é possível estimar o erro cometido em cada tipo de apro-

ximação.

A expansão em série de Taylor do valor de u em , em torno do valor de u em

, é dada por:

(7.5)

Figura 7.1 – Interpretação gráfica das aproximações por diferenças finitas

Que pode ser reescrita na forma:

(7.6)

Tornando-se:

(7.7)

ix x h

ix x

i i i i

2 2 3 3 n n

i i 2 3 nx x x x x x x x

du h d u h d u h d ux(x h) u(x ) h ...

dx 2 6 2ndx dx dx

x +hixi

y

xx - hi

h h

y= u(x)

tg

i i i i

2 2 3 n ni i

2 3 nx x x x x x x x

u(x h) u(x )du h d u h d u h d u...

dx h 2 6 2ndx dx dx

i

i i

x x

u(x h) u(x )du

dx h


68

Que é a aproximação de diferença progressiva. Pode-se notar que foram desprezados os

termos relativos à derivada de segunda ordem e de ordem superior. Como h é um valor pe-

queno, maior termo desprezado é igual a uma constante multiplicada por h, ou seja, é da

ordem de h.

Analogamente,

(7.8)

Que dá a aproximação:

(7.9)

Para se obter a aproximação de diferença central, subtrai-se (5.8) de (5.5), chegando-se a:

(7.10)

Ou,

(7.11)

O erro cometido pelo esquema de diferença central é menor que o de diferença progressiva

ou regressiva.

7.1.1. Derivadas de ordem superior

Para se obter aproximações por diferenças finitas para derivadas de ordem superior, pode-

se empregar séries de Taylor ou uma maneira mais intuitiva, através da aplicação repetida

da aproximação (6.2) a (6.4). Tomando-se, por exemplo, a derivada de segunda ordem, po-

de-se somar as expressões (6.5) e (6.8), obtendo-se:

(7.12)

i i i i

2 2 3 3 n n

i i 2 3 nx x x x x x x x

du h d u h d u h d ux(x h) u(x ) h ...

dx 2 6 2ndx dx dx

i

i i

x x

u(x ) u(x h)du

dx h

i i i

3 3 n n

i i 3 nx x x x x x

du h d u h d ux(x h) u(x h) 2h ...

dx 3 ndx dx

i

i i

x x

u(x h) u(x h)du

dx 2h

i i i

2 4 4 n n2

i i i 2 4 n

x x x x x x

d u h d u h d uu(x h) u(x h) 2u x h ...

12 3ndx dx dx

Profª Celme Torres

69

Reescrevendo a expressão (7.12), desprezando-se as derivadas de ordem superior, temos:

(7.13)

Que é a aproximação de diferença central para a derivada de segunda ordem.

Outra opção é aproximar a derivada segunda como uma diferença entre derivadas de pri-

meira ordem:

(7.14)

(7.15)

Cujo resultado é igual a equação (7.13). Os dois processos acima podem ser estendidos pa-

ra a determinação de qualquer ordem de derivada.

Exemplo 7.1

Seja o problema descrito pela equação e condições de contorno abaixo:

Resolver a equação diferencial pelo método das diferenças finitas, com aproximação de

diferenças centrais.

A equação aproximada de diferenças centrais é dada por:

i

2i i i

2 2

x x

u(x h) 2u(x ) u(x h)d u

dx h

i i

i

2x x h / 2 x x h / 2

2

x x

du du

dx dxd u

hdx

i

i i i i2

2

x x

u x h u(x ) u x u(x h)

d u h h

hdx

2

2

d uu x 0 (1 x 0)

dx

0 em x 0u

0 em x 1


70

Reescrevendo os termos, temos:

(i)

Adotando-se , o problema irá apresentar 4 incógnitas, correspondentes ao valor

de u nos pontos:

É necessário então, gerar 4 equações, que são obtidas aplicando a equação (i) nesses 4

pontos.

Introduzindo as condições de contorno e arrumando em forma matricial:

Dessa forma, obtemos os valores de u nos pontos determinados:

A solução exata do problema é dada por:

Comparando graficamente a solução dada por diferenças finitas com aproximação central

com a solução exata da ED, temos:

i i ii i2

u x h 2u(x ) u(x h)u(x ) x 0

h

2 2i i i iu x h (h 2)u(x ) u(x h) h x

h 0,2

x 0, 2

x 0, 4

x 0,6

x 0,8

i

i

i

i

x 0,2 0,4u 1,96u(0,2) 0,0u 0,008

x 0,4 0,6u 1,96u(0,4) 0,2u 0,016

x 0,6 0,8u 1,96u(0,6) 0,4u 0,024

x 0,8 1,0u 1,96u(0,8) 0,6u 0,032

1,96 1 0 0 0,2u 0,008

1 1,96 1 0 0,4u 0,016

0 1 1,96 1 0,6u 0,024

0 0 1 1,96 0,8u 0,032

u(0,2) 0,0362

u(0,4) 0,0630

u(0,6) 0,0713

u(0,8) 0,0527

sinxu

sin1 x

Profª Celme Torres

71

Figura 7.2 – Comparação gráfica da solução exata e da solução por diferenças finitas com aproximação central.

Exemplo 7.2

Resolver a Equação Diferencial do exemplo anterior, com a condição de contorno natural,

dada por:

O problema agora possui 5 incógnitas, já que também não se conhece o valor de u no pon-

to . É necessário, então, obter mais uma equação pela aproximação de diferença re-

gressiva da condição de contorno natural.

Substituindo esta aproximação e a condição de contorno em no sistema de equa-

ções, temos:

du1 em x 1

dx

x 1

x 1

du 1,0u 0,8u1 u(1,0) 0,2 u(0,8)

dx 0,2

x 0


72

Uma alternativa à solução acima é a de se aplicar a equação (i) também no ponto ,

envolvendo o valor de u no ponto fictício e utilizar uma aproximação de diferen-

ças centrais para a condição de contorno natural. Nesta caso, o sistema final terá 5 equa-

ções e 5 incógnitas, na forma:

A solução da ED considerando o primeiro sistema através da condição de contorno natu-

ral com diferença finita regressiva é:

Considerando a condição natural com diferenças finitas aproximadas, temos:

A solução exata do problema é dada por:

Comparando graficamente a solução dada por diferenças finitas com aproximação regres-

siva com a solução dada por diferenças centrais e a solução exata da ED:

1,96 1 0 0 0,2u 0,008

1 1,96 1 0 0,4u 0,016

0 1 1,96 1 0,6u 0,024

0 0 1 0,96 0,8u 0,232

x 1

x 1,2

0,2u 0,0081,96 1 0 0 0

0,4u 0,0161 1,96 1 0 0

0,6u 0,0240 1 1,96 1 0

0,8u 0,0320 0 1 1,96 1

1,0u 0,4400 0 0 2 1,96

u(0, 2) 0,4415

u(0, 4) 0,8573

u(0,6) 1, 2229

u(0,8) 1,5155

u(1,0) 1,7155

u(0, 2) 0,5423

u(0, 4) 1,0550

u(0,6) 1,5094

u(0,8) 1,8794

u(1,0) 2,1422

2sinxu

cos1 x

Profª Celme Torres

73

Figura 7.3 – Comparação entre a solução exata e a solução por diferenças finitas com aproximação regressiva e aproximação central.

7.1.2. Exercícios

(1) Seja o problema descrito pela EDO sujeito as condições de contorno apresentadas:

2

2

2

d H10x 0

dx 0 x 1

H(0) 0

H(1) 0

Resolva a EDO pelo método das diferenças finitas com aproximação de diferença central.

(2) Considere a seguinte equação diferencial de transferência de calor:


2

2

d T100 0

dx 0 x 10

T(0) 0

T(10) 0


74

Encontre a solução numérica pelo método das diferenças finitas com aproximação central,

aproximação regressiva e aproximação progressiva. Compare graficamente as soluções a-

proximadas com a solução exata da EDO.

(3) Aplique o método das diferenças finitas com aproximação progressiva para encontrar

a solução aproximada da equação diferencial no intervalo. Mostre graficamente as so-

luções exata e aproximada utilizando um programa gráfico (maple ou Excel).

y' 4 2x y(1) 5 0 x 1 considere h = 0,25

(4) Aplique o método das diferenças finitas com aproximação regressiva para encontrar a

solução aproximada da equação diferencial no intervalo. Mostre graficamente as solu-

ções exata e aproximada utilizando um programa gráfico (maple ou Excel).

y' cos(x) 1 y(0) 1 0 x 120 considere h =

(5) Aplique o método das diferenças finitas com aproximação central para encontrar a so-

lução aproximada da equação diferencial no intervalo. Mostre graficamente as solu-

ções exata e aproximada utilizando um programa gráfico (maple ou Excel).

2y' yx y y(0) 1 0 x 2 considere h =

(6) Aplique o método de diferenças finitas com aproximação central para encontrar a so-

lução aproximada das equações diferenciais nos intervalos definidos.

a) y'' 2y ' y x 0 x 1 y(0) 2 e y'(1) 0 (use h =

0,25)

b) y'' 2y ' x 0 x 6 y(0) 4 e y'(6) 2 (use h =

1,50)

Profª Celme Torres

75

8. Métodos dos Elementos Finitos

O Método dos Elementos Finitos (MEF) é uma técnica de análise numérica para obtenção

de soluções aproximadas em problemas de valor de contorno, sendo aplicados uma grande

variedade de problemas de engenharia. Uma das principais publicações sobre MEF ilustra

sua aplicação na análise de tensões em estruturas aeronáuticas, mas sua utilização esten-

deu-se rapidamente para outros campos da engenharia (Huebner, 1995).

O MEF consiste em dividir (discretizar) o domínio de soluções do contínuo em uma quan-

tidade finita de subdomínios simples, denominados elementos finitos. Esses elementos

são representados matematicamente quase sempre por expansões polinomiais, denomina-

das funções de forma do elemento. O conjunto desses elementos, que podem assumir as

mais variadas formas, é denominado malha. A solução global é obtida pelo somatório das

soluções locais de cada elemento.

Um elemento finito pode ser unidimensional, bidimensional ou tridimensional. A Figura

8.1 ilustra alguns tipos de elementos finitos.

Elemento unidimensional linear de dois nós.

Elemento triangular de três nós.

Elemento quadrilátero de quatro nós.

Elementos quadrilátero de oito nós.


76

Elemento tridimensional de oito nós

Figura 8.1 – Tipos de elementos finitos

A fundamentação do método dos elementos finitos conta com três etapas distintas: sua

concepção matemática, aplicações e análise física e sua utilização na engenharia.

A formulação do MEF requer a existência de uma equação integral, de modo que seja pos-

sível substituir a integral sobre um domínio complexo (de volume V) por um somatório de

integrais de subdomínios de geometria simples (de volume Vi).

(8.1)

Onde se pressupõe que:

(8.2)

Se for possível calcular todos as integrais dos subdomínios Vi, basta efetuar o somatório

correspondente ao segundo membro de (1) para se obter a integral estendida a todo o do-

mínio. Cada subdomínio (Vi) corresponde a um elemento finito de geometria simples (e.g.,

segmento de reta, triângulo, quadrilátero, tetraedro, paralelepípedo).

Neste módulo apresentaremos problemas de engenharia em uma dimensão, nos quais o

método dos elementos finitos pode ser aplicado. Como nesses casos só existe uma variável

independente, as equações que regem estes problemas são equações diferenciais ordinárias.

Uma vez conhecida a equação diferencial que rege o problema e discretizado o domínio

(ou geometria) em elementos finitos, nós devemos encontrar funções que se aproximem da

solução verdadeira (geralmente desconhecida) dentro de cada elemento finito. Essas fun-

ções de aproximação são também conhecidas como funções de forma.

i

n

i 1v v

fdV fdV

n

i

i 1

V V

Profª Celme Torres

77

Para exemplificar, vamos considerar um problema unidimensional no qual h(x) é a variável

dependente. Então as funções de forma são usadas para interpolar o valor da variá-

vel h entre os nós dos elementos, da seguinte maneira:

(8.3)

Onde o superscrito e significa que a função só é válida dentro de um elemento e não

em todo o domínio. são os valores que a variável h assume nos nós do elemento e n é o

número de nós de cada elemento.

Escrevendo a equação acima em forma matricial:

(8.4)

Onde é a matriz linha que contém as funções de aproxi-

mação e é o vetor que contém o valor da variável h nos nós do elemento.

A maior dificuldade que os métodos aproximados estudados apresentam, como é o caso do

método dos elementos finitos, é a escolha da função de aproximação a ser empregada em

um problema complexo. Essas funções, em geral, têm que satisfazer exatamente as condi-

ções de contorno essenciais do problema e, além disso, representar adequadamente a geo-

metria e propriedades físicas do meio. Em termos computacionais, para poder-se usar esses

conceitos de forma sistemática, empregam-se funções de aproximação locais ao invés das

globais.

Desta maneira, discretiza-se (divide-se) inicialmente o domínio contínuo em um certo nú-

mero de elementos e aplica-se uma aproximação local da variável do problema em cada

um desses elementos. Fazendo-se, posteriormente, a junção de todos os elementos, obtém-

se um sistema de equações relativo ao modelo discreto, com um número finito de incógni-

tas.

iN (x)

ne e

i i

i 1

h (x) N (x)h

eh (x)

eih

e eh (x) [N]h

1 2 n[N] [N (x) N (x) ...N (x)]

1

2e

3

h

hh

h


78

8.1.1. Sistemas discretos e funções de base

8.1.2. Funções de base

As funções de base representam, em princípio, funções que são definidas e possuem valo-

res diferentes de zero em todo o domínio. A formulação clássica não fornece um método

sistemático para definição das funções de base, existindo para cada caso uma infinidade de

funções possíveis. Levando-se em conta que a precisão e qualidade dos resultados depende

fortemente da escolha das funções de base, isto se constitui uma séria desvantagem. A situ-

ação é ainda mais complicada quando se considera domínios em 2 e 3 dimensões com ge-

ometrias complexas. Uma escolha inadequada das funções de base pode levar a um sistema

de equações mal-condicionado e de difícil solução ou até mesmo impossível de resolver

numericamente. Por estes motivos, o Método de Galerkin possui uso restrito na solução de

casos práticos.

O MEF pode ser considerado uma modificação do Método de Galerkin onde as funções de

base são definidas de forma sistemática e levando a sistemas de equações numericamente

estáveis e fáceis de resolver, evitando desta forma as dificuldades mencionadas. No MEF

são utilizadas funções de base que só tem valores diferente de zero em uma pequena parte

do domínio. Conforme será visto, a escolha deste tipo de funções simplifica consideravel-

mente a obtenção da solução aproximada do problema. Devido à sua forma particular, elas

são também conhecidas como funções piramidais. O uso deste tipo particular de funções

constitui-se, desta forma, numa das características principais do MEF, sendo que a sua efi-

cácia está ligada à esta forma particular das funções de base. Exemplos de funções de base

piramidais em uma dimensão são mostrados na Figura 7.2.

Profª Celme Torres

79

Figura 8.2 –Subdivisão do domínio e funções de base: (a) subdivisão do domínio, (b) funções de base, (c) derivada das funções de base e (d) forma da função de aproxima-

ção.

Com o objetivo de definir conjunto de funções piramidais, como as mostradas na Figura

8.2, é necessário primeiramente dividir o intervalo [0,1] em um conjunto de n subinterva-

los, os quais por simplicidade serão todos do mesmo tamanho. Este conjunto de intervalos

será denominado genericamente de malha. Considerando, portanto, intervalos de igual ta-

manho resulta para a dimensão de cada um dos intervalos a relação, Figura 8.2(a):

(8.5)

Adotando-se a nomenclatura comum no MEF, os subintervalos são denominados de ele-

mentos e os pontos comuns de cada subintervalo são denominados de nós. Assim, cada e-

lemento é limitado pelos seguintes nós:

1 0 1h

n n


80

(8.6)

De acordo com a Figura 8.2(b), as funções de base são definidas como segue:

(8.7)

As funções lineares definidas por (7.7) são chamadas de funções de base definidas por tre-

chos. A sua característica principal é o fato de que elas são diferentes de zero apenas em

um pequena parcela do domínio, no caso considerado o domínio é o intervalo [0,1]. Por

este motivo elas são por vezes denominadas de funções de base local. Esta característica

particular simplifica enormemente a solução do problema por métodos numéricos. As deri-

vadas das funções de base são mostradas na figura 1c e dadas matematicamente por:

(8.8)

Com as definições anteriores, a solução aproximada será composta de segmentos de

retas, ou seja ela será to tipo linear por trechos, conforme mostra a figura 7.2(d). Observa-

se ainda que função de aproximação é contínua, não havendo saltos nos pontos de transi-

ção entre um elemento e outro. Esta característica é em geral desejada na solução. Haverá

entretanto uma descontinuidade na derivada da função, conforme mostra a figura 7.2(e).

Isto não causa, todavia, problemas. Caso seja necessário continuidade, também na derivada

de primeira ordem,da função de aproximação, poderão ser escolhidas funções de base de

ordem maior, por exemplo um polinômio de segunda ordem. O procedimento delineado

kx k h k 0,1,2,3,...n

k 1

k 1k 1 k

j

k 1k k 1

k 1

0 0 x x

x x x x x

h(x)

x x x x x

h

0 x x 1

k 1

k 1 kj

k k 1

k 1

0 0 x x

1 x x x

d (x) h

1dx x x x

h

0 x x 1

nu (x)

Profª Celme Torres

81

não se altera ao ser alterada a função de base. As matrizes do sistema no entanto possuirão

um número maior de elementos diferentes de zero.

8.1.3. Sistemas discretos

8.1.3.1. Mola elástica

Uma mola elástica é um elemento linear discreto onde a força-deslocamento é expresso

pela relação:

(8.9)

(8.10)

(8.11)

onde F é a força (N), é o deslocamento (m) e k é a constante elástica da mola (N/m).

A relação entre as forças ) e os deslocamentos de um elemento elástico tí-

pico (Figura 7.3) pode ser considerado um sistema discreto.

Figura 7.3 –Sistema linear elástico discreto

A força aplicada no nó 1é igual a constante de deslocamento multiplicada pelo deslo-

camento relativo do nó 1 com respeito ao nó 2, :

(8.12)

Igualmente, a força aplicada no nó 2 é igual a:

F ma

2 2

2 2

dv d x d xa m k

dt dt dt

F k

e e1 2(F ,F e e

1 2( , )

e1F

e e1 2

e e e e e1 e 1 2 e 1 e 2F k k k


82

(8.13)

Observe que a força de equilíbrio, , é automaticamente satisfeita. A equação

acima pode ser escrita em forma matricial, como:

(8.14)

A equação 7.12 é aplicada para qualquer elemento elástico onde a relação força-

deslocamento é linear.

8.1.3.2. Torção em barra circular

Outro problema que pode ser considerado como um sistema discreto é a torção de uma bar-

ra circular. Os vetores momento, T e T' designam-se momentos de torção ou momentos

torsores (torques).

Figura 7.4 – Torção em uma barra cilíndrica

Equação de equilíbrio de momentos em torno do eixo, é dada por :

(8.15)

A distribuição de tensões tangenciais é indeterminada, é por isso necessário analisar as de-

formações. A distribuição de tensões tangenciais não pode considerar-se uniforme no eixo,

ao contrário do que acontecia na barra com carregamento axial.

e e e e e2 e 2 1 e 1 e 2F k k k

e e1 2F F 0

e e1 1

e e e2 2

F1 1k

1 1 F

T dF dA

e e1 1T , e e

2 2T ,

he

Profª Celme Torres

83

No caso de mecânica dos sólidos, o ângulo de torção da seção circular de um corpo ci-

líndrico é descrita pelo torque T em torno do eixo, sendo expresso na forma:

(8.16)

Onde J é o momento, L é o comprimento e G é o momento de torção. A equação 6.14 pode

ser usada para descrever as relações entre os torques e os ângulos de um

elemento cilíndrico de comprimento . Escrevendo na forma matricial, temos:

(8.17)

8.1.3.3. Fluxo de um fluido em um tubo

Outro exemplo é dado para o caso de um fluido incompressível através de tubo circular. A

velocidade sendo laminar o fluxo do fluido viscoso através de um tubo circular é dada por:

(8.18)

Onde é o gradiente de pressão, d é o diâmetro do tubo e é a viscosidade do flui-

do.

GJT

L

e e1 2T ,T e e

1 2,

eh

e ee e1 1

e ee 2 2

T1 1G J

1 1h T

2

x

1 dP 2rv 1

4 dx d

dP / dx

ee 4

e

128 hR

d


84

Figura 8.5 – Fluxo viscoso através de um tubo cilíndrico

O volume da taxa de fluxo, Q, á obtida através da integração de em relação a área toda

seção transversal do tubo. Assim a relação entre Q e o gradiente de pressão é dada

pela equação:

(8.19)

O sinal negativo indica que o fluxo ocorre na direção do gradiente de pressão negativo.

A equação 7.19 pode ser aplicada para relação entre o valor nodal da taxa de volume de

fluxo, , a pressão, , em um tubo de comprimento e diâmetro . A

taxa de volume que entra no nó 1 é dada por:

(8.20)

Igualmente a taxa de volume que entra no 2 é dada por:

(8.21)

Logo, na forma matricial, temos:

A constante é chamada de resistência do tubo.

he1 2

xv

dP / dx

4d dPQ

128 dx

e e1 2Q ,Q e e

1 2P ,P eh ed

4

e e ee1 2 1

e

dQ P P

128 h

4

e e ee2 1 2

e

dQ P P

128 h

e e41 1

e ee 2 2

P Q1 1d

1 1128 h P Q

ee 4

e

128 hR

d

e1P e

2P

e1Q e

2Q

Profª Celme Torres

85

Exemplo 8.1

Seja o problema descrito pela equação e condições de contorno abaixo:

Condições de contorno:

Resolver a equação diferencial pelo método dos elementos finitos.

O domínio será discretizado em 5 elementos de mesma dimensão ( ), conectados

por meio de nós onde o valor da incógnita u será calculada.

Onde no elemento 1, considerou-se implicitamente que , condição de contorno do

problema.

A equação básica pelo Método de Galerkin para a equação diferencial, pode ser escrita

na seguinte forma discreta:

(ii)

O passo seguinte é a substituição das variações de u e , em cada elemento, no lado es-

querdo da equação (ii).

Tomando o elemento 2, temos:

2

2

d uu x 0 (1 x 0)

dx

u 0 em x 0

du1 em x 1

dx

h 0,2

u 0

i

i

x 15

x 1i 1 x

du d u(u x) u dx u

dx dx

u


86

Efetuando as integrações acima, chega-se a:

Repetindo-se as operações acima para todos os elementos e agrupando termo a termo,

chega-se à equação:

onde:

Como as variações são arbitrárias, a identidade acima só é válida se os termos que

multiplicam cada , dos dois lados da equação, forem idênticos simultaneamente. In-

troduzindo o valor de , temos o sistema global, dado por:

A solução do sistema é dada por:

Comparando graficamente as soluções:

2h

3 2 3 2 3 2 3 22

h

1 x x x xu u u u 1 u 2 u x 1 u 2 u dx

h h h hh

2

3 2 3 2 3 3 2 2 2 3 3 2 3 2

1 h h hu u u u u u u u u u u u 5 u 4 u

h 3 6 6

2 3 2

2 3 4 3

3 4 5 4

4 5 6 5

5 6 6 6

au bu c u

bu au bu 2c u

bu au bu 3c u

bu au bu 4c u

a 7u u c u u

2 3

2

2

3

a 4h 12

b h 6

c 6h

iu

iu

h 0,2

2

3

4

5

6

u 0,04811,84 6,04 0 0 0

u 0,0966,04 11,84 6,04 0 0

u 0,1440 6,04 11,84 6,04 0

0,1920 0 6,04 11,84 6,04 u

1,3120 0 0 6,04 11,84 u

u(0,2) 0,5335

u(0,4) 1,0379

u(0,6) 1,4851

u(0,8) 1,8494

u(1,0) 2,1085

Profª Celme Torres

87


88

Note que a matriz do sistema é simétrica e em banda, o que permite a utilização, em com-

putador, de algoritmos eficientes de armazenamento e solução.

Em termos computacionais, para uma melhor sistematização do programa, deve-se consi-

derar um elemento genérico (ao invés de se definirem funções de interpolação específicas

para cada elemento) e introduzir as condições de contorno do problema e fazer a expansão

das matrizes dos elementos para obtenção da matriz global. Conforme apresentado no Mé-

todo de Rayleigh-Ritz.

Exemplo 8.2

A parede apresentada na figura consiste de três materiais. De um lado da parede a tempe-

ratura é de 200 ºC e do outro lado a temperatura é de 50 ºC, com um coeficiente de con-

vecção de ( ). A equação que governa o problema é dada por:

Determine a temperatura na parede.

Condições de contorno:

100w 2m K

2

2

d TkA 0

dx 0 x L

h1h2

h3

1

h1

2 43

Material 1 Material 2 Material 3

T = 200 C0

o T = 50 Cf

o

1

2

3

1

2

3

2

k 70W /(m K)

k 40W /(m K)

k 20W /(m K)

h 2 cm

h 2,5 cm

h 4 cm

10W /(m K)

Profª Celme Torres

89

Onde A é a área da seção transversal (que pode ser considerada igual a 1) e k é a condu-

tividade.

Solução

Para uma malha com 3 elementos lineares, considerando a matriz de rigidez dos elemen-

tos, temos:

A matriz global do sistema é dada por:

As condições de contorno e as condições de equilíbrio, são dadas por:

Omitindo a primeira linha e a primeira coluna, temos:

o

f

x L

T(0) T

dTkA A T T 0

dx

1 1

1

2 2

2

3 3

3

1 1 1 1 3500 3500k A 70.1[K]

1 1 1 1 3500 3500h 0,02

1 1 1 1 1600 1600k A 40.1[K]

1 1 1 1 1600 1600h 0,025

1 1 1 1 500 500k A 20.1[K]

1 1 1 1 500 500h 0,04

[K] U Q

111

1 22 2 1

2 33 2 1

34

2

QU3500 3500 0 0

U Q Q3500 3500 1600 1600 0

U0 1600 1600 500 500 Q Q

0 0 500 500 U Q

1

2 11 2

2 32 1

32 4 f

U 200

Q Q 0

Q Q 0

Q A U T


90

A solução é dada por:

ºC

ºC

ºC

ºC

A quantidade de calor por unidade de área é dada por:

W/m2 na parte externa (lado esquerdo)

W/m2 na parte interna (lado direito)

Exemplo 8.3

A estrutura representada na figura é unidimensional, tem quatro nós (1 a 4) e quatro bar-

ras (A a D). Cada barra tem as suas características, nomeadamente, o módulo de Young

(E), a área da secção transversal (A) e o comprimento (L). Em cada nó existe um único

grau de liberdade. Em correspondência com os quatro graus de liberdade existem quatro

deslocamentos nodais (a) e quatro forças nodais equivalentes à ação exterior (F). Cada

barra tem dois graus de liberdade (um em cada extremidade).

2

3

44

U5100 1600 0 3500 200

1600 2100 500 U 0

0 500 510 10U 500U

52

3

4

U5100 1600 0 70 10

1600 2100 500 U 0

0 500 510 500U

1U 200

2U 199,58

3U 198,67

4U 195,76

11 1 2Q 3500U 3500U 1457,6

32 4Q 10 1 U 50 1457,6

2,0 m 2,5 m 2,5 m

Dados:

E = 2500 Pa;

A = 20cm x 40cm

Forças: F1: 4 N

F2: 3 N

F3: 5 N

F4: 2 N

Profª Celme Torres

91

Usando o MEF encontre os deslocamentos nodais (a1, a2, a3 e a4) das barras da estrutura

sabendo que a Equação diferencial que rege o problema é dada por:

d du

EA 0dx dx

Solução

# Matriz de rigidez das barras:

Barra A: A 1 1 1 1 100 100EA 2500 0,08[K]

1 1 1 1 100 100L 2,0

Barra B: B 1 1 1 1 80 80EA 2500 0,08[K]

1 1 1 1 80 80L 2,5

Barra C: C 1 1 1 1 80 80EA 2500 0,08[K]

1 1 1 1 80 80L 2,5

Barra D: D 1 1 1 1 40 40EA 2500 0,08[K]

1 1 1 1 40 40L 5,0

# Montagem da matriz de rigidez global:

RELEMBRANDO: Mecânica dos Sólidos

O módulo de Young ou módulo de elasticidade é um parâmetro mecânico que proporciona uma

medida da rigidez de um material sólido. Obtém-se da razão entre a tensão (ou pressão) exerci-

da e a deformacão unitária sofrida pelo material. Isto é,

Tensão F/ A FLE

Deformação a / L Ax

onde (em unidades do SI):

E é o módulo de Young, medido em pascal (1 Pa ≡ 1 N/m²);

F é a força medida em Newton;

A é a secção através da qual é exercida a tensão, e mede-se em metros quadrados;

a é a deformação, medido em metros;

L é o comprimento natural medido em metros.


92

A A11 12

A A B D B D21 22 11 11 12 12

B B C C21 22 11 12

D C C D21 21 22 22

K K 0 0

K K K K K K[K]

0 K K K K

0 K K K K

100 100 0 0 100 100 0 0

100 100 80 40 80 40 100 220 80 40[K]

0 80 80 80 80 0 80 160 80

0 40 80 80 40 0 40 80 120

A formulação do MEF, uma vez que a relação envolvendo todos os graus de liberdade da

estrutura, é dada por:

[K]a F

# Montagem do vetor de força F

O vetor de força é a soma algébrica de todas as forças que atuam no nós da estrutu-

ra,sendo portanto:

1

2

3

4

F 4,0

F 3,0F

F 5,0

F 2,0

# Montagem do vetor dos deslocamentos em todos os graus de liberdade da estrutura a

1

2

3

4

a

aa

a

a

# Solução do sistema

1

2

3

4

a100 100 0 0 4,0

a100 220 80 40 3,0

a0 80 160 80 5,0

a0 40 80 120 2,0

1

2

3

4

a 0,66

a 0,70

a 0,64

a 0,64

Profª Celme Torres

93

Aplicando o Maple temos:

A matriz de rigidez global do sistema é dada por:

O vetor de carga nodal é:

Os deslocamentos são dados pelo vetor a, onde:


94

9. Métodos dos Elementos Finitos Bidimensional

A análise bidimensional de elementos finitos envolve alguns conceitos básicos, assim co-

mo foi descrito para problemas unidimensionais. A análise do problema é algumas vezes

complicada pelo fato dos problemas bidimensionais serem descritos por equações diferen-

ciais parciais sobre uma região geometricamente complexa. O contorno () de um domínio

() bidimensional é em geral uma curva. Assim, em um problema bidimensional nos não

só procuramos uma solução aproximada de um dado problema em um domínio, como tam-

bém aproximamos o domínio de uma malha de elementos finitos adequada. Conseqüente-

mente, iremos ter erros de aproximação devido a aproximação da solução, como também,

erros de discretização devido a aproximação do domínio por elementos finitos.

A malha de elementos finitos (discretização) consiste de simples elementos bidimensio-

nais, tal como, triângulos, retângulos ou quadriláteros, que permite uma única derivação

das funções de interpolação.

No desenvolvimento de um modelo de elementos finitos bidimensional, é necessário seguir

alguns passos:

Discretização do domínio em um conjunto de elementos finitos;

Formulação ponderada da equação diferencial;

Encontrar as funções de interpolações dos elementos finitos;

Desenvolvimento do modelo usado às funções de ponderação;

Agrupamento dos elementos finitos para obtenção de um sistema global de equa-

ções algébricas;

Imposição das condições de contorno;

Solução das equações;

Pós-processamento da solução graficamente.

Profª Celme Torres

95

9.1.1. Modelos de elementos finitos

A formulação que descreve o modelo de elementos finitos bidimensional é dada por:

n

e e e eij j i i

j 1

K u f Q

(i 1,2,3,...,n) (8.1)

onde,

e

e e e ee ej j j je e ei i

ij 11 12 21 22 00 i j

N N N NN NK a a a a a N N dxdy

x x y y x y

(8.2)

e

e ei if fN dxdy

(8.3)

e

e ei n iQ q N ds

(8.4)

Em notação matricial, a equação (8.4) fica na forma:

e e e eK u f Q

(8.5)

Note que e eij jiK K (i.e., [K

e] é uma matriz simétrica de ordem n n , somente quando

12 21a a

9.1.2. Funções de forma

Se u(x, y) é a variável dependente de um determinado problema genérico bidimensional,

então as funções de forma iN (x, y) são usadas para interpolar o valor da variável u entre

os nós dos elementos da seguinte forma:

ne e

i i

i 1

u (x, y) N (x, y)u

(8.6)


96

Onde o supercrito e significa que a função eu (x, y) só é válida dentro do elemento e não

em todo domínio; eiu são os valores que a variável u assume nos nós do elemento e n é o

número de nós de cada elemento.

A equação (8.6) pode ser escrita em forma matricial, como segue:

e eu (x, y) N u (8.7)

onde, 1 2 n[N] N (x, y) N (x, y) ... N (x, y) é a matriz linha que contém as funções de

forma, ou funções de aproximação; e eu é o vetor que contém os valores da variável u

nos nós do elemento, sendo dado por:

1

2e

n

u

uu

u

(8.8)

9.1.3. Elemento triangular linear

Considerando o elemento triangular linear apresentado na figura 8.1.

Figura 8.1 – Elemento triangular linear

h

y

x

ui

uk

uj

i

(xi, yi) k

(xk, yk) j

(xj, yj)

Profª Celme Torres

97

O polinômio de interpolação é dado por:

1 2 3u x y (8.9)

O polinômio de interpolação deve obedecer as seguintes condições nodais:

i i i

j j i

k k k

u u x x ; y y

u u x x ; y y

u u x x ; y y

(8.10)

O qual pode ser resolvido para 1 2 3, e :

j k k j i k i i k j i j j i k

1

x y x y u x y x y u x y x y u

2A

(8.12)

j k i k i j i j k

2

y y u y y u y y u

2A

(8.13)

k j i i k j j i k

1

x x u x x u x x u

2A

(8.14)

Com,

i i

j j

k k

1 x y

2A 1 x y

1 x y

(8.15)

Onde A é a área do elemento triangular.

Substituindo 1 2 3, e no polinômio e rearranjando os termos:

i i j j k kh N u N u N u (8.16)

Onde, Ni, Nj e Nk são as funções de forma, dadas por:

i i ii

j j j

j

k k kk

a b x c yN

2A

a b x c yN

2A

a b x c yN

2A

(8.17)

Os coeficientes a, b e c foram definidos nas equações (8.12), (8.13) e (8.14), sendo:


98

i j k k j

i j k

i k j

a x y x y

b y y

c x x

j k i i k

j k i

j i k

a x y x y

b y y

c x x

k i j j i

i i j

i j i

a x y x y

b y y

c x x

(8.18)

Para derivação e exemplificação do método dos elementos finitos aplicado a problemas

bidimensionais no estado permanente, será usada a seguinte equação diferencial:

x y

u ug g Q 0

x x y y

(8.19)

onde, u é a variável desconhecida, gx e gy são coeficientes que regem o problema descrito

pela ED(constantes) nas direções x e y e Q é um termo que representa forças ou cargas a-

plicadas.

De acordo com o Método dos Resíduos Ponderados de Galerkin, para problemas bidimen-

sionais, temos:

e eT T T

x yA A A

u u[N] R(x, y)dA [N] g g dA [N] QdA 0

x x y y

(8.20)

9.1.4. Teorema de Green-Gaus

Segundo o Teorema de Green-Gaus em 2D, dado um vetor q que é o produto de um escalar

por um vetor p, e uma área A, cujo contorno delimitador é C, o qual tem um vetor unitá-

rio normal n, como mostrado na Figura 8.2.

Figura 8.2 – Teorema de Green-Gaus.

Área A

Contorno C

j

ix

y

dA

dC

q

n

Profª Celme Torres

99

A C

q dA q n dC (8.21)

onde,

i jx y

e x yn n i n j (8.22)

sendo, i e j os vetores unitários nas direções x, y.

Mas, 𝑞 = 𝛽𝑝, portanto,

A C

p dA p ndC (8.23)

Sabendo que:

( p) p p (8.24)

Logo,

A A C

pdA pdA p ndC (7.25)

ou

A C A

pdA p ndC pdA (8.26)

A equação (8.26) define o Teorema de Green- Gauss.

Voltando a integral (8.20), o primeiro termo a direita, dado por:

e e

Tx y

A

u u[N] g g dA

x x y y

(8.27)

Pode ser escrito na forma:

e

xT

eA

y

ug

x[N] dA

ug

y

(8.28)

Aplicando o Teorema de Green-Gaus


100

e

xT

e

y

ug

x[N] p

ug

y

(8.29)

Substituindo no teorema de Green-Gaus, temos:

e e

x xe eT T T

x x y ye e

A C Ay y

(a)

(b)

u ug g

u ux x[N] dA [N] g n g n dC [N] dA

x yu ug g

y y

(8.30)

Analisando o primeiro termo da equação (8.30), identifica-se que o termo (a) está relacio-

nado com a forças ou cargas aplicadas no contorno do domínio, assim o termo (a) da equa-

ção (8.30) fica:

T

C[N] qdC (8.31)

Sendo, q o módulo ou valor das forças ou cargas aplicadas no contorno do elemento. Por

convenção, é adotado o valor positivo para as forças ou cargas aplicadas entrando no do-

mínio.

A integral (8.31) é válida para qualquer elemento bidimensional e pode ser determinada

desde que as funções de forma [N] estejam determinadas. Esta integral representa a con-

tribuição das condições de contorno naturais do problema.

9.1.5. Imposição das condições de contorno

Existem duas condições de contorno, as chamadas condições de contorno essenciais ou

geométricas e as condições de contorno naturais.

As condições de contorno naturais se referem à derivada da função que está sendo anali-

sada e as condições de contorno geométricas ou essenciais, são valores especificados da

função ao longo do contorno ou domínio. As condições de contorno geométricas são in-

corporadas diretamente no vetor das incógnitas u, como será mostrado a seguir.

Profª Celme Torres

101

9.1.6. Determinação da matriz de rigidez

Voltando a equação (8.31), temos:

eT

qC[N] qdC f (8.32)

que é chamado de vetor de carga do elemento.

Analisando o termo (b) da equação (8.30)

e

xT

eA

y

ug

x[N] dA

ug

y

, temos:

e e

xx

e ey

y

u ug

g 0x x

0 gu ug

y y

(8.33)

Como, eh (x, y) [N(x, y)]

e e1 2 nh N N ... N h (8.34)

e1 2 n

e e

e1 2 n

N N Nu...

x x x xh [N]h

N N Nu ...y y yy

(8.35)

Chamando, x

y

g 0[D]

0 g

, temos:

e

xT T e

eA A

y

ug

x[N] dA [N] [D] [N]dA f

ug

y

(8.36)

A integral T e

A[N] [D] [N]dA f

é chamada Matriz de Rigidez, portanto:


102

e

xT e e

n n n 1eA

y

ug

x[N] dA [K ] f

ug

y

9.1.7. Determinação do vetor de carga nodal

Analisando o segundo termo do lado direito da equação (8.20), temos:

T

A[N] QdA (8.38)

Essa integral representa a contribuição das forças e cargas no Vetor de Carga Nodal, lo-

go:

T eQ

A[N] QdA f (8.39)

Portanto, a equação matricial (8.5) para solução de problemas bidimensionais através do

método de elementos finitos é composta por uma matriz de rigidez [K] um vetor de incóg-

nitas e o vetor de carga nodal que incorpora, em sua derivação, as condições de contorno

do problema.

e e e eq QK u f f 0

(8.40)

9.1.8. Formulação de elementos finitos

Fazendo,

e e eq Qf f f (8.41)

Temos,

e e eK u f

(8.42)

O sistema apresentado na equação (8.42) é um sistema de equações lineares no qual ue é

o vetor das incógnitas, resultante da formulação do problemas transcrito por uma equação

diferencial parcial bidimensional.

Profª Celme Torres

103

Exemplo 9.1

Elementos finitos aplicados à análise bidimensional de problemas de fluxo hídrico

subterrâneo.

A equação diferencial que rege o problemabidimensional de fluxo hídricos subterrâneo em

estado permanente é dada por:

x y

h hk k 0

x x y y

Considere um aqüífero no qual o potencial hidráulico é conhecido no contorno do domí-

nio, sendo dado por h 10 m no lado esquerdo, h 5 m no lado direito e h 0,0 m no

contorno impermeável. Sabendo que a condutividade hidráulica do aqüífero é

5x yk k 25 10 m/s, determine o potencial hidráulico no interior do aqüífero.

Solução

O domínio foi dividido em 4 elementos e 5 nós.

Área dos elementos: A 1,0 m2

# Elemento 1 i = 1; j = 2; k = 5

11

11

x 0

y 0

12

12

x 0

y 0

15

15

x 1

y 1

1

1

b 0 1 1

c 1 2 1

2

2

b 1 0 1

c 0 1 1

5

5

b 0 0 0

c 2 0 2

Matriz de rigidez do elemento 1:

//////////////////////////////////////////

/////////////////////////////////////////

2,0 m h = 10 m h = 5 m

h = 0 m

h = 0 m

impermeável

2,0 m

y

x

5

3

2

4

1

1

2

4

3


104

1 5

12,5 0 12,5

[K ] 0 12,5 12,5 10

12,5 12,5 25

# Elemento 2 i = 1; j = 5; k = 4

1

1

b 1 2 1

c 0 1 1

5

5

b 2 0 2

c 0 0 0

4

4

b 0 1 1

c 1 0 1

Matriz de rigidez do elemento 2:

2 5

12,5 12,5 0

[K ] 12,5 25 12,5 10

0 12,5 12,5

# Elemento 3 i = 2; j = 5; k = 3

3 2 5

12,5 12,5 0

[K ] [K ] 12,5 25 12,5 10

0 12,5 12,5

# Elemento 4 i = 4; j = 5; k = 3

4 1 5

12,5 0 12,5

[K ] [K ] 0 12,5 12,5 10

12,5 12,5 25

# Montagem da matriz global do sistema

5

25 0 0 0 25

0 25 0 0 25

[K] 100 0 25 0 25

0 0 0 25 25

25 25 25 25 100

# Vetor de carga nodal e vetor das incógnitas

Profª Celme Torres

105

O fluxo especificado no contorno do domínio é q = 0, devido ao contorno impermeável. O

nó 5 é um nó interno e deve atender a equação da continuidade, sendo portanto q5 = 0,

logo o fluxo que entra, necessariamente, é igual ao fluxo que sai do sistema.

Os valores das incógnitas (h), potencial hidráulico nos nós 1, 2, 3 e 4 são conhecidos, logo

o sistema global é descrito por:

1

2

53

4

5 5

q1025 0 0 0 25

q00 25 0 0 25

[K] 10 0 q0 0 25 0 25

100 0 0 25 25 q

h25 25 25 25 100 q

Tomando a última linha do sistema encontramos o valor da incógnita h5 5h 5,0 m

Graficamente, temos a variação do potencial hidráulico, ou a superfície piezométrica, do

exemplo em estudo:

O mapa dos vetores de fluxo pode ser observado na figura abaixo, indicando que o movi-

mento da água ocorre do maior para o menor potencial hidráulico. No contorno superior

e inferior, como não existe fluxo (h = 0) os vetores são adjacentes ao contorno.

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

4.4

4.8

5.2

5.6

6

6.4

6.8

7.2

7.6

8

8.4

8.8

9.2

9.6

10

Carga hidráulica h (m)


106

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Profª Celme Torres

107

Referências Bibliográficas

AKIM, J.E. Finite Element for Analysis and Design. Editora Academic Press. 1994.

REDDY, J.N. An Introduction to the Finite Element Analysis. Editora Mc Graw Hill.

1984.

COOK, R.D. Concepts and Aplications os the Finite Element Analysis. Editora Jonh Wiley

& Song. 1989.

ZIENKIEWICZ, O. C.; TAYLOR, R. L.; ZHU, J. Z. The Finite Element Method: Its Basis

And Fundamentals.

ASSAN, A. E. Método dos Elementos Finitos: Primeiros Passos. 2ª Edição. Editora Uni-

camp. 2003.

Engineering

Apostila metodos numericos_formatada_vf