• Author / Uploaded
• idiazm

• Categories
• Algoritmos
• Álgebra
• Análisis numérico
• Objetos matemáticos
• Análisis

Citation preview

IEG 3110 Elementos Finitos Lineales Clase 05

a ´olic

ile h eC

d

t

Daniel Hurtado, Ph.D. Profesor Asistente

i

n U a i

a dC

a

id s r ve

Departamento de Ingenier´ıa Estructural y Geot´ ecnica Pontificia Universidad Cat´ olica de Chile

©

4 201

el i n a

do a t ur

fic i t n o3/Junio/2014 ,P

H

D

Daniel Hurtado

IEG 3110

1 / 11

Integraci´ on num´erica

ile h eC

d a c i Recordamos que las componentes de la matriz de rigidez y vector de fuerzas t´ol de un elemento est´ a expresado en t´ ermino de integrales a CZ Z Z d a iddv + N (x)¯t ds, k = B (x)CB (x) dv y f = N (x)b s r ive n Para el caso de Tri3, la integral se puede resolver a Uanal´ıticamente (B es constante, N es i c lineal) fi nti Quad4, Quad8, Quad9 las expresiones son Para el caso de elementos isoparam´ eotricos P bastante m´ as complicadas , o d Nos interesa aproximar elta valor de las integrales involucradas mediante m´ etodos simples, r u pero controlando elH error num´ erico: Integraci´ on Num´ erica el i n Da 4 201 e

eT

e

e

Ωe

eT

eT

∂Ωe t

Ωe

e

e

1

© 1

tambi´ en conocida como Cuadratura Num´ erica Daniel Hurtado

IEG 3110

2 / 11

Integraci´ on num´erica: f´ ormula de Newton-Cotes Sea f (x) una funci´ on continua en [a, b]. Definimos

ile h eC

b

Z I(f ) :=

f (x) dx

a ´olic

a

d

t

a dC

Como f (x) puede ser dif´ıcil de integrar, consideramos la alternativa de aproximar f (x) por un polinomio pn (x) con n el orden del polinomio, y luego integrar el polinomio aproximante, es decir, Z

a d i s r I(f ) ≈ I (f ) := p i(x) vedx n a U Usando interpolaci´on Lagrangeana Lo anterior se conoce como la formula de Newton-Cotes. i c para p (x), llegamos a ifi t n o X w f (x ) IP(f ) = , do a t urcuadratura y las constantes w son los pesos de cuadratura. Para donde x son los puntos de H estudiar el error num´ eelrico asociado a la cuadratura, definimos el error como i n Z X Da 4 E (f ) = f (x) − w f (x ) 1 20 b

n

h

a

n

n

h

k

k

k=0

k

k

n

b

©

n

k

a

Daniel Hurtado

k

k=0

IEG 3110

3 / 11

Integraci´ on num´erica: Regla trapezoidal y Simpson

ile h eC

Regla de Simpson

Regla Trapezoidal

a ´olic

d

t

a dC

a

i

Z

b

a

b−a f (x) dx ≈ [f (a) + f (b)] 2

do a t ur

n U a i

c tifi

Z

on ,P

©

H

D 4 Error en la estimaci´ on 1 20

f (x) dx ≈

a

b−a a+b [f (a) + 4f ( ) + f (b)] 6 2

Error en la estimaci´ on

(b − a)3 |En (f )| = max f (3) (ξ) ξ∈[a,b] 12 Daniel Hurtado

b

Puntos y pesos de cuadratura xk wk b−a a 6 a+b 2 (b − a) 2 3 b−a b 6

Puntos y pesos de cuadratura xk wk b−a a 2 b−a b 2

el i n a

id s r ve

|En (f )| = IEG 3110

(b − a)4 max f (4) (ξ) 196 ξ∈[a,b] 4 / 11

Integraci´ on num´erica: Cuadratura de Gauss La cuadratura de Gauss entrega f´ ormulas para integrar en forma exacta un polinomio de orden 2N − 1 con s´ olo N puntos de cuadratura.

a ´olic

Por ejemplo, consideramos el caso N = 2. Sea

ile h eC

d

t

f (x) = a3 x3 + a2 x2 + a1 x + a0

a dC

a d i s f (x) dx = w f (x ) + w efr(x ) niv U Entonces, evaluando y juntando t´ erminos llegamosaal sistema de ecuaciones fici i 2 =w + w (a ) t n o 0 =w x +, wPx (a ) o d 2/3 =w (a ) rtax + w x u (a ) l H0 =w x + w x e i Resolviendo obtenemos an que para integrar en forma exacta un polinomio c´ubico los puntos y pesos D de cuadratura son 4 201 k x√ w

y queremos que

Z

1

1

1

2

2

−1

1

2

0

1 1

2 2

2 1 1 3 1 1

2 2 2 3 2 2

©

1

2 3

k

1 2

Daniel Hurtado

k

−1/√3 +1/ 3

1 1

IEG 3110

5 / 11

Integraci´ on num´erica: Cuadratura de Gauss

a ´olic

ile h eC

d

t

N 1 2

©

4 201

el i n a

wk 2 1 1

i

n U a i

−c fi0 i t n o q ,P 3

do a t ur

xk 0√ −1/√3 +1/ q 3 3 5

3 5

a dC

a

id s r ve

5 9 8 9 5 9

H

D

Daniel Hurtado

IEG 3110

6 / 11

Integraci´ on num´erica: Cuadratura de Gauss 2D y 3D Para una integral 2D en dominio isoparam´ etrico utilizamos el teorema de Fubini y aplicamos cuadratura num´ erica sobre una integracion iterada Z

+1

Z

Z

+1

f (ξ, η)dξ dη ≈

f (ξ, η)dv = ˆe Ω

−1

−1

≈



NX ×N

a ´olic

f (ξk , ηl )wk wl

a

f (ξq , ηq )Wq

i

n U a i

Algunos ejemplos

id s r ve

c

Regla 1 × 1

do a t ur

(ξq , ηq )

el i n a

H

(0, 0)

ifi t n o ,WP

(ξq , ηq )

Wq

Regla 2 × 2 −1 (√ ,

q

3 ( √1 , 3 ( √1 , 3 ( √1 , 3

4

D 4 1 0

Daniel Hurtado

d

t

a dC

l=1 k=1

q=1

©2

N X N X

ile h eC

−1 √ ) 3 −1 √ ) 3 √1 ) 3 √1 ) 3

1 1 1 1

IEG 3110

7 / 11

Cuadratura num´erica en elementos finitos Usando cuadratura num´ erica, aproximamos las integrales involucradas en el c´ alculo de la matriz de rigidez y el vector de fuerzas, Z Z ke = BeT (x)CBe (x) dv = BeT (ξ)CBe (ξ) | det J(ξ)|dξ ˆe Ω

Ωe

≈

NX ×N

fe =

t

NeT (x)b dv =

q=1

i

id s r ve

ˆ eT (xi)b dξ N

n U a i

c

ifi t n o ,P

ˆ eT (ξq , ηq )b| det J(ξq , ηq )| Wq N

do a t ur

Algunas observaciones,

el i n a

Z ˆe Ω

NX ×N

a dC

a

Ωe

≈

d

BeT (ξq , ηq ) C Be (ξq , ηq ) | det J(ξq , ηq )| Wq

q=1

Z

a ´olic

ile h eC

H

Para elemento Quad4, la regla 1 × 1 no es suficiente, y entrega matrices con m´ as valores propios nulos que modos de cuerpo r´ıgido. Este fen´ omeno es conocido como subintegraci´ on, y debe evitarse.

©

4 201

D

Para el elemento Quad4, una regla 2 × 2 es adecuada (aunque no necesariamente exacta) Para elementos Quad8 y Quad9 se utiliza cuadratura Gaussiana de 3 × 3.

Daniel Hurtado

IEG 3110

8 / 11

Modos incompatibles

a ´olic

ile h eC

d

t

a dC

Recordamos que en flexi´ on, el elemento Quad4 tiene un mejor desempe˜ no que el Tri3.

a

id s r ve

Sin embargo, Quad4 presenta errores en la estimaci´ on del corte - fen´ omeno de corte paras´ıtico

i

n U a i

El corte paras´ıtico puede eliminarse usando elementos Quad8/9, pero con el costo de aumentar sustancialmente los grados de libertad

c

ifi t n o ,P

Una t´ ecnica muy utilizada en estos casos de flexi´ on predominante son los elementos con modos incompatibles. En particular, SAP2000 tiene este elemento en su librer´ıa.

do a t ur

La idea es enriquecer las funciones de forma para el campo de desplazamiento de Quad4, pero sin agregar nodos adicionales.

©

4 201

el i n a

H

D

Daniel Hurtado

IEG 3110

9 / 11

Modos incompatibles: Elemento de Wilson Q6 El punto de partida es la formulaci´ on del elemento Quad4. Entonces, agregamos funciones burbuja al campo de desplazamientos u ˆeh (ξ, η) =

4 X

ˆ e (ξ, η)ua + P1 (ξ, η)α1 + P2 (ξ, η)α3 N a

a=1

vˆeh (ξ, η) =

4 X

⇒

a=1

i

n U a i



P1 (ξ, η) P (ξ) = 0 e

fi(ξ,c η) i P t on 0 ,P

0 P2 (ξ, η)

do a t ur

1

d

t

a dC

ˆ eh (ξ) = Ne (ξ)ue + Pe (ξ)αe u

a

ˆ e (ξ, η)va + P1 (ξ, η)α2 + P2 (ξ, η)α4 N a

 donde P1 (ξ, η) = 1 − ξ 2 y P2 (ξ, η) = 1 − η 2 , αe = α1

a ´olic

ile h eC

id s r ve α2

α3

0 P2 (ξ, η)

α4

T

P1 (ξ, η) 0

0 P2 (ξ, η)



Luego, las deformaciones est´ an dadas por

donde

©

4 201

el i n a

D

H

ε(ueh ) = Be ue + Ge αe  ∂P

1

∂x

 Ge =  0

∂P1 ∂y

Daniel Hurtado

0 ∂P2 ∂y ∂P1 ∂x

IEG 3110

∂P2 ∂x

0 ∂P2 ∂y

0



∂P2  ∂y  ∂P2 ∂x

10 / 11

Modos incompatibles: Elemento de Wilson Q6 Finalmente del equilibrio del elemento obtenemos la ecuaci´ on matricial    e   e kuu kuα u f = kαu kαα αe 0 R donde kuu es la matriz de rigidez del Quad4,kuα = kαu = Ωe BeT CGe dv y R kαα = Ωe GeT CGe dv.

a ´olic

ile h eC

d

t

a dC

a

id s r ve

Entonces, aplicando condensaci´ on est´ atica obtenemos la matriz de rigidez del elemento Q6

i

n U a i

ke = kuu − kuα k−1 αα kαu

c

ifi t n o ,P

Adicionalmente, a la formulaci´ on anterior se le aplica un esquema de integraci´ on selectiva, de manera de que la formulaci´ on capture campos de deformaci´ on constante. El elemento con esta u ´ltima modificaci´ on se conoce como el elemento de Wilson-Taylor QM6.

do a t Para el ejemplo de viga en flexi´ ur on de la clase anterior tenemos H Deformaci´ on vertical el i n Quad4, 10 nodos, 4 elementos 0.607 a D QM6, 10 nodos, 4 elementos 1.009 4 1 0 Quad4, 27 nodos, 16 elementos 0.876 2

©

QM6, 27 nodos, 16 elementos

1.024

Error 41.1% 2.0% 15.0% 0.6%

Y vemos que la formulaci´ on de modos incompatibles mejora notablemente el comportamiento en flexi´ on. Daniel Hurtado

IEG 3110

11 / 11