Introducción al Método de los Elementos Finitos• Para llevar el problema de Stokes a la forma...

Alberto Cardona, Víctor Fachinotti

Cimec (UNL/Conicet), Santa Fe, Argentina

25/Oct/16

Introducción al Método

de los Elementos Finitos

Parte 5

Algunas aplicaciones del MEF a problemas elípticos

v2

v12

1

5

S

4

Introducción al Método de los Elementos Finitos 2

Algunas aplicaciones del MEF a problemas elípticos

Problema de Elasticidad

• Consideremos un cuerpo elástico, isótropo y homogéneo B, que ocupa el dominio acotado W3, con frontera

G=GtGu tal que GtGu= y área Gu >0.

• El cuerpo B está sometido a una fuerza volumétrica f,

y a una fuerza superficial t aplicada sobre Gt.

• Se supone B fijo a lo largo de Gu.

• Se busca determinar

– desplazamientos u.

– deformaciones e, dependientes de los desplazamientos de acuerdo a la

cinemática de la deformación.

– tensiones s, dependientes de las deformaciones de acuerdo a la ley

constitutiva del material.

WW Gu

Gt


Problema de Elasticidad

• Ecuaciones de clausura

– Cinemáticas: asumiendo pequeñas deformaciones:

– Constitutivas: asumiendo comportamiento elástico lineal (ley de Hooke):

• Nota: en adelante, usaremos las siguientes convenciones de notación:

- derivada parcial: - sumatoria:

div en Ecuación de equilibrio

sobre CB Dirichlet (despl. impuesto)

sobre CB Neumann (tracción impuesta)

u

t

W

G G

σ f 0

u 0

σn t

WW Gu

Gt

n

1

2

jiij

j i

uu

x xe

: módulo de Elasticidadcon , ,

: coeficiente de Poisson1 (1 )(1 2 )

EE E

,i

i j

j

uu

x

3

1

ij j ij j

j

n ns s

3

1

, , , ctes. de Lamé. ij kk ij ij

k

s e e


Forma variacional del problema de Elasticidad

• Dado

se llega a

haciendo

3

1

Hallar V / a( , ) L( ), V

( )

con V= : H ( ) y en u

V

W G

u u v v v

v v v 0

, 0 en ,( , 1,2,3) ij j if i jD s W

,

,

, ,

, ,

,

0

0

02

02

( ) 0

( ) ( )

t

t

t

ij j i i i

ij i j ij j i i i

ij i j ji j i

i i i i

i j j i

ij i i i i

ij ij i i i i

i i ij ij ij

v dx f v dx

v dx n v ds f v dx

v vdx t v ds f v dx

v vdx t v ds f v dx

dx t v ds f v dx

u

s

s s

s s

s

s e

e e

W W

W G W

W G W

W G W

W G W

v

u v 0t

i i i idx t v ds f v dxW G W

Cinemática

de pequeñas

deformaciones

T. de Green

Ley de Hooke

CB

Simetría de s

L( )va( , )u v


Forma variacional del problema de Elasticidad (cont.)

• Se puede demostrar que la forma lineal es continua,

i.e.,

• Se puede demostrar que la forma bilineal

– es simétrica, i.e.,

– continua, i.e.,

– V-elíptica, i.e., con

Demo.:

,a( , ) ( ) ( ) div div ( ) ( )k k ij ij ij ij iju dx dx e e e eW W

u v u v u v u v

1 1

+

H ( ) H ( )a( , ) , , V, .

W W u v u v u v

1 1

32 2

H ( ) H ( )1

.i

i

vW W

v1

2 +

H ( )a( , ) , V, ,

W v v v v

a( , ) a( , ), , V. u v v u u v

L( )t

i i i if v dx t v dsW G

v

VL( ) , V, . v v v

1

2

2

H ( )

a( , ) div ( ) ( )

( ) ( ) , .

ij ij

ij ij

dx dx

dx c c

e e

e e

W W

WW

v v v v v

v v v

Desigualdad de Korn


MEF aplicado al problema de Elasticidad

• Consideremos el problema de Elasticidad en W2.

• Sea Th={K} una malla de triángulos de W. Definimos el espacio de EF

• El MEF aplicado al problema de Elasticidad consiste en

• La solución uhVh satisface

2

h 1KV : V y P (K) , K Th

v v v

1 2H ( ) H ( )h ChW W

u u u

Hallar V / a( , ) L( ), Vh h h h u u v v v


Funciones de base para el triángulo lineal

Toda función puede representarse

Las funciones de base Ni = i ( coord de área del triángulo K) resultan:

con:

1

31

1 1 2 3 2

31 2 3 2

1 3

3

( , )( , ) 0 0 0

( , ) , ( , ) K.( , ) 0 0 0

( , )

x

y

x

i i xx i

yy y

i i xi

y

V

VN x y V

v x y N N N Vx y x y

v x y N N N VN x y V

V

V

v

( , )i i i iN x y a b x c y

2 3 3 21

K2

x y x ya

A

2 3

1

K2

y yb

A

3 2

1

K2

x xc

A

1

3

1

2

1

1

3

1

2

2

1

3

1

2

3

3 1 1 32

K2

x y x ya

A

3 1

2

K2

y yb

A

1 3

2

K2

x xc

A

1 2 2 13

K2

x y x ya

A

1 2

3

K2

y yb

A

2 1

3

K2

x xc

A

2

1P (K) v



Notar:

( ) 0 ( )

a( , ) ( ) ( ) ( ) 0 ( )

2 ( ) 2 ( )0 0

2

T

xx xx h

h jj ii h ij ij h yy yy h

xy xy h

dx dx

e e

e e e e e e

e eW W

v u

v u v u v u v u

v u

131 2

1

231 2

2

3 31 1 2 2 3

3

0 0 0

( )

( ) 0 0 0

2 ( )

x

xy

xx xy

yy y

xy x

yx

y

Vv NN N

Vx x x x

v VNN N

y y y y V

v N NN N N N Vv

y x y x y xy x V

e

e

e

v

ε v v

v

B V

1 2 3

1 2 3

1 1 2 2 3 3

0 0 0

0 0 0

b b b

c c c

c b c b c b

B

El vector de deformación :

Para el triángulo lineal, la matriz es constante :B

h hε u BU



Luego:

con

En consecuencia, la matriz de rigidez elemental :

Notar:

K Ka ( , ) T T T T

h h hK

dx A v u V B DBU V B DBU

Ka( , ) a ( , )h h

K

v u v u

K

T

K AA B DB

3

1

0

0

20 02

xx xx

ij kk ij ij yy yy

k

xy xy

s e

s e e s e

s e

σ Dε

0 1 / 1 01

0 / 1 1 01 1 2

1 20 0 0 02 2 1

E

D

(estado plano de deformación)



Problema de Stokes

• Consideremos las ecuaciones de Stokes para el flujo estacionario de un fluido Newtoniano incompresible encerrado en un dominio W3, sometido a una

fuerza volumétrica f :

• Definimos el espacio de funciones de prueba

• Luego, podemos llevar el problema de Stokes a la forma variacional

,

,

0 en , Balance de cant. de movto. : velocidad

2 ( ) en , Ley const. de fluido Newtoniano : tensión

0 en , Condición de incompresibilidad : presión

0 sobre , CB Dirichlet : viscos

ij j i

ij ij ij

i i

i

f

p

u p

u

s

s e

W

W

W

G

u

u σ

idad

, en , Balance de cant. de movto. p/fluido Newtonianoi i iu p f W

3

1

0V : H ( ) y div 0 en W W v v v

Hallar V / a( , ) L(( ), V)V u u v v v


Forma variacional del problema de Stokes

• Para llevar el problema de Stokes a la forma variacional hacemos

• Dado que 0, se demuestra (ídem problema de Poisson) que a(.,.) es simétrica,

contínua y V-elíptica.

• Se demuestra también (ídem problema de Poisson) que L(.) es continua.

• Nota: al adoptar un espacio de velocidades de divergencia nula, la formulación variacional no involucra la presión.

a( , )u v

,

,

,

0 00

i i i

i i i i i i

ii i i i i i i i i

i i i i

f u p

f v dx u v dx p v dx

uf v dx u v dx v ds pv dx pn v ds

n

f v dx u v dx

W W W

W W G W G

W W

L( )v


Elemento P5 - MEF aplicado al problema de Stokes

Consideremos el problema de Stokes en W2. Luego:

• Si W es simplemente conexo (i.e., no contiene agujeros), div v=0 en W si y solo si

o sea:

• Adoptamos luego un subespacio Wh de dimensión finita de (usamos por

ej. el elemento finito C1-continuo ya visto) y definimos

• Se formula el MEF remplazando V por VhV en la formulación variacional. La

solución uhVh satisface

21 1 2

1 2 0

1 2

V : ( , ) H ( ) y 0 en v v

v vx x

W W

v v

2 1

, rot para alguna función .

: función de corriente del campo de velocidades .

x x

v

v

2

0V rot , H ( ). Wv v

2

0H ( )W

V : rot , W .h h v v

1 5

4

H ( ) H ( )h ChW W

u u u


, ,i i i j i ju v d u v d x x

3

0

rot rot 0

u e

,rot ijk k jtet

3 ,i ij ju e

,2

,1

0

u

, 3 ,i j ik kju e

, , 3 , 3 , 33 3 3 , , , , ,3 ,3i j i j ik kj il lj kl k l kj lj lj lj j ju v e e

, , ,1 ,1 ,2 ,2i j i j j j j ju v d d x x


Elemento P5 - Desarrollo del MEF para Stokes

2 2

00 10 01 20 11 02

0 5

( , ) , , ( , ) K.i j T

ij

i j

x y c x y c c x c y c x c x y c y x y x y

p c

2 2

00 10 01 20 11 02

0 5

( , ) , 1,2,3.i j k

k k ij k k k k k k k k

i j

x y c x y c c x c y c x c x y c y k

0 0

0 0

2 2 3

00 10

; ; : longitud característica de elemento (ver adelante)

, : centroide del elemento

, 1T

x x y yx y L

L L

x y

x y x y x x y y x

c c

p

c 01 20 11 02 30 21 12 03 40 31 22 13 23 14 05

Tc c c c c c c c c c c c c c c

2 2 3 2 2 3 4 3 2 2 3 4 5 4 3 2 2 3 4 51

1,2,3.

k

k k k k k k k k k k k k k k k k k k k k k k k k k k k k k kx y x x y y x x y x y y x x y x y x y y x x y x y x y x y y

k

c

Desarrollaremos un EF triangular C1-continuo. El campo será aproximado por

polinomios de grado 5:

( , )x y

con:

Para el cálculo de los coeficientes c expresamos el valor del campo en los nodos

Obtenemos así tres ecuaciones con 21 incógnitas:


Calculando luego el valor de la derivada respecto de x en los nodos, logramos

tres ecuaciones más:

Repitiendo el proceso para la evaluación de todos los grados de libertad nodales,

obtenemos un conjunto de 18 ecuaciones con 21 incógnitas.

Este proceso puede hacerse usando un programa de manipulación simbólica,

como se indica en el siguiente slide.


( 1)

, 10 20 11 ,

0 5

1 2 1( , ) , 1,2,3.i j k

x k k ij k k k k x

i j

ix y c x y c c x c y k

L L L L

2 2 3 2 2 3 4 3 2 2 3 4

,0 1 0 2 0 3 2 0 4 3 2 0 5 4 3 2 0

1,2,3.

k

k k k k k k k k k k k k k k k k k k k k xx y x x y y x x y x y y x x y x y x y y L

k

c



clear

syms C x y x1 y1 x2 y2 x3 y3 L

C = 1;

for i=1:5

for j=0:i

C = [C (x/L)^(i-j)*(y/L)^(j)];

end

end

Cx = diff(C,x);

Cy = diff(C,y);

Cxx = diff(C,x,2);

Cxy = diff(Cx,y);

Cyy = diff(C,y,2);

% NODO 1

% phi(x1,y1) = phi_1

A(1,:) = subs(C ,{x,y},{x1,y1});

% dphi/dx (x1,y1) = dphix_1

A(2,:) = subs(Cx ,{x,y},{x1,y1});

% dphi/dy (x1,y1) = dphiy_1

A(3,:) = subs(Cy, {x,y},{x1,y1});

% d2phi/dx2 (x1,y1) = dphixx_1

A(4,:) = subs(Cxx,{x,y},{x1,y1});

% d2phi/dxdy (x1,y1) = dphixy_1

A(5,:) = subs(Cxy, {x,y},{x1,y1});

% d2phi/dy2 (x1,y1) = dphiyy_1

A(6,:) = subs(Cyy,{x,y},{x1,y1});

% NODO 2

A(7,:) = subs(C, {x,y},{x2,y2});

A(8,:) = subs(Cx, {x,y},{x2,y2});

A(9,:) = subs(Cy, {x,y},{x2,y2});

A(10,:) = subs(Cxx, {x,y},{x2,y2});

A(11,:) = subs(Cxy, {x,y},{x2,y2});

A(12,:) = subs(Cyy, {x,y},{x2,y2});

% NODO 3

A(13,:) = subs(C, {x,y},{x3,y3});

A(14,:) = subs(Cx, {x,y},{x3,y3});

A(15,:) = subs(Cy, {x,y},{x3,y3});

A(16,:) = subs(Cxx,{x,y},{x3,y3});

A(17,:) = subs(Cxy,{x,y},{x3,y3});

A(18,:) = subs(Cyy,{x,y},{x3,y3});

Programa para calculo simbólico coeficientes


2 2 3 2 2 3 4 3 2 2 3 4 5 4 3 2 2 3 4 5

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

1

1

0 1 0 2

x y x x y y x x y x y y x x y x y x y y x x y x y x y x y y

x

2 2 3 2 2 3 4 3 2 2 3 4

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

2

1 1 1 1 1

0 3 2 0 4 3 2 0 5 4 3 2 0

0 0 1 0 2 0 2

y x x y y x x y x y y x x y x y x y y

x y x x y

2 3 2 2 3 4 3 2 2 3 4

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

2 2

1 1 1 1 1 1

3 0 2 3 4 0 2 3 4 5

0 0 0 2 0 0 6 2 0 0 12 6 2

y x x y x y y x x y x y x y y

x y x x y y

3 2 2 3

1 1 1 1 1 1

2 2 3

1 1 1 1 1 1 1

0 0 20 12 6 2 0 0

0 0 0 0 1 0 0 2 2 0 0 3 4 3 0 0 4

x x y x y y

x y x x y y x

2 2 3

1 1 1 1 1

2 2 3 2 2 3

1 1 1 1 1 1 1 1 1 1 1 1

6 6 4 0

0 0 0 0 0 2 0 0 2 6 0 0 2 6 12 0 0 2 6 12 20

x y x y y

x y x x y y x x y x y y

2 2 3 2 2 3 4 3 2 2 3 4 5 4 3 2 2 3 4 5

2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2

2

1

0 1 0 2


x

2 2 3 2 2 3 4 3 2 2 3 4

2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1

2

2 2 2 2 2

0 3 2 0 4 3 2 0 5 4 3 2 0

0 0 1 0 2 0 2


x y x x y

2 3 2 2 3 4 3 2 2 3 4

2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

2 2

2 2 2 2 2 2

3 0 2 3 4 0 2 3 4 5

0 0 0 2 0 0 6 2 0 0 12 6 2


x y x x y y

3 2 2 3

2 2 2 2 2 2

2 2 3

2 2 2 2 2 2 2

0 0 20 12 6 2 0 0

0 0 0 0 1 0 0 2 2 0 0 3 4 3 0 0 4

x x y x y y

x y x x y y x

2 2 3

2 2 2 2 2

2 2 3 2 2 3

2 2 2 2 2 2 2 2 2 2 2 2

6 6 4 0

0 0 0 0 0 2 0 0 2 6 0 0 2 6 12 0 0 2 6 12 20

x y x y y


2 2 3 2 2 3 4 3 2 2 3 4 5 4 3 2 2 3 4 5

3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

3

1

0 1 0 2


x

2 2 3 2 2 3 4 3 2 2 3 4

3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

2

3 3 3 3 3

0 3 2 0 4 3 2 0 5 4 3 2 0

0 0 1 0 2 0 2


x y x x y

2 3 2 2 3 4 3 2 2 3 4

3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

2 2

3 3 3 3 3 3

3 0 2 3 4 0 2 3 4 5

0 0 0 2 0 0 6 2 0 0 12 6 2


x y x x y y

3 2 2 3

3 3 3 3 3 3

2 2 3

3 3 3 3 3 3 3

0 0 20 12 6 2 0 0

0 0 0 0 1 0 0 2 2 0 0 3 4 3 0 0 4

x x y x y y

x y x x y y x

2 2 3

3 3 3 3 3

2 2 3 2 2 3

3 3 3 3 3 3 3 3 3 3 3 3

6 6 4 0

0 0 0 0 0 2 0 0 2 6 0 0 2 6 12 0 0 2 6 12 20

x y x y y


1

1

,

1

,

2 1

,

2 1

,

2 1

,

2

2

,

2

,

2 2

,

2 2

,

2 2

,

3

3

,

3

,

2 3

,

2 3

,

2 3

,

x

y

xx

xy

yy

x

y

xx

xy

yy

x

y

xx

xy

yy

L

L

L

L

L

L

L

L

L

L

L

L

L

L

L

c

Tenemos así 18 ecuaciones y 21 incógnitas:


Elemento P5 - Desarrollo del MEF para StokesLas tres ecuaciones faltantes se obtienen calculando la derivada normal en los puntos

medios de los lados:

Matricialmente, el proceso realizado puede escribirse:

2 2 3 2 2 3 4 3 2 2 3 4

12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12

12

0 1 0 2 0 3 2 0 4 3 2 0 5 4 3 2 0

0 0 1 0

x

y

n x y x x y y x x y x y y x x y x y x y y

n

2 2 3 2 2 3 4 3 2 2 3 4

12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12

2

23 23 23 23 23

2 0 2 3 0 2 3 4 0 2 3 4 5

0 1 0 2 0 3 2x

x y x x y y x x y x y y x x y x y x y y

n x y x x

2 3 2 2 3 4 3 2 2 3 4

23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23

2 2

23 23 23 23 23 23 23

0 4 3 2 0 5 4 3 2 0

0 0 1 0 2 0 2 3 0y

y y x x y x y y x x y x y x y y

n x y x x y y

3 2 2 3 4 3 2 2 3 4

23 23 23 23 23 23 23 23 23 23 23 23 23 23

2 2 3 2

31 31 31 31 31 31 31 31 31 31 31

2 3 4 0 2 3 4 5

0 1 0 2 0 3 2 0 4 3 2x

x x y x y y x x y x y x y y

n x y x x y y x x y x y

2 3 4 3 2 2 3 4

31 31 31 31 31 31 31 31 31 31

2 2 3 2 2 3

31 31 31 31 31 31 31 31 31 31 31 31 31

0 5 4 3 2 0

0 0 1 0 2 0 2 3 0 2 3 4 y

y x x y x y x y y

n x y x x y y x x y x y y

12

,

23

,

31

,

4 3 2 2 3 4

31 31 31 31 31 31 31 310 2 3 4 5

n

n

n

L

L

L

x x y x y x y y

c

1 Ac LΦ c A LΦ

1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 12 23 31

, , , , , , , , , , , , , , , , , ,

T

x y xx xy yy x y xx xy yy x y xx xy yy n n n Φ

00 10 01 20 11 02 30 21 12 03 40 31 22 13 23 14 05

Tc c c c c c c c c c c c c c c c cc

2 2 2 2 2 2 2 2 2diag [1 1 1 ]L L L L L L L L L L L L L L L L L LL


El parámetro de escala L se elige de forma de condicionar adecuadamente la

matriz A y así evitar problemas al invertirla.


10-4

10-2

100

102

104

102

104

106

108

1010

1012

1014

1016

1018

Fact

cond(A

)

2 2 2 2 2 2

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

1max , ,

2L x x y y x x y y x x y y

La figura muestra la variación de la

condición de A para un triángulo arbitrario

(coordenadas aleatorias) en función de un

parámetro de escala igual a Fact*L



i i i iv u dx v f dxW W

2

1

1 21

2

2

22

2 22

u

x yxu

u

y yy

u

x xxu

u

y x y

u

2 2 2 2 2 2 2 2 2 2

22 2

20 30 212 2 2 2

22 2 3 2 2 36 6 20 62 2 12 2 12 2

2

2

1 2 6 2

0 0 0 0 0 0 0 0 0 0 0

i j

ij

T

L L L L L L L L L L

i ic x y c c x c x y

x L L L L L

x y x x y y x x y x y yx

x y

pc c

2 2 2 2 2 2 2 2 2 2

2 2 2 2

22 2 3 2 2 33 3 6 61 2 2 4 4 4

2

2

62 2 2

2

0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0

T

L L L L L L L L L L

L L L L

x y x x y y x x y x y yx y

x yy

pc c

2 2 2 2 2 2

22 2 3 2 2 36 6 2012 2 12

2 0 0

T

L L L L L Lx x y y x x y x y y

y

pc c

2 2

00 10 01 20 11 02

0 5

( , ) , , ( , ) K.i j T

ij

i j


p c



2 2 2 2 2 2

1 1 2 2 2 2 2 22 Wi i hv u dx v u v u dx dx

x x x y x y y y

W W W

2 2 2 2 2 21 1

2 2 2 2

18 18

no nulos últimos 18 18 (primeros tres términos, constante, lineal en x,y, nulos)

2T T T

T Tdxx x x y x y y y

W

0 0p p p p p p

K LA A L LA A L0 P

2 2 2

2 2 2

2 2 2

2 2 2

2 2 2

T T T

T T T

T T T

x x x

x y x y x y

y y y

p pc Ψ LA

p pc Ψ LA

p pc Ψ LA

2 2 21

2 2 2

2 2 21

2 2 21

2 2 2

T T

T T

T T

x x x

x y x y x y

y y y

p pc A LΦ

p pc A LΦ

p pc A LΦ

Usando:

y teniendo en cuenta que los parámetros son arbitrarios, obtenemos:Ψ


Si introducimos la hipótesis cinemática siguiente (cálculo de derivadas en el

punto medio a partir de las derivadas en los nodos vértice):

podemos luego expresar la matriz de rígidez de 18x18 (con grados de libertad

sólo en los vértices) como:

siendo:


12 12 12 121 1 2 2

, , , ,12

, 12 12 12 12

23 23 23 2323 2 2 3 3

, , , , , 2

31

,31 31 31 313 3 1 1

, , , ,

2 2 2 20 0 0 0 0 0 0 0 0 0 0 0 0 0

10 0 0 0

2 2 2 2 2

2 2 2 2

x y x y

x y x y

n x x y y

x y x y

n x y x y

nx y x y

x y x y

n n n n

n n n nn n n n

n

n n n n

3 23 23 23 1:18

31 31 31 31

12

,

23

, 1:18

31

,

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0

x x y y

x x y y

n

n

n

n n n

n n n n

Φ

DΦ

18 181

(4:21,1:21)18 18 3 18

IB A

D

TK B P B

Ejemplo: Flujo en una Cavidad Cuadrada


-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.60

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.40

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1) Elemento sin escalado ni traslacion

Malla 10x10x2 triángulos Malla 30x30x2 triángulos



-0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.40

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.40

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1) Elemento sin escalado ni traslacion




-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.60

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.40

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

2) Elemento c/escalado




-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.40

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.40

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1


2) Elemento c/escalado

Función línea de corriente (70 x 70)



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

-2

-1

0

1

2

3

4

5

6

7

8

x 10-4


Flexión de placas elásticas• Consideremos una delgada placa elástica P,

cuya superficie media está dada por el dominio W2, sujeta a carga transversal f.

• Se desea determinar

– deflexión transversal u

– momentos Mij, i,j=1,2.

• Los momentos están definidos:

• Convención de signos para rotaciones

x1

x2

fdx

W

u

t

n

Gl (lib

re) G

c (empotrado)

Gs (simpl. apoyado)

2

2

h

ij ijhM z dzs

11 22,M M Momentos flectores

12 21,M M Momentos de “twist”

1 2 ,1

2 1 ,2

ˆ

ˆ

u

u

Flexión de placas elásticas (cont.)

• El problema está gobernado por


,

c

s

l

en Ecuación de equilibrio

0 sobre CB empotrado

0 sobre CB simpl. apoyado

( ) 0 sobre CB libre

ij ij

nn

nn

M f

uu

n

u M

M Q

W

G

G

GM

,( ) Fuerza de corte o transversalntij j i

MQ M n

t

M


Flexión de placas elásticas (cont.)

• Ecuaciones de clausura

– Asumiendo pequeñas deflexiones y material elástico lineal, la ecuación

constitutiva (ley de Hooke) toma la forma

– Las constantes y dependen del módulo de elasticidad E y del coef. de

Poisson , así como del espesor de la placa h, de acuerdo a

2

,

, constantes

curvaturaij ij

i j

uu

x x

ij ij ijM u

3 3

212(1 ) 12(1 )

Eh Eh


Formulación variacional del problema de flexión de placas elásticas

1. Adoptamos el espacio

2. Hacemos

2

c sV : H ( ), 0 en , 0 en v

v v v vn

W G G

,

, , ,

, , ,

,( )

( )

ij ij

ij j i ij j i

ij ij ij j i ij j i

ij ij nn nt ij j i

ntij ij nn ij

fvdx M vdx

fvdx M v dx M n vds

fvdx M v dx M n v ds M n vds

v vfvdx M v dx M ds M ds M n vds

n t

Mvfvdx M v dx M ds vds M

n t

W W

W W G

W W G G

W W G G G

W W G G

,

,

0 0

( )

( ) ( ) ( )

j i

ntij ij nn ij j i

ij ij ij ij ij

n vds

Mvfvdx M v dx M ds M n vds

n t

fvdx u v dx u v u v dx

G

W W G G

W W W

,ij j i ij i j ij i j

nn nt

v vM n v M n n M t n

n tv v

M Mn t

,i i i

v vv n t

n t

Integración

por partes

(con G suave)

a( , )u vL( )v


Formulación variacional del problema de flexión de placas elásticas

(cont.)

• La forma variacional del problema de flexión de placas elásticas resulta

– La forma bilineal a(.,.) es en general simétrica y continua.

– Además, a(.,.) es V-elíptica si Gc0, i.e., si la placa está empotrada a lo largo

de una parte de su borde.

– La forma lineal L(.) es continua.

2

c s

Hallar V / a( , ) L( ), V

con: a( , )= ( ) ( )

L( )=

V : H ( ), 0 en , 0 en

( )

ij ij

u u v v v

u v u v u v dx

v fvdx

vv v v v

n

V

W

W

W G G

Formulación variacional del problema de flexión de placas

elásticas (cont.)

• Definiendo:

expresamos la forma bilineal matricialmente:

• Ahora se puede formular el MEF para el problema de flexión de placas

elásticas usando el elemento C1-continuo ya descrito.


2

2

1

2

2

2

2

1 2

2

u

x

uu

x

u

x x

κ

3

2

0 1 1 0

0 1 1 012 1

0 00 0 22

Eh

D

a( , )=u v u v dxW

κ Dκ


MEF no conforme aplicado al problema de Stokes

Consideremos el problema de Stokes en W2. Luego:

• Si W es simplemente conexo (i.e., no contiene agujeros), div v=0 en W si y solo si

o sea:

• Anteriormente, usamos un subespacio Wh de dimensión finita de (usamos

el elemento finito C1-continuo) y definimos

• Probaremos ahora un elemento cuadrático NO CONFORME (o sea, no verifica

que Vh !!!!

21 1 2

1 2 0

1 2

V : ( , ) H ( ) y 0 en v v

v vx x

W W

v v

2 1

, rot para alguna función .

: función de corriente del campo de velocidades .

x x

v

v

2

0V rot , H ( ). Wv v

2

0H ( )W

V : rot , W .h h v v

2

0H ( )W


Elemento de Morley - Desarrollo del MEF para Stokes

2 2

00 10 01 20 11 02

0 2

( , ) , , ( , ) K.i j T

ij

i j


p c

2 2

00 10 01 20 11 02

0 5

( , ) , 1,2,3.i j k

k k ij k k k k k k k k

i j

x y c x y c c x c y c x c x y c y k

2 2

00 10 01 20 11 02, 1TT x y x y x x y y c c c c c c p c

2 21 1,2,3.k

k k k k k kx y x x y y k c

Consiste en un EF triangular con polinomios cuadráticos. El campo será

aproximado por

( , )x y

con:

Para el cálculo de los coeficientes c expresamos primero el valor del campo en los nodos


0 0

0 0

; ; : longitud característica de elemento

, : centroide del elemento

x x y yx y L

L L

x y

2 2 2 2 2 2

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

1max , ,

2L x x y y x x y y x x y y



Matricialmente, el proceso realizado puede escribirse:

1 Ac LΦ c A LΦ

1 2 3 12 23 31

, , , 00 10 01 20 11 02

TT

n n n c c c c c c Φ c

i i i iv u dx v f dxW W

2

1

1 21

2

2

22

2 22

u

x yxu

u

y yy

u

x xxu

u

y x y

u

Las ecuaciones restantes se obtienen calculando la derivada normal en los puntos medios de

los lados:

12 12 12 12 12 12

23 23 23 23 23 23

31 31 31 31 31 31

0 1 0 2 0 0 0 1 0 2

0 1 0 2 0 0 0 1 0 2

0 1 0 2 0 0 0 1 0 2

x y

x y

x y

n x y n x y

n x y n x y

n x y n x y

12

,

23

,

31

,

n

n

n

L

L

L

c

diag [1 1 1 ]L L LL



2

2

202 2 2

2

112 2

2

022 2 2

1 2 20 0 0 0 0

1 10 0 0 0 0

2 20 0 0 0 0

i j

ij

i ic x y c

x L L L L

cx y L L

cy L L

c

c

c

2 2 2 2 2 2


x x x y x y y y

W W W

2 2

00 10 01 20 11 02

0 2

( , ) , , ( , ) K.i j T

ij

i j


p c

1

1

1

( ', ')

1 20 0 0

,1 2

0 0 0

ji

ij

i j

ij

x y

i x yc x y

Ly L L Lx y

i x yc x y

L L L Lx

G

u c G A LΦ

Las velocidades se calculan

en post-tratamiento (no son

variables primales). Ejemplo:

en el centro del triángulo




2 2 2 2 2 2


x x x y x y y y

W W W

2 2 2 2 2 21 1

2 2 2 2 2

no nulos últimos 3 3 términos

4 0 02

0 2 0

0 0 4

T T TT T

Kd Ax x x y x y y y L

W

0 0

p p p p p pK LA x A L LA A L

0

2 2 2

2 2 2

2 2 2

2 2 2

2 2 2

T T T

T T T

T T T

x x x

x y x y x y

y y y

p pc Ψ LA

p pc Ψ LA

p pc Ψ LA

2 2 21

2 2 2

2 2 21

2 2 21

2 2 2

T T

T T

T T

x x x

x y x y x y

y y y

p pc A LΦ

p pc A LΦ

p pc A LΦ

Usando:

y teniendo en cuenta que los parámetros son arbitrarios, obtenemos:Ψ

donde es el área del elemento.KA

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.60

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1


Ejemplo elemento de Morley – Cavidad cuadrada

Malla 10x10x2 triángulos



-0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.40

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Malla 70x70x2 triángulos



Malla 70x70x2 triángulos – Función línea de corriente0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

-9

-8

-7

-6

-5

-4

-3

-2

-1

0

1

x 10-4

Desarrollo subelemento 10gdl


3k

1

ˆ( ) ( ) k

kj j j

x Fξ ξ x

11 2 31

21 2 32

3

ˆ( )

ˆ ˆ( )( )

ˆ

F x x x

F y y y

ξF x

ξ

3

k

1

ˆˆ( ) ( ) , Kk

k

x F ξ ξ a ξ

1 1

1 1 1 21 3 2 3

1 2 3 1 2 31 2 2 2

1 2 3 1 2 3 1 3 2 32 2 1 2

1 2 3 3

1 2

ˆ ˆ

1 0ˆ ˆ

( ) ( ) 0 1

1 1ˆ ˆ

F F

x x x xx x x x x x

F F y y y y y y y y y y

FJ ξ ξ

ξ

Aproximación de la geometría: polinomios lineales

1 1

2 2

3 1 2

ˆ

ˆ

ˆ 1

2 3 2 3

1 11

1 3 1 32 2

1 1( )

2 2

y y x x b a

b ay y x x

ξJ ξ

xdonde

(doble del área del

triángulo)

2 det J

Notar: 1 3 2 3 2 3 1 3

2 1 1 22 det x x y y x x y y a b a b J

Desarrollo subelemento 10gdl (cont)


2 2

1 2 00 10 1 01 2 20 1 11 1 2 02 2 1 2 1 2

0 3

ˆ( ) , , ( , ) Ki j T

ij

i j

w c c c c c c c

ξ p c

Aproximación del campo: polinomio cúbico completo

2 2 3

1 2 1 2 1 1 2 2 1

00 10 01 20 11 02 30 21 12 03

, 1T

Tc c c c c c c c c c

p

c

2 2

1 2 00 10 1 01 2 20 1 11 1 2 02 2

0 3

( ) ( ) ,

1,2,3.

i jk k k k k k k k k

k k ij

i j

w w c c c c c c c w

k

x ξ

1

2

3

1 1 0 1 0 0 1 0 0 0

1 0 1 0 0 1 0 0 0 1

1 0 0 0 0 0 0 0 0 0

w

w

w

c

Para el cálculo de los coeficientes c expresamos el valor del campo en los nodos


donde1 2 3

1 0 0; ; .

0 1 0

ξ ξ ξ

Calculando luego el valor de la derivada respecto de x en los

nodos, logramos tres ecuaciones más:1

, 1 2

2

1 1( )

2 2kk k k

T T

k

x

bw w ww b b

bx

x xξ ξ ξ

ξx

ξ ξ ξ


1

1 2 10 20 1 11 2

0 31 1

1 2

1

1 2 01 11 1 02 2

0 32 2

( ) 2

ˆ, ( , ) K

ˆ( ) 2

Ti j

ij

i j

Ti j

ij

i j

wi c c c c

wj c c c c

pξ c

px c

2 21 1 2 1 1 2 2

2 2

1 2 1 1 2 2

2

0 1 0 2 0 3 2 0

0 0 1 0 2 0 2 3

w

w

c

Evaluando ambas derivadas (respecto de x y de y) en los nodos:

1 2 31 0 0

; ; .0 1 0

ξ ξ ξ

Obtenemos seis ecuaciones con 10 incógnitas:


1

,1 2

1

,1 2

2

,1 2

2

,1 2

3

,1 2

3

,1 2

0 0 0 0 0 1 0 2 0 0 3 0 0 0

0 0 0 0 0 0 1 0 1 0 0 1 0 0

0 0 0 0 0 1 0 0 1 0 0 0 1 01

0 0 0 0 0 0 1 0 0 2 0 0 0 32

0 0 0 0 0 1 0 0 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0 0 0

x

y

x

y

x

y

wb b

wa a

wb b

wa a

wb b

wa a

c

1

,1 2 1 2 1 2

1

,1 2 1 2 1 2

2

,1 2 1 2 1 2

2

,1 2 1 2 1 22 1 1 2

3

,1 2

3

,1 2

0 2 0 3 0 0

0 2 0 3 0 0

0 0 2 0 0 31

0 0 2 0 0 3

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

x

y

x

y

x

y

wb b b b b b

wa a a a a a

wb b b b b b

wa a a a a aa b a b

wb b

wa a

c


Finalmente, evaluando la derivada respecto de la normal en el punto

intermedio entre los nodos 1 y 2:

12 12

12 12

, 1 2 1 2 1 1 2 2 ,

1 1( )

2 2n x y x y x y n

w ww n b b n a a n b n a n b n a w

ξ ξ

xξ ξ

1

1 12 ,2 2

1 3 1 10 1 0 1 0 0

2 4 2 4

1 1 1 30 0 1 0 1 0

2 4 2 4

w

w

c

2 22 1

12 2 1 2 1

1 212

1x

y

n y yx x y y

n x x

n

Obtenemos así una décima ecuación:

121 2 1 2 1 21 2 1 2 1 2 ,

2 1 1 2

1 3 30

2 4 2 4 4 2 4n

d d d d d dd d d d d d w

a b a b

c

donde:

1 1 1 2 2 2; .x y x yd n b n a d n b n a


Reordenando ecuaciones, tenemos el sistema:

1 2 1 2 1 2

1 2 1 2 1 2

1 2 1 2 1 2

1 2 1 2 1 2

1 2

1 2

1 2 1 2 1 21 2 1 2 1 2

2 2 0 2 0 0 2 0 0 0

0 2 0 3 0 0

0 2 0 3 0 0

2 0 2 0 0 2 0 0 0 2

0 0 2 0 0 31

0 0 2 0 0 32

2 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

3 30

2 4 2 4 4 2 4

b b b b b b

a a a a a a

b b b b b b

a a a a a a

b b

a a

d d d d d dd d d d d d

A c

1

1

,

1

,

2

2

,

2

,

3

3

,

3

,

12

,

x

y

x

y

x

y

n

w

w

w

w

w

w

w

w

w

w

c

La matriz A puede invertirse analíticamente y calcular así las

funciones de forma: 1

1 2 1 2, ,T φ p A


Se obtiene:

2

1 1 2 3 1 1 3 1 2 3

2

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

2

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

2

4 1 2 3 2 2 3 1 2 3

2

5 1 2 3 2 1 3 3 1 3 3 2 1 2 3

2

6 1 2 3 2 1 3 3 1 3 3

, , 3 2 6

, ,

, ,

, , 3 2 6

, ,

, ,

a a a a

b b b b

a a a a

b b b

2 1 2 3

2

7 1 2 3 3 3

2

8 1 2 3 3 2 1 1 2

2

9 1 2 3 3 2 1 1 2

10 1 2 3 3 1 2 3

, , 3 2

, ,

, ,

, , 4

b

a a

b b

H

donde: 1

3

1 2

23

1 2

3

1 2

;

;

1

d

d d

d

d d

Hd d

y:

3 1 21

Se adjunta a continuación código de

manipulación simbólica para cálculo

de estas funciones de forma.

clear

syms x1 x2 x3 y1 y2 y3 real

syms xi1 xi2 positive

phih1 = xi1;

phih2 = xi2;

phih3 = 1 - xi1 - xi2;

AA = [x1 x2 x3

y1 y2 y3];

Dphih = [eval(diff (phih1,xi1)) eval(diff (phih1,xi2));

eval(diff (phih2,xi1)) eval(diff (phih2,xi2));

eval(diff (phih3,xi1)) eval(diff (phih3,xi2)); ];

J = AA * Dphih;

J1 = inv(J);

% b1/(a1*b2-a2*b1) <- J1(1,1);

% a1/(a1*b2-a2*b1) <- J1(1,2);

% b2/(a1*b2-a2*b1) <- J1(2,1);

% a2/(a1*b2-a2*b1) <- J1(2,2);

syms b1 a1 b2 a2 n1 n2 real

p = 1;

for i=1:3

for j=0:i

p = [p xi1^(i-j)*xi2^(j)];

end

end

%

dpxi1 = subs(pxi1,{xi1,xi2},{0,1});


C(5,:) = b1/(a1*b2-a2*b1)*dpxi1 + b2/(a1*b2-a2*b1)*dpxi2;



C(8,:) = b1/(a1*b2-a2*b1)*dpxi1 + b2/(a1*b2-a2*b1)*dpxi2;

%

% dw/dx2 (ak1,ak2) = [a1 a2]* dw/dxh (xi1^k,xi2^k) = dwx1_k k=1,3

%

pxi1 = diff(p,xi1);

pxi2 = diff(p,xi2);



C(3,:) = a1/(a1*b2-a2*b1)*dpxi1 + a2/(a1*b2-a2*b1)*dpxi2;







%

% dw/dn = [n1 n2]*[ dw/dx1(a^{23}) dw/dx2(a^{23}] =

% l23 = sqrt((y3-y2)^2 +(x3-x2)^2);

% n = [(y3-y2); -(x3-x2)]/l23;

dpxi1 = subs(pxi1,{xi1,xi2},{1/2,1/2});

dpxi2 = subs(pxi2,{xi1,xi2},{1/2,1/2});

C(10,:) = [0, c1, c2, c1, (c1+c2)/2, c2, 3*c1/4, c1/2 + c2/4, c1/4 + c2/2, 3*c2/4];

Cinv = inv(C);

phi1 = p*Cinv(:,1);

phi2 = p*Cinv(:,2);

phi3 = p*Cinv(:,3);

phi4 = p*Cinv(:,4);

phi5 = p*Cinv(:,5);

phi6 = p*Cinv(:,6);

phi7 = p*Cinv(:,7);

phi8 = p*Cinv(:,8);

phi9 = p*Cinv(:,9);

phi10 = p*Cinv(:,10);

syms mu3 lam3 xi3 b3 a3 H3

phi1Fel = xi1^2*(3-2*xi1) + 6*mu3*xi1*xi2*xi3;

phi1Fel = subs(phi1Fel,mu3,c1/(c1+c2));

phi1Fel = subs(phi1Fel,xi3,1-xi1-xi2);

errphi1 = simplify(phi1-phi1Fel)

phi2Fel = - xi1^2*(a3*xi2-a2*xi3) - (a3*mu3-a1)*xi1*xi2*xi3;

phi2Fel = subs(phi2Fel,a3,-a1-a2);

phi2Fel = subs(phi2Fel,mu3,c1/(c1+c2));



phi3Fel = xi1^2*(b3*xi2-b2*xi3) + (b3*mu3-b1)*xi1*xi2*xi3;

phi5Fel = - xi2^2*(a1*xi3-a3*xi1) + (a3*lam3-a2)*xi1*xi2*xi3;

phi5Fel = subs(phi5Fel,a3,-a1-a2);

phi5Fel = subs(phi5Fel,lam3,c2/(c1+c2));



phi6Fel = xi2^2*(b1*xi3-b3*xi1) - (b3*lam3-b2)*xi1*xi2*xi3;

phi6Fel = subs(phi6Fel,b3,-b1-b2);

phi6Fel = subs(phi6Fel,lam3,c2/(c1+c2));



phi7Fel = xi3^2*(3-2*xi3);



phi8Fel = - xi3^2*(a2*xi1-a1*xi2);



phi9Fel = xi3^2*(b2*xi1-b1*xi2);



phi10Fel = 4*H3*xi1*xi2*xi3;


Date post:	06-Feb-2020
Category:	Documents
Upload:	others
View:	2 times
Download:	0 times

Introducción al Método de los Elementos Finitos• Para llevar el problema de Stokes a la forma...

Documents