Post on 06-Feb-2020
transcript
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
Introducción al Método de los Elementos Finitos 3
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
Introducción al Método de los Elementos Finitos 4
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
Introducción al Método de los Elementos Finitos 5
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
Introducción al Método de los Elementos Finitos 6
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
Introducción al Método de los Elementos Finitos 7
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
Introducción al Método de los Elementos Finitos 8
MEF aplicado al problema de Elasticidad
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
Introducción al Método de los Elementos Finitos 9
MEF aplicado al problema de Elasticidad
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)
Introducción al Método de los Elementos Finitos 11
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
Introducción al Método de los Elementos Finitos 12
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
Introducción al Método de los Elementos Finitos 13
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
Introducción al Método de los Elementos Finitos 14
, ,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
Introducción al Método de los Elementos Finitos 15
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:
Elemento P5 - Desarrollo del MEF para Stokes
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.
Introducción al Método de los Elementos Finitos 16
( 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
Introducción al Método de los Elementos Finitos 17
Elemento P5 - Desarrollo del MEF para Stokes
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
Introducción al Método de los Elementos Finitos 18
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 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
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
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
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
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
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
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
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 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
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
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
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
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
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
x y x x y y x 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:
Introducción al Método de los Elementos Finitos 19
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
Elemento P5 - Desarrollo del MEF para Stokes
El parámetro de escala L se elige de forma de condicionar adecuadamente la
matriz A y así evitar problemas al invertirla.
Introducción al Método de los Elementos Finitos 20
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
Introducción al Método de los Elementos Finitos 21
Elemento P5 - Desarrollo del MEF para Stokes
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
x y c x y c c x c y c x c x y c y x y x y
p c
Introducción al Método de los Elementos Finitos 22
Elemento P5 - Desarrollo del MEF para Stokes
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:Ψ
Elemento P5 - Desarrollo del MEF para Stokes
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:
Introducción al Método de los Elementos Finitos 23
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
Introducción al Método de los Elementos Finitos 24
-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
Ejemplo: Flujo en una Cavidad Cuadrada
Introducción al Método de los Elementos Finitos 25
-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
Malla 70x70x2 triángulos Malla 100x100x2 triángulos
Ejemplo: Flujo en una Cavidad Cuadrada
Introducción al Método de los Elementos Finitos 26
-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
Malla 10x10x2 triángulos Malla 30x30x2 triángulos
Ejemplo: Flujo en una Cavidad Cuadrada
Introducción al Método de los Elementos Finitos 27
-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
Malla 70x70x2 triángulos Malla 100x100x2 triángulos
2) Elemento c/escalado
Función línea de corriente (70 x 70)
Introducción al Método de los Elementos Finitos 28
Ejemplo: Flujo en una Cavidad Cuadrada
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
Introducción al Método de los Elementos Finitos 29
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
Introducción al Método de los Elementos Finitos 30
,
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
Introducción al Método de los Elementos Finitos 31
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
Introducción al Método de los Elementos Finitos 32
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
Introducción al Método de los Elementos Finitos 33
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.
Introducción al Método de los Elementos Finitos 34
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κ
Introducción al Método de los Elementos Finitos 35
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
Introducción al Método de los Elementos Finitos 36
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
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
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
Obtenemos así tres ecuaciones con 6 incógnitas:
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
Introducción al Método de los Elementos Finitos 37
Elemento de Morley - Desarrollo del MEF para Stokes
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
Introducción al Método de los Elementos Finitos 38
Elemento de Morley - Desarrollo del MEF para Stokes
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
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
00 10 01 20 11 02
0 2
( , ) , , ( , ) 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
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
Introducción al Método de los Elementos Finitos 40
Elemento de Morley - Desarrollo del MEF para Stokes
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 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
Introducción al Método de los Elementos Finitos 41
Ejemplo elemento de Morley – Cavidad cuadrada
Malla 10x10x2 triángulos
Introducción al Método de los Elementos Finitos 42
Ejemplo elemento de Morley – Cavidad cuadrada
-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
Introducción al Método de los Elementos Finitos 43
Ejemplo elemento de Morley – Cavidad cuadrada
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
Introducción al Método de los Elementos Finitos 55
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)
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
Obtenemos así tres ecuaciones con 10 incógnitas:
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
ξ ξ ξ
Desarrollo subelemento 10gdl (cont)
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:
Desarrollo subelemento 10gdl (cont)
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
Desarrollo subelemento 10gdl (cont)
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
Desarrollo subelemento 10gdl (cont)
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
Desarrollo subelemento 10gdl (cont)
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});
dpxi2 = subs(pxi2,{xi1,xi2},{0,1});
C(5,:) = b1/(a1*b2-a2*b1)*dpxi1 + b2/(a1*b2-a2*b1)*dpxi2;
dpxi1 = subs(pxi1,{xi1,xi2},{0,0});
dpxi2 = subs(pxi2,{xi1,xi2},{0,0});
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);
dpxi1 = subs(pxi1,{xi1,xi2},{1,0});
dpxi2 = subs(pxi2,{xi1,xi2},{1,0});
C(3,:) = a1/(a1*b2-a2*b1)*dpxi1 + a2/(a1*b2-a2*b1)*dpxi2;
dpxi1 = subs(pxi1,{xi1,xi2},{0,1});
dpxi2 = subs(pxi2,{xi1,xi2},{0,1});
C(6,:) = a1/(a1*b2-a2*b1)*dpxi1 + a2/(a1*b2-a2*b1)*dpxi2;
dpxi1 = subs(pxi1,{xi1,xi2},{0,0});
dpxi2 = subs(pxi2,{xi1,xi2},{0,0});
C(9,:) = 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));
phi2Fel = subs(phi2Fel,xi3,1-xi1-xi2);
errphi2 = simplify(phi2-phi2Fel)
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));
phi5Fel = subs(phi5Fel,xi3,1-xi1-xi2);
errphi5 = simplify(phi5-phi5Fel)
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));
phi6Fel = subs(phi6Fel,xi3,1-xi1-xi2);
errphi6 = simplify(phi6-phi6Fel)
phi7Fel = xi3^2*(3-2*xi3);
phi7Fel = subs(phi7Fel,xi3,1-xi1-xi2);
errphi7 = simplify(phi7-phi7Fel)
phi8Fel = - xi3^2*(a2*xi1-a1*xi2);
phi8Fel = subs(phi8Fel,xi3,1-xi1-xi2);
errphi8 = simplify(phi8-phi8Fel)
phi9Fel = xi3^2*(b2*xi1-b1*xi2);
phi9Fel = subs(phi9Fel,xi3,1-xi1-xi2);
errphi9 = simplify(phi9-phi9Fel)
phi10Fel = 4*H3*xi1*xi2*xi3;
phi10Fel = subs(phi10Fel,xi3,1-xi1-xi2);