PROBLEMAS RESUELTOS: MÉTODO DE LOS ELEMENTOS FINITOS Guillermo Rus Carlborg 1 , Esther Puertas García 2 Enero de 2008

1 Profesor

Contratado Doctor. Departamento de Mecánica de Estructuras. Universidad de Granada. 2 Profesor Ayudante. Departamento de Mecánica de Estructuras. Universidad de Granada.

c copyright 2007: Guillermo Rus Carlborg, Esther Puertas García ° Editor: Departamento de Mecánica de Estructuras, Universidad de Granada

G. Rus, E. Puertas

Problema 1 Se considera la viga empotrada en un extremo y sometida a axil p(x) representada en la figura. Empleando una discretización de dos elementos lineales y una discretización de un elemento cuadrático y suponiendo que la carga es constante p(x) = p0 y variable p(x) = pL0 x.

Se pide: 1. Plantear el problema teórico. 2. Discretizar y hallar las funciones de forma. 3. Obtener la matriz de rigidez. 4. Obtener el el vector de fuerzas externas. 5. Obtener el desplazamiento en el centro y extremo de la viga, comparando los resultados obtenidos para L = 10m, E = 0.1M P a, A = 0.01m2 , p0 = 0.1N/m.

1

G. Rus, E. Puertas

Solución 1 1. Planteamiento teórico del Problema Para obtener el problema a resolver basta con aplicar las ecuaciones de equilibrio en una rebanada de la viga:

N + dN + p(x)dx − N = 0 dN + p(x) = 0 dx Sabemos N =

Z

t/2

σxx b(z) dz = EA

−t/2

du dx

sustituyendo en la expresión anterior: d2 u dN + p(x) = EA 2 p(x) = 0 dx dx En consecuencia, la formulación fuerte del problema se puede escribir: Hallar u(x); x ∈ [0, L] tal que EA u,xx (x) + p(x) = 0; x ∈ (0, L) u(0) = 0; u,x (0) = 0 La formulación débil del problema consiste en aplicar un desplazamiento virtual v definido en [0, L] con las mismas condiciones de contorno e integrar en el dominio: Z L Z L − EA u,xx (x)v(x)dx = p(x)v(x)dx 0

0

Integrando por partes el primer miembro 1 Z L Z L u,xx (x)v(x)dx = u,x (x)v(x)|0 − 0

1

2

L

u,x (x)v,x (x)dx

0 L

Para las condiciones de contorno del problemam el producto u,x (x)v(x)|0 se anula.

G. Rus, E. Puertas Se deduce Z

L

EA u,x (x)v,x (x)dx =

0

EAu,x (x)v(x)|L0

+

Z

L

p(x)v(x)dx

0

En consecuencia, la formulación débil del problema, equivalente el Principio de los Trabajos Virtuales2 es: u ∈ H01 (0, L) tal que Z L L p(x)v(x)dx ∀v(x) ∈ H01 (0, L) EA u,x (x)v(x)dx = EAu,x (x)v(x)|L0 + Hallar

Z

0

0

Dada una partición uniforme de [0, L] en n intervalos de igual longitud, el problema discreto asociado sobre el espacio de elementos finitos construido sobre esta partición a partir de las funciones de forma Hi se define: Hallar u(x) ≃ H1 (x)u1 + H2 (x)u2 + · · · + HN (x)uN tal que Z x(e) n Z x(e) X j j ¯x(e) j (e) (e) ¯ EA u,x (x)v,x (x)dx = EAu,x (x)v(x) (e) + p(x)v(x)dx e=1

xi

(e)

xi

(e)

xi

Este planteamiento es análogo a resolver el sistema de ecuaciones:

Ku = f ¯ ¯ R1 (e) R L(e) (e) (e) donde Kij = 0 Bi C (e) Bj ¯J (e) ¯ A(e) dx′ ; fj = Fj + 0 p(x) Hj dx siendo K la matriz de rigidez, u el vector de desplazamientos de los nodos y f el ¯L (e) (e) ¯ vector de fuerzas externas, Fj = EAHj,x u Hj 0 . Para el caso general, el problema discreto consiste en XZ e

Ve

e Hallar ui (x1 , x2 , x3 ) = Hin (x1 , x2 , x3 )uni tal que Z X XZ e e n V e e Cijkl Bijnc Bklmd dV uc = fd Hdm dV + fdS Hdm dS e

Ve

e

Se

2. Discretización y funciones de forma Dos elementos lineales La discretización empleada mediante dos elementos lineales se recoge en la figura 1. Al definirse 3 nodos, existen 3 grados de libertad, de éstos los nodos 2 y 3 están definidos en desplazamientos, el único grado de libertad en fuerzas se define para el nodo 1. 2 En esta expresión el primer término representa el trabajo virtual interno que realizan las tensiones reales en la viga sobre las deformaciones virtuales. Y el segundo término es el trabajo virtual de las fuerzas exteriores sobre los desplazamientos virtuales

3

G. Rus, E. Puertas

Figure 1: Discretización de la estructura: (a) Nodos globales, (b) Nodos locales Las funciones de forma lineales se caracterizan porque toman el valor unidad en el nodo y cero en el resto de nodos. Para un elemento lineal, definido en coordenadas locales, se tiene: H1 (x) = 1 − x′ ;

H2 (x) = x′

Las derivadas en coordenadas locales de las funciones de forma son: dH1 = −1; dx′

dH2 =1 dx′

El jacobiano de la transformación entre coordenadas locales y globales para cada uno de los elementos es: J (e) =

dx = L(e) dx′

Las derivadas de las funciones de forma en coordenadas globales: dH1 dx′ 1 dH1 = = − (e) B1 = ′ dx dx dx L dH2 dH2 dx′ 1 B2 = = = (e) ′ dx dx dx L

Un elemento cuadrático Al tratarse de un elemento de tres nodos, las funciones de forma serán de tipo cuadrático. H1 (x) se tomará de forma que toma el valor 1 en el punto 1 y se anula en el resto; H2 (x) es la función que toma el valor unidad en el nodo 2 y H3 (x) la que tiene valor unidad en el nodo 3. 4

G. Rus, E. Puertas

(a)

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

1

2

3

4

5

(b)

1

2

3

(1)

4

5

(1)

Figure 2: Funciones de forma lineales para L = 10 (a) H1 ; (b) H2

1 H1 (x) = x′ (x′ − 1) 2 H2 (x) = (x′ + 1) (1 − x′ ) 1 H3 (x) = x′ (x′ + 1) 2 Las derivadas de las funciones de forma respecto a las coordenadas locales son: 1 dH1 ′ = x − dx′ 2 dH2 = −2x′ ′ dx dH3 1 = x′ + ′ dx 2 El jacobiano de la transformación entre coordenadas locales y globales para el elementos considerado es: J=

dx =L dx′

Las derivadas de las funciones de forma respecto a las coordenadas globales se obtendrán al realizar el producto entre las derivadas en coordenadas 5

G. Rus, E. Puertas 1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2 0.2

(a)

2

4

6

8

10

(b)

4

2

6

8

10

1 0.8 0.6 0.4 0.2

2

(c)

4

6

8

10

Figure 3: Funciones de forma cuadráticas para L = 10 (a) H1 ; (b) H2 ; (c) H3 . locales por el jacobiano de la transformación: ¶ µ dH1 dH1 dx′ 1 1 ′ B1 = = = x − dx dx′ dx 2 L ′ dH2 dx 1 dH2 = = −2x′ B2 = ′ dx dx dx µ ¶L dH3 dH3 dx′ 1 1 ′ B3 = = = x + dx dx′ dx 2 L

3. Matriz de Rigidez Dos elementos lineales Para el cálculo de la matriz de rigidez hay que obtener, en primer lugar, las matrices de rigidez elementales teniendo en cuenta los miembros ¯ (e) ¯ (e)de la R 1 que (e) (e) matriz se obtienen a partir de la integral Kij = 0 Bi C Bj ¯J ¯ A dx′ . Al ser los intervalos de igual longitud, con igual rigidez y área, las matrices de rigidez elementales son iguales para ambos elementos. K (e) 6

à R 1 B C (e) B = R01 1 (e) 1 B2 C B1 0

¯ (e) ¯ (e) ′ R 1 (e) ¯J ¯ A dx ¯ (e) ¯ (e) ′ R01 B1 C (e) B2 ¯J ¯ A dx B2 C B2 0

¯ (e) ¯ (e) ′ ! ¯J ¯ A dx ¯ (e) ¯ (e) ′ ¯J ¯ A dx

G. Rus, E. Puertas ¯ ¯ donde C (e) = E;A(e) = A;¯J (e) ¯ = L2 . K

(1)

=K

(2)

2EA = L

µ

1 −1 −1 1

¶

Una vez calculadas las matrices de rigidez elementales, se obtiene la matriz de rigidez global mediante el proceso de ensamblaje.   (1)   (1) K12 0 K11 1 −1 0 2EA   (1) (1) (2) (2)  −1 2 −1  K =  K12 K22 + K11 K12  = L (2) (2) 0 −1 1 0 K12 K22 Un elemento cuadrático El procedimiento a seguir es análogo. En este caso, al tener un único elemento, la matriz de rigidez de dimensión 3 × 3 se obtiene directamente, es decir, no es necesario el proceso de ensamblaje.   R1 R1 R1 ′ ′ ′ B C B |J| A dx B C B |J| A dx B C B |J| A dx 3 2 1 1 R0 1 R0 1   R0 K =  01 B2 C B1 |J| A dx′ 01 B2 C B2 |J| A dx′ 01 B2 C B3 |J| A dx′  R1 R1 R1 B3 C B1 |J| A dx′ 0 B3 C B2 |J| A dx′ 0 B3 C B3 |J| A dx′ 0 K=



EA  L

7 3 −8 3 1 3

−8 3 16 3 −8 3

1 3 −8 3 7 3

 

4. Vector de fuerzas Dos elementos lineales El vector de fuerzas elemental es elementos, se calcula ¯ para ¯ ambos R 1equivalente (e) (e) ¯ ′ v ′ ¯ dx . Se distinguen los dos mediante las ingrales fj = 0 f (x ) Hj J casos de carga planteados en el enunciado. Ã R ¯ (e) ¯ ′ ! 1 v ′ ¯J ¯ dx f (x ) H 1 ¯ (e) ¯ ′ f (e) = R01 v ′ ¯J ¯ dx f (x ) H 2 0 Carga constante (p(x) = p0 ) f

(1)

=f

(2)

p0 L = 4

µ

1 1

¶ 7

G. Rus, E. Puertas Mediante el proceso de ensamblaje se obtiene el vector de fuerzas global.     (1) f1 1 p0 L    (1) (2)  2 f =  f2 + f1  = 4 (2) 1 f2 ¡ Carga variable p(x) = f

(1)

p0 x L

p0 L = 24

¢

µ

1 2

¶

;

f

(2)

p0 L = 24

µ

4 5

¶

El vector de fuerzas global es 

   (1) f1 1 p0 L    (1) (2)  6 f =  f2 + f1  = 24 (2) 5 f2

Un elemento cuadrático En este caso se obtiene directamente el vector de fuerzas global.  R1  v ′ ′ f (x ) H |J| dx 1  R0  f =  01 f v (x′ ) H2 |J| dx′  R1 v ′ f (x ) H3 |J| dx′ 0 Carga constante (p(x) = p0 )

 1 p0 L   4 f= 6 1 ¡ ¢ Carga variable p(x) = pL0 x   0 p0 L   2 f= 6 1 

5. Desplazamiento de los nodos El desplazamiento de los nodos se obtiene directamente de la resolución del sistema de ecuaciones Ku = f . En primer lugar se simplifica el sistema ya que sabemos que el desplazamiento en el nodo inicial es nulo. 8

G. Rus, E. Puertas Dos elementos lineales Carga constante El sistema resultante es:  1     u 1 −1 0 1 p0 L   2EA  −1 2 −1   u2  = 2 L 4 3 u 0 −1 1 1

Se sabe que el desplazamiento de nodo del empotramiento u1 es nulo, por lo que se reduce el sistema eliminando la primera fila y columna: ¶µ 2 ¶ µ µ ¶ 2EA p0 L 2 u 2 −1 = u3 −1 1 1 L 4 Sustituyendo por los datos dados para las variables: ¶ ¶µ 2 ¶ µ µ 0.50 u 40 −20 = 0.25 u3 −20 20 Resolviendo el sistema se obtiene: u2 = 0.0375;

u3 = 0.05

Carga variable Para la carga variable, el sistema sólo cambia en su término independiente.  1     u 1 −1 0 1 p0 L   2EA  −1 2 −1   u2  = 6 L 24 3 u 0 −1 1 5

Eliminando la primera fila y colmuna se tiene ¶µ 2 ¶ µ µ ¶ 2EA p0 L 6 u 2 −1 = u3 −1 1 5 L 24 Sustituyendo por los valores de las variables ¶ ¶µ 2 ¶ µ µ 0.25 u 40 −20 = 0.208333 u3 −20 20 Resolviendo: u2 = 0.0229167;

u3 = 0.0333333 9

G. Rus, E. Puertas Un elemento cuadrático Para el caso en el que se emplee un elemento cuadrático, la resolución se realiza de forma análoga. Carga constante  7 −8 1   1    u 1 3 3 3 EA  −8 16 −8   2  p0 L   u 4 = 3 3 3 L 6 1 −8 7 3 1 u 3 3 3 µ 16 −8 ¶ µ 2 ¶ µ ¶ EA p0 L 4 u 3 3 = −8 7 1 u3 L 6 3 3 ¶ ¶µ 2 ¶ µ µ 0.666667 u 53.3333 −26.6667 = 0.166667 u3 −26.6667 23.3333 u2 = 0.0375;

u3 = 0.05

Carga variable 7 3 −8 3 1 3

1 3 −8 3 7 3

   0 u1 EA    u2  = p 0 L  2  L 6 1 u3 µ 16 −8 ¶ µ 2 ¶ µ ¶ EA p0 L 2 u 3 3 = −8 7 3 u 1 L 6 3 3 ¶ ¶µ 2 ¶ µ µ 0.333333 u 53.3333 −26.6667 = 0.166667 u3 −26.6667 23.3333 

−8 3 16 3 −8 3



u2 = 0.0229167;

u3 = 0.0333333

6. Conclusiones La solución análitica para el problema dado por su formulación fuerte puede obtenerse fácilmente. Las expresiones para los dos casos de carga considerados son: Carga constante Solución analítica. p 0 x2 Lp0 u(x) = − + x; 2EA EA Solución MEF mediante dos elementos lineales. ( (1) u2 H 2 si 0 ≤ x ≤ L2 u(x) = (2) (2) u2 H 1 + u3 H 2 si L2 ≤ x ≤ L 10

G. Rus, E. Puertas 0.05 0.04 0.03 0.02 0.01

2

4

6

8

10

Figure 4: Solución para carga constante Solución MEF mediante un elemento cuadrático. u(x) = u2 H2 + u3 H3 Carga variable

Lp0 p 0 x3 + x; 6EAL 2EA Solución MEF mediante dos elementos lineales. ( (1) u2 H 2 si 0 ≤ x ≤ L2 u(x) = (2) (2) u2 H 1 + u3 H 2 si L2 ≤ x ≤ L u(x) = −

Solución MEF mediante un elemento cuadrático. u(x) = u2 H2 + u3 H3 Sustituyendo los valores dados para las variables y comparando con los resultados obtenidos para los problemas planteados mediante el método de los elementos finitos podemos observar que se obtiene la solución real para los nodos. Si bien la solución es la real en el caso de carga constante empleando elementos cuadráticos pero no es válido para el caso de carga variable, ya que la solución real es un polinomio de tercer grado. Bastaría con emplear un elemento cúbico para comprobar que se obtiene la solución real.

11

G. Rus, E. Puertas

0.03 0.025 0.02 0.015 0.01 0.005 2

4

6

8

10

Figure 5: Solución para carga variable

12

G. Rus, E. Puertas

Problema 2 Resuelva por el método de los elementos finitos el problema de una viga a axil de longitud L = 3m empotrada por el extremo izquierdo y con una carga R = 100kN en el extremo derecho, sin cargas distribuidas f v , cuyo módulo elástico es E = 210GP a y su área es de A = 0.01m2 en su mitad izquierda (0 < x < L/s) y de A = 0.02m2 en su mitad derecha (L/2 < x < L). Utilícese para resolverlo una discretización de dos (2) elementos de igual longitud, y funciones de forma lineales. Se pide: 1. Definir los grados de libertad. 2. Obtener la matriz de rigidez de toda la barra. 3. Obtener el alargamiento total de la barra. 4. Obtener la ley completa de desplazamientos y de tensiones para 0 < x < L. Dibujar ambas leyes. (Ejercicio evaluado de diciembre de 2005. 1 hora.)

13

G. Rus, E. Puertas

Solución 2 1. Discretización y grados de libertad

14

G. Rus, E. Puertas

15

G. Rus, E. Puertas

16

G. Rus, E. Puertas

Problema 3 Se define un elemento sólido 2D de tensión plana. El módulo elástico es E = 210GP a y el de Poisson ν = 0.3. Se utiliza una discretización de un (1) elemento rectangular de un nodo en cada esquina, cuyas funciones de forma son lineales. La geometría es un rectángulo de 0.4 × 0.3m (horizontal × vertical), con un origen de coordenadas globales en el lado inferior izquierdo y un espesor unitario. Las condiciones de contorno son de empotramiento del lado izquierdo y de una tracción uniforme de t = 2GP a en el lado derecho. Se pide: 1. Indicar gráficamente los grados de libertad de la estructura. 2. Obtener la matriz de rigidez del elemento. 3. Obtener el vector de cargas. 4. Obtener los desplazamientos de los nodos.

Elemento 1

Nota: Considérese que el tensor B que relaciona tensiones con deformaciones según εij = Bijnc unc es,

Bijnc

  B11n1 B11n2 = B22n1 B22n2  B12n1 B12n2

según

   H1n,1 ε11 ε22  =  0 1 ε12 H 2 1n,2

(Examen final del 9-II-2006. 1 hora.)

0 H2n,2 1 H 2 2n,1



µ n¶  u1n u2

17

G. Rus, E. Puertas

Solución 3

18

G. Rus, E. Puertas

19

G. Rus, E. Puertas

20

G. Rus, E. Puertas

21

G. Rus, E. Puertas

Problema 4 Se define un elemento sólido 2D en tensión plana. El módulo elástico es E = 210GP a y el de Poisson ν = 0.3. Se utiliza una discretización de dos (2) elementos rectangulares de un nodo en cada esquina, cuyas funciones de forma son lineales. La geometría es un rectángulo de 0.4 × 0.2m (horizontal × vertical), con un origen de coordenadas globales en la esquina inferior izquierda y un espesor unitario. Las condiciones de contorno son de empotramiento del lado izquierdo, deslizaderas horizontales en los lados superior e inferior, y una tracción normal uniforme de t = 2GP a en el lado derecho. Se pide: 1. Indicar gráficamente los grados de libertad de la estructura, distinguiendo cuáles son de desplazamientos y cuáles de fuerzas. 2. Obtener la matriz de rigidez de un elemento. 3. Obtener el vector de cargas. 4. Obtener los desplazamientos de los nodos. Nota: Defínanse las funciones de forma, y hágase la integración en coordenadas naturales. Considérese que el tensor B que relaciona tensiones con deformaciones según εij = Bijnc unc es,

Bijnc

  B11n1 B11n2 = B22n1 B22n2  B12n1 B12n2

según

Elemento 1

22

   H1n,1 ε11 ε22  =  0 1 ε12 H 2 1n,2

Elemento 2

0 H2n,2 1 H 2 2n,1



µ n¶  u1n u2

G. Rus, E. Puertas

Solución 4

23

G. Rus, E. Puertas

24

G. Rus, E. Puertas

25

G. Rus, E. Puertas

Problema 5 Se define un muro de contención de hormigón (E=20 GPa, ν=0.3) de 20 metros de altura, 8 de base y 4 de coronación, cuya geometría y condiciones de contorno se describen en el dibujo. El muro se carga con una presión uniforme en el lado izquierdo de valor p=0.098 MPa. Considerando una sección en tensión plana de espesor unitario, y una discretización de un elemento cuadrado lineal, y un origen de coordenadas globales en la esquina inferior izquierda, se pide: 1. Indicar gráficamente los grados de libertad de la estructura, distinguiendo cuáles son de desplazamientos y cuáles de fuerzas. 2. Obtener las derivadas de las funciones de forma. 3. Obtener la matriz de rigidez de un elemento. 4. Obtener el vector de cargas. 5. Obtener los desplazamientos de los nodos, y esbozar gráficamente la deformada.

P

26

G. Rus, E. Puertas

Solución 5

27

G. Rus, E. Puertas

28

G. Rus, E. Puertas

29

G. Rus, E. Puertas

Problema 6 Se define un muro de contención de hormigón (E=20 GPa, ν=0.3) de 20 metros de altura, 8 de base y 4 de coronación, cuya geometría y condiciones de contorno se describen en el dibujo. El muro se carga con una presión hidrostática procedente del peso del agua que baña la cara izquierda, de modo que en la coronación tiene valor 0 y en la base valor 2p=0.196 MPa. Considerando una sección en tensión plana de espesor unitario, y una discretización de un elemento cuadrado lineal, y un origen de coordenadas globales en la esquina inferior izquierda, se pide: 1. Indicar gráficamente los grados de libertad de la estructura, distinguiendo cuáles son de desplazamientos y cuáles de fuerzas. 2. Obtener las derivadas de las funciones de forma. 3. Obtener la matriz de rigidez de un elemento. 4. Obtener el vector de cargas. 5. Obtener los desplazamientos de los nodos, y esbozar gráficamente la deformada.

2P

30

G. Rus, E. Puertas

Solución 6

31

G. Rus, E. Puertas

Problema 7 Se define un muro de contención de hormigón (E=20 GPa, ν=0.3) de 20 metros de altura, 8 de base y 4 de coronación, cuya geometría y condiciones de contorno se describen en el dibujo. El muro se carga con una presión hidrostática procedente del peso del agua que baña la cara izquierda, de modo que en la coronación tiene valor 0 y en la base valor 2p=0.196 MPa. Considerando una sección en tensión plana de espesor unitario, y una discretización de dos elementos cuadrados lineales, tal y como se describe en el dibujo, y un origen de coordenadas globales en la esquina inferior izquierda, se pide: 1. Indicar gráficamente los grados de libertad de la estructura, distinguiendo cuáles son de desplazamientos y cuáles de fuerzas. 2. Obtener las derivadas de las funciones de forma. 3. Obtener la matriz de rigidez de un elemento. 4. Obtener el vector de cargas. 5. Obtener los desplazamientos de los nodos, y esbozar gráficamente la deformada.

e=1

e=2 2P

32

G. Rus, E. Puertas

Solución 7

33

G. Rus, E. Puertas

34

G. Rus, E. Puertas

35

G. Rus, E. Puertas

Problema 8 Se define un muro de contención de hormigón (E=20 GPa, ν=0.3) de 20 metros de altura, 8 de base y 4 de coronación, cuya geometría y condiciones de contorno se describen en el dibujo. El muro tiene como única carga la gravitatoria, debida al peso propio originado por la densidad del hormigón ρ=2700 kg/m3 . Considerando una sección en tensión plana de espesor unitario, y una discretización de un elemento cuadrado lineal, y un origen de coordenadas globales en la esquina inferior izquierda, se pide: 1. Indicar gráficamente los grados de libertad de la estructura, distinguiendo cuáles son de desplazamientos y cuáles de fuerzas. 2. Obtener las derivadas de las funciones de forma. 3. Obtener la matriz de rigidez de un elemento. 4. Obtener el vector de cargas. 5. Obtener los desplazamientos de los nodos, y esbozar gráficamente la deformada.

g

36

G. Rus, E. Puertas

Solución 8

37

G. Rus, E. Puertas

Problema 9 Se aisla un bloque cuadrado de un estrato de terreno (E=2 GPa, ν=0.4) de 20 metros de lado, cuya geometría y condiciones de contorno se describen en el dibujo. El terreno tiene como única carga la gravitatoria, debida al peso propio originado por su densidad ρ=2300 kg/m3 . Considerando una sección en tensión plana de espesor unitario, y una discretización de un elemento cuadrado lineal, y un origen de coordenadas globales en la esquina inferior derecha, se pide: 1. Indicar gráficamente los grados de libertad de la estructura, distinguiendo cuáles son de desplazamientos y cuáles de fuerzas. 2. Obtener las derivadas de las funciones de forma. 3. Obtener la matriz de rigidez de un elemento. 4. Obtener el vector de cargas. 5. Obtener los desplazamientos de los nodos, y esbozar gráficamente la deformada.

g

38

G. Rus, E. Puertas

Solución 9

39

G. Rus, E. Puertas

40

G. Rus, E. Puertas

41

G. Rus, E. Puertas

Problema 10 Se define un sistema compuesto por la interacción de un bloque cuadrado de un estrato de terreno (E=2 GPa, ν=0.4) de 20 metros de lado (izquierda), contenido por un muro de hormigón (E=20 GPa, ν=0.3) de 20 metros de altura, 8 de base y 4 de coronación (derecha), cuya geometría y condiciones de contorno se describen en el dibujo. El sistema tiene como única carga la gravitatoria, debida al peso propio originado por la densidad del hormigón ρ=2700 kg/m3 y del terreno ρ=2300 kg/m3 . Considerando una sección en tensión plana de espesor unitario, y una discretización de un elemento cuadrado lineal, se pide: 1. Indicar gráficamente los grados de libertad de la estructura, distinguiendo cuáles son de desplazamientos y cuáles de fuerzas. 2. Obtener las derivadas de las funciones de forma. 3. Obtener la matriz de rigidez de un elemento. 4. Obtener el vector de cargas. 5. Obtener los desplazamientos de los nodos, y esbozar gráficamente la deformada.

g

42

g

G. Rus, E. Puertas

Solución 10

43

G. Rus, E. Puertas

44

G. Rus, E. Puertas

Problema 11 Resuelva por el Método de los Elementos Finitos la estructura de la figura, cuyo módulo elástico es E = 2 1011 P a y el coeficiente de Poisson ν = 0.3. La estructura está empotrada en tres de sus cuatro lados, estando el cuarto (el derecho) sometido a una tracción uniforme p = 1GP a. Utilícese para su resolución elementos de dos nodos y cuadrados de cuatro nodos con funciones de forma lineales y tenga en cuenta que se considera en tensión plana.

Se pide: 1. Discretizar la estructura y definir los grados de libertad. 2. Obtener las funciones de forma. 3. Definir la matriz de rigidez global y calcular las matrices de rigidez elementales necesarias para el cálculo de los desplazamientos de los nodos. 4. Definir el vector de cargas global y calcular los miembros necesarios para la obtención de los desplazamientos de los nodos. 5. Obtener el desplazamiento de los nodos.

45

G. Rus, E. Puertas

Solución 11 1. Discretización La estructura se discretiza siguiendo el siguiente esquema:

Figure 6: Discretización de la estructura: nodos globales Existen dos grados de libertad por nodo. En total serán 12 grados de libertad, de los cuales sólo los dos correspondientes al nodo 4 corresponden a desplazamientos.

2. Funciones de forma Para un elemento lineal, las funciones de forma en coordenadas locales son (dibujar las funciones de forma): H 1 = 1 − x′ H 2 = x′ Para un elemento cuadrado de cuatro nodos, las funciones de forma en coordenadas locales son (dibujar las funciones de forma): H1 H2 H3 H4 46

= = = =

1 (1 4 1 (1 4 1 (1 4 1 (1 4

+ x′1 )(1 + x′2 ) − x′1 )(1 + x′2 ) − x′1 )(1 − x′2 ) + x′1 )(1 − x′2 )

G. Rus, E. Puertas

Figure 7: Discretización de la estructura: nodos locales

3. Matriz de rigidez global Los elementos de las matrices de rigidez para un elemento plano 2D se obtienen a partir de la integración:

(e) kncmd

=

Z

V(e)

(e)

(e)

Cijkl Bijnc Bklmd dV

La matriz de rigidez global se obtiene mediante el proceso de ensamblaje:



    K=    

(1)

k11 (1) k21 (1) k31 (1) k41 0 0

(1)

(1)

(1)

k12 k13 k14 (1) (1) (1) k22 k23 k24 (1) (2) (3) (3) (1) (1) (2) k32 k33 + k22 + k11 k34 + k21 + k12 (1) (1) (2) (3) (1) (2) (3) k42 k43 + k12 + k21 k44 + k11 + k22 (2) (2) k41 0 k42 (2) (2) k31 0 k32

0 0 (2) k24 (2) k14 (2) k44 (2) k34

0 0 (2) k23 (2) k13 (2) k43 (2) k33

         

47

G. Rus, E. Puertas En consecuencia, el sistema a resolver será  (1) (1) (1) (1) k11 k12 k13 k14  (1) (1) (1) (1)  k21 k22 k23 k24  (1) (1) (1) (2) (3) (1) (2) (3)  k31 k32 k33 + k22 + k11 k34 + k21 + k12  (1) (1) (1) (2) (3) (1) (2) (3)  k  41 k42 k43 + k12 + k21 k44 + k11 + k22  (2) (2) k41 0 k42  0 (2) (2) 0 0 k32 k31

0 0 (2) k24 (2) k14 (2) k44 (2) k34

0 0 (2) k23 (2) k13 (2) k43 (2) k33



         

U1 U2 U3 U4 U5 U6





      =      

4. Desplazamientos de los nodos Como todos los nodos están fijos a excepción del nodo 4, eliminamos las filas y columnas correspondientes a estos nodos (1, 2, 3, 5, 6) y nos queda un sistema de dos ecuaciones con dos incógnitas. i i h ¤ h (1) (2) (3) (1) (2) (3) £ k44 + k11 + k22 U 4 = f4 + f1 + f2 Ã

(1)

(2)

(3)

k4141 + k1111 + k22 (1) (2) k4241 + k1211

(1)

(2)

k4142 + k1112 (1) (2) k4242 + k1212

!µ

u41 u42

¶

=

Ã

(1)

(2)

(3)

f41 + f11 + f21 (1) (2) (3) f42 + f12 + f22

!

La simetría de la estructura nos permite concluir que u42 = 0. Teniendo en cuenta esta condición eliminamos las fila y columna segunda y nos quedamos con una única ecuación con una incógnita: ´ ³ (1) (2) (3) (1) (2) (3) k4141 + k1111 + k22 u41 = f41 + f11 + f21 Elemento 1 (1) k4141

=

Z

(1)

V(1)

(1)

Cijkl Bij41 Bkl41 dV

siendo la matriz C la definida para tensión plana σ3 = 0:   1 ν 0 Ez   ν 1 0 C= 1 − ν2 0 0 2 (1 − ν) El jacobiano de la tansformación es: Ã dx dx ! J=

48

1

2

dx′1 dx1 dx′2

dx′1 dx2 dx′2

=

µ

2 0 0 1.5

¶

f1 f2 f3 f4 f5 f6

       

G. Rus, E. Puertas cuyo determinante es: |J| = 3. El tensor B que relaciona tensiones con deformaciones según εij = Bijnc unc es:       µ n¶ H1n,1 0 ε11 B11n1 B11n2 u1       0 H2n,2 ε22 = según Bijnc = B22n1 B22n2 un2 1 1 H H ε12 B12n1 B12n2 1n,2 2n,1 2 2 Por tanto,

Bij41 = Bkl41 Se tiene

    H14,1 B1141 = 0 =  B2241  =  1 B1241 H 2 14,2 

1 (1 8

 − x′2 )  0 1 ′ (1 − x2 ) 16

 1 ν 0 ¢ zE ¡ 1 (1) 1   ν 1 (1 − x′2 ) 0 16 (1 − x′2 ) 0 k4141 = 8 1 − ν2 −1 −1 0 0 2 (1 − ν)   1 (1 − x′2 ) 8  3 dx′1 dx′2  0 1 ′ (1 − x2 ) 16 Z

1

Z



1

(1)

k4141 = 63.8991 106 N/m

Elemento 2 El elemento 2 es simétrico con el 1, por lo que el valor será el mismo al anteriormente indicado. Elemento 3 Para el elemento lineal tenemos la integración: Z L EA (3) (3) (3) k2121 = EAH2,1 H2,1 dx1 = = 500 106 N/m L 0 Vector de cargas (1)

f41 = (2)

f11 =

R

R

(1)

S(1)

f1S H14 dS =

S(2)

f1S H11 dS =

(2)

R1

R−1 1 −1

dx2 dx′2 dx′2 2 H1 |x′ =1 dx dx′2 dx′2 1

p H4 |x′ =1

= 15 106 N

p

= 15 106 N

1

f14 = 30. 106 N 49

G. Rus, E. Puertas Sistema de ecuaciones Se tiene entonces el sistema: ¡ ¢ 2 · 63.8991 106 + 500 106 u41 = 30. 106 u41 = 47.786 mm

50

G. Rus, E. Puertas

Problema 12 Se define la sección de un azud de hormigón (E=20 GPa, ν=0.3) de 4 × 1 metros sometida a una presión constante de valor p = 44.145 103 N/m tal y como se indica en la figura. Considerando una sección en tensión plana de espesor unitario, y una discretización de dos elementos cuadrados lineales, tal y como se describe en el dibujo, y un origen de coordenadas globales en la esquina inferior izquierda, se pide: 1. Indicar gráficamente los grados de libertad de la estructura, distinguiendo cuáles son de desplazamientos y cuáles de fuerzas. 2. Hallar las funciones de forma y sus derivadas. 3. Describir la matriz de rigidez. 4. Obtener el vector de cargas. 5. Obtener los desplazamientos de los nodos.

Solución 12 Grados de libertad La discretización de la sección se realiza mediante dos elementos cuadrados lineales, tal y como se refleja en la figura 8 51

G. Rus, E. Puertas

(a)

(b)

Figure 8: Discretización de la estructura: (a) Nodos globales; (b) Nodos locales Funciones de forma Las funciones de forma para un elemento cuadrado lineal en coordenadas locales son: H1 = 14 (1 + x′1 )(1 + x′2 ) H2 = 14 (1 − x′1 )(1 + x′2 ) H3 = 14 (1 − x′1 )(1 − x′2 ) H4 = 14 (1 + x′1 )(1 − x′2 )

Las derivadas de las funciones de forma respecto a las coordenadas locales son: à ∂H ! µ ¶ 1 1 (1 + x′2 ) ∂x′1 4 = 1 ∂H1 (1 + x′1 ) ′ ∂x à ∂H2 ! µ 4 ¶ 2 − 41 (1 + x′2 ) ∂x′1 = 1 ∂H2 (1 − x′1 ) ′ ∂x à ∂H2 ! µ 4 ¶ 3 − 41 (1 − x′2 ) ∂x′1 = ∂H3 − 1 (1 − x′1 ) ′ ∂x 2 à ∂H ! µ 4 ¶ 4 1 (1 − x′2 ) ∂x′1 4 = ∂H4 − 14 (1 + x′1 ) ∂x′ 2

El jacobiano de la transformación se obtiene teniendo en cuenta la relación entre las coordenadas locales y globales: xi = Hn xni 52

G. Rus, E. Puertas ∂Hn n ∂xi = x = J= ′ ∂xj ∂x′j i

Ã

∂x1 ∂x′1 ∂x1 ∂x′2

∂x2 ∂x′1 ∂x2 ∂x′2

!

Cada término de la matriz se obtiene: ∂x1 ∂x′1 ∂x1 ∂x′2 ∂x2 ∂x′1 ∂x2 ∂x′2

= = = =

∂H1 ∂x′1 ∂H1 ∂x′2 ∂H1 ∂x′1 ∂H1 ∂x′2

x11 + x11 + x12 + x12 +

∂H2 ∂x′1 ∂H2 ∂x′2 ∂H2 ∂x′1 ∂H2 ∂x′2

x21 + x21 + x22 + x22 +

∂H3 ∂x′1 ∂H3 ∂x′2 ∂H3 ∂x′1 ∂H3 ∂x′2

x31 + x31 + x32 + x32 +

∂H4 ∂x′1 ∂H4 ∂x′2 ∂H4 ∂x′1 ∂H4 ∂x′2

x41 x41 x42 x42

Para el primer elemento se tiene: 1 1 1 1 ∂x1 = 4 (1 + x′2 ) − 2 (1 + x′2 ) − 2 (1 − x′2 ) + 4 (1 − x′2 ) = 1 ′ ∂x1 4 4 4 4 El jacobiano para ambos elementos es igual debido a que son de la misma forma: µ ¶ 1 0 (1) (2) J =J = 0 21 y su determinante |J| = ¡

J

1 2

La inversa del jacobiano es inmediata:

¢ (1) −1

¢−1 ¡ = J (2) =

µ

1 0 0 2

¶

Las derivadas de las funciones de forma en coordenadas globales se obtienen al multiplicar la derivada en coordenadas locales por la inversa del jacobiano de la transformación. Éstas serán las mismas para ambos elementos. ¶ ¶ Ã ∂x′1 ∂x′2 ! Ã ∂H1 ! µ 1 µ ′ H1,1 (1 + x ) ∂x′1 2 ∂x1 ∂x1 4 = H1,i = = 1 ∂H1 ∂x′1 ∂x′2 (1 + x′1 ) H1,2 2 ∂x′ ∂x2

H2,i =

µ

H2,1 H2,2

¶

H3,i =

µ

H3,1 H3,2

¶

H4,i =

µ

H4,1 H4,2

¶

∂x2

2

=

Ã

∂x′1 ∂x1 ∂x′1 ∂x2

∂x′2 ∂x1 ∂x′2 ∂x2

!Ã

∂H2 ∂x′1 ∂H2 ∂x′2

!

=

Ã

∂x′1 ∂x1 ∂x′1 ∂x2

∂x′2 ∂x1 ∂x′2 ∂x2

!Ã

∂H3 ∂x′1 ∂H3 ∂x′2

!

=

Ã

∂x′1 ∂x1 ∂x′1 ∂x2

∂x′2 ∂x1 ∂x′2 ∂x2

!Ã

∂H4 ∂x′1 ∂H4 ∂x′2

!

=

µ

− 14 (1 + x′2 ) 1 (1 − x′1 ) 2

¶

=

µ

− 41 (1 − x′2 ) − 12 (1 − x′1 )

¶

=

µ

1 4

(1 − x′2 ) 1 − 2 (1 + x′1 )

¶ 53

G. Rus, E. Puertas Matriz de rigidez global El montaje de la matriz de rigidez:  (1) (1) k k1•2• 0  1•1• (1) (1) (2) (2)  k2•1• k2•2• + k1•1• k1•2•  (2) (2)  0 k2•1• k2•2• K= (2) (2)  0 k3•1• k3•2•   (1) (1) (2) (2)  k3•1• k3•2• + k4•1• k4•2• (1) (1) k4•1• k4•2• 0

(1)

0

(1)

k1•3• k1•4• (1) (2) (1) k2•3• + k1•4• k2•4• (2) k4•4• 0 (2) k3•4• 0 (1) (2) k3•3• + k4•4• 0 (1) 0 k4•4•

(2)

k1•3• (2) k2•3• (2) k3•3• (2) k4•3• 0

Por tanto, el sistema a resolver será:



(1) k1•1• (1) k2•1•

(1) k1•2•

0

 (1) (2)  k2•2• + k1•1•  (2)  0 k2•1•  (2)  0 k3•1•   (1) (1) (2)  k3•1• k3•2• + k4•1• (1) (1) k4•2• k4•1•

(2)

k1•2• (2) k2•2• (2) k3•2• (2) k4•2• 0

0 (2)

(1)

(1) k1•3•

(2)

(1) k1•4• (1) k2•4•

k1•3• k2•3• + k1•4• (2) (2) k2•3• k4•4• 0 (2) (2) k3•3• k3•4• 0 (2) (1) (2) k4•3• k3•3• + k4•4• 0 (1) 0 0 k4•4•





        

U1 U2 U3 U4 U5 U6



         



      =      

f1 f2 f3 f4 f5 f6

       

Las condiciones de contorno del problema indican que los desplazamientos en ambas direcciones son nulos para los nodos 1, 3, 4 y 6. De ahí que el sistema se reduzca: !µ Ã ¶ µ 2 ¶ (1) (2) (1) (2) f U2 k2•2• + k1•1• k2•3• + k1•4• = 5 (1) (2) (1) (2) f5 U k3•2• + k4•1• k3•3• + k4•4• El sistema completo será:     

(1)

(2)

k2121 + k1111 (1) (2) k2221 + k1211 (1) (2) k3121 + k4111 (1) (2) k3221 + k4211

(1)

(2)

k2122 + k1112 (1) (2) k2222 + k1212 (1) (2) k3122 + k4112 (1) (2) k3222 + k4212

(1)

(2)

k2131 + k1141 (1) (2) k2231 + k1241 (1) (2) k3131 + k4141 (1) (2) k3231 + k4241

(1)

(2)

k2132 + k1142 (1) (2) k2232 + k1242 (1) (2) k3132 + k4142 (1) (2) k3232 + k4242



  2 u21 f   2   12 u   2   f2  5  =  5 f1  u1 5 u2 f25

Por simetría u21 y u51 son cero. Con lo que el sistema se reduce a dos ecuaciones con dos incógnitas. !µ Ã ¶ µ 2 ¶ (1) (2) (1) (2) u22 k2222 + k1212 k2232 + k1242 f2 = 5 (1) (2) (1) (2) u2 f25 k3222 + k4212 k3232 + k4242 54

   

G. Rus, E. Puertas Cada una de las componentes a calcular de la matriz de rigidez se obtienen: Z (e) (e) (e) kncmd = Cijkl Bijnc Bklmd dV V(e)

siendo la matriz C la definida para tensión plana σ3 = 0:   1 ν 0 Ez   ν 1 0 C= 1 − ν2 0 0 2 (1 − ν)

Y el tensor B que relaciona tensiones con deformaciones según εij = Bijnc unc es:       µ n¶ H1n,1 0 ε11 B11n1 B11n2 u1       0 H2n,2 ε22 = según Bijnc = B22n1 B22n2 un2 1 1 ε12 B12n1 B12n2 H H 1n,2 2n,1 2 2 Por tanto, para la segunda dirección coordenada se tendrá:       0 B11n2 0 B =  B22n2  =  H2n,2  =  12 (1 − x′1 )  1 B12n2 − 18 (1 + x′2 ) H 2 2n,1

A continuación describimos el cálculo de uno de los ocho términos necesarios de la matriz de rigidez. Z Z 1Z 1 (1) (1) (1) (1) (1) k2222 = Cijkl Bij22 Bkl22 |J| dx′1 dx′2 Cijkl Bij22 Bkl22 dV = V(1)

−1

−1

  0 1 ν 0 ¢ Ez ¡ (1)   H2,2  1 dx′1 dx′2  ν 1 0 H2,2 12 H2,1 0 k2222 = 1 − ν2 2 1 −1 −1 H2,1 0 0 2 (1 − ν) 2 µ ¶ Z 1Z 1 1 1 ′ ′ Ez (1) H2,2 H2,2 + (1 − ν) H2,1 H2,1 dx dx k2222 = 2 2 2 1 2 −1 −1 1 − ν Z

1

Z

1



Vector de cargas Únicamente hay que calcular el vector de cargas en el nodo 2, ya que en el nodo 5 no hay ninguna carga actuando. Para el cálculo hay que tener en cuenta la carga correspondiente a cada elemento: (1) (2) f22 = f22 + f12 Detallamos el cálculo para el primer elemento: Z 1 Z dx1 ′ (1) (1) S f2 H22 dS = p H2 |x′ =1 f22 = dx 2 dx′1 1 −1 S (1) 55

G. Rus, E. Puertas Desplazamiento de los nodos El sistema de ecuaciones resultante es: ¶ ¶µ 2 ¶ µ µ u2 88290. 3.18681 101 0 −2.8022 101 0 = K= u52 0 −2.8022 101 0 3.18681 101 0 Resolviendo se obtiene que el desplazamiento en los nodos 2 y 5 es: µ 2 ¶ µ ¶ u2 1.22148 10−5 = u52 1.07406 10−5

56

G. Rus, E. Puertas

Problema 13 Se define un pequeño azud de hormigón (E=20 GPa, ν=0.3) de 3 metros de altura y sección 1 × 4m , cuya geometría y condiciones de contorno se describen en el dibujo. El azud está sometido a una presión hidrostática procedente del peso del agua, de modo que en la coronación tiene valor 0 y en la base valor p=29.43 kPa.

Considerando una discretización de dos elementos hexaédricos lineales de ocho nodos, tal y como se describe en el dibujo, se pide: 1. Indicar gráficamente los grados de libertad de la estructura, distinguiendo cuáles son de desplazamientos y cuáles de fuerzas. 2. Obtener las funciones de forma. 3. Definir la matriz de rigidez. 4. Obtener el vector de cargas. 5. Obtener los desplazamientos de los nodos.

Solución 13 Grados de libertad La discretización de la estructura se toma tal y como indica la siguiente figura (figura 9) 57

G. Rus, E. Puertas

(a)

(b)

Figure 9: Discretización de la estructura: (a) Nodos globales; (b) Nodos locales La estructura es tridimensional, por lo que hay que considerar tres grados de libertad por nodo, que en total serán 36 grados de libertad. Teniendo en cuenta las condiciones de contorno y la simetría de la estructura, únicamente se definen dos grados de libertad en desplazamientos, correspondientes con los nodos globales 9 y 12. Funciones de forma El elemento hexaédrico recto más sencillo de clase C 0 es el de ocho nodos que se muestra en la figura siguiente (figura 10). Las coordenadas locales de cada uno de los nodos se definen en la siguiente tabla: Nodo 1 2 3 4 5 6 7 8 58

x′1 -1 1 1 -1 -1 1 1 -1

x′2 -1 -1 1 1 -1 -1 1 1

x′3 -1 -1 -1 -1 1 1 1 1

G. Rus, E. Puertas

Figure 10: Elemento hexaédrico de 8 nodos Las funciones de forma de un nodo se obtienen, como producto de las tres funciones de una sola variable correspondientes a cada una de las tres direcciones x′1 , x′2 , x′3 en ese nodo. La expresión general de la función de forma de un nodo cualquiera i, se obtiene mediante la siguiente expresión: Hn(e) =

1 ′ ′n ′ ′n ′ (1 + x′n 1 x1 ) (1 + x2 x2 ) (1 + x3 x3 ) 8

siendo x′n j la coordenada j-ésima del nodo n. Para el nodo 5, la función de forma se puede obtener fácilmente tal y como se indica en la figura 11. (e)

H5 =

1 1 1 (1 − x′1 ) (1 − x′2 ) (1 + x′3 ) 2 2 2

Así, las funciones de forma para un elemento hexaédrico lineal de ocho nodos son: H1 = 18 (1 − x′1 ) (1 − x′2 ) (1 − x′3 ) H2 = 81 (1 + x′1 ) (1 − x′2 ) (1 − x′3 ) H3 = 81 (1 + x′1 ) (1 + x′2 ) (1 − x′3 ) H4 = 81 (1 − x′1 ) (1 + x′2 ) (1 − x′3 ) H5 = 81 (1 − x′1 ) (1 − x′2 ) (1 + x′3 ) H6 = 81 (1 + x′1 ) (1 − x′2 ) (1 + x′3 ) H7 = 81 (1 + x′1 ) (1 + x′2 ) (1 + x′3 ) H8 = 81 (1 − x′1 ) (1 + x′2 ) (1 + x′3 ) 59

G. Rus, E. Puertas

Figure 11: Obtención de la función de forma para el nodo 5 Las derivadas de las funciones de forma respecto a las coordenadas locales:  dH1     dH5   1  − 18 (1 − x′2 ) (1 − x′3 ) − 8 (1 − x′2 ) (1 + x′3 ) dx′1 dx′1  5   1  dH1   1 − 8 (1 − x′1 ) (1 − x′3 )   dH − 8 (1 − x′1 ) (1 + x′3 )   dx′2  = dx′2  = 1 dH1 dH5 − 18 (1 − x′1 ) (1 − x′2 ) (1 − x′1 ) (1 − x′2 ) 8 dx′3  dx′3    dH  1  1   dH6 2 ′ ′ ′ ′ ′ (1 − x ) (1 − x ) (1 − x ) (1 + x ) dx′1 dx 2 3 2 3 1  6   81  dH2   81 − 8 (1 + x′1 ) (1 − x′3 )   dH − 8 (1 + x′1 ) (1 + x′3 )   dx′2  = dx′2  = 1 dH2 dH6 − 18 (1 + x′1 ) (1 − x′2 ) (1 + x′1 ) (1 − x′2 ) 8 dx′3  dx′3    dH  1  1   dH7 3 ′ ′ ′ ′ ′ (1 + x ) (1 − x ) (1 + x ) (1 + x ) dx′1 dx 2 3 2 3 1  7   81  dH3   81 (1 + x′1 ) (1 − x′3 )   dH (1 + x′1 ) (1 + x′3 )   dx′2  = dx′2  = 8 8 1 dH3 dH7 − 81 (1 + x′1 ) (1 + x′2 ) (1 + x′1 ) (1 + x′2 ) 8 dx′3  dx′3    dH    1  1 dH8 4 − 8 (1 + x′2 ) (1 − x′3 ) − 8 (1 + x′2 ) (1 + x′3 ) dx′1 dx′1  8   1  dH4   1 (1 − x′1 ) (1 − x′3 )   dH (1 − x′1 ) (1 + x′3 )   dx′2  = dx′2  = 8 8 1 1 ′ ′ dH4 dH8 (1 − x′1 ) (1 + x′2 ) − 8 (1 − x1 ) (1 + x2 ) 8 dx′ dx′ 3

3

El jacobiano de la transformación se halla teniendo en cuenta la relación entre las coordenadas globales y locales:  ∂x1 ∂x2 ∂x3   J =

∂x′1 ∂x1 ∂x′2 ∂x1 ∂x′3

∂x′1 ∂x2 ∂x′2 ∂x2 ∂x′3

∂x′1 ∂x3 ∂x′2 ∂x3 ∂x′3

 

A modo de ejemplo calculamos un elemento de la matriz: ∂H1 1 ∂H2 2 ∂H3 3 ∂H4 4 ∂H5 5 ∂H6 6 ∂H7 7 ∂H8 8 1 ∂x1 = x+ x+ x+ x+ x+ x+ x+ x = ′ ∂x1 ∂x′1 1 ∂x′1 1 ∂x′1 1 ∂x′1 1 ∂x′1 1 ∂x′1 1 ∂x′1 1 ∂x′1 1 2 60

G. Rus, E. Puertas Obteniendo que el jacobiano es igual para ambos elementos:   1 0 0 2 J (1) = J (2) =  0 1 0  0 0 32

y su determinante es |J| = 34 . La inversa se obtiene:   2 0 0 ¡ (1) ¢−1 ¡ (2) ¢−1 J = J = 0 1 0  0 0 32

Una vez calculado el jacobiano de la transformación, las derivadas en coordenadas globales de las funciones de forma se obtienen multiplicando éste por las derivadas en coordenadas locales.   ∂x′1 ∂x′2 ∂x′3   ∂Hn     Hn,1 2 Hn,1 ∂x′1 ∂x1 ∂x1 ∂x1  ′1 ∂x′2 ∂x′3   ∂Hn   Hn,2  Hn,i =  Hn,2  =  ∂x = ′ ∂x2 ∂x2 ∂x2   ∂x2  ′ ′ ′ 2 ∂H ∂x1 ∂x2 ∂x3 n Hn,3 Hn,3 3 ∂x′ ∂x3

∂x3

∂x3

3

Matriz de rigidez global La matriz de rigidez, teniendo en cuenta que es simétrica y sólo se representan los términos de la diagonal y por encima de ésta:  (1) (1) (1) (1) (1) (1) (1) k1•1• k1•2• k1•3• 0 0 k1•4• k1•5• k1•6• k1•7• 0 (1) (1) (1) (1) (1) (1)  k2•2• k2•3• 0 0 k2•4• k2•5• k2•6• k2•7• 0   (1) (1) (1) k3•7• k3•4• k3•3•  (2) (1) (1) (2) (2) k2•7• k3•5• k3•6• k2•3• k2•4•  (2) (2) (2)  +k2•6• +k2•1• +k2•2•  (2) (2) (2) (2) (2)  k3•3• k3•4• k3•1• 0 0 k3•6• k3•7•  (2) (2) (2) (2)  k4•4• k4•1• 0 0 k4•6• k4•7•   (1) (1)  k4•7• k4•4• (2) (1) (1)  k1•7• k4•5• k4•6• (2) (2)  +k1•6• +k1•1•  (1) (1) (1)  k5•5• k5•6• k5•7• 0   (1) (1) k6•6• k6•7• 0   (1)  k7•7• (2)  k6•7• (2)  +k 6•6•  (2)  k7•7•      61

0 0 (2)

k2•8• (2)

k3•8• (2) k4•8• (2)

k1•8• 0 0 (8)

k6•8• (2)

k7•8• (2) k8•8•

(

k1 ( k2 ( k3 +k ( k3 ( k4 ( k4 +k ( k5 ( k6 ( k7 +k ( k7 ( k8 ( k8 +k

G. Rus, E. Puertas Como el desplazamiento en los nodos 1, 2, 3,4, 5, 6,7, 8, 10 y 11 es nulo. El sistema se simplifica: !µ Ã ¶ µ 9 ¶ (1) (2) (1) (2) f U9 k7•7• + k6•6• k7•8• + k6•5• = 12 (1) (2) f 12 U k8•8• + k5•5• El sistema completo será entonces, teniendo en cuenta que la matriz de rigidez es simétrica:  (1) (2) (1) (2) (1) (2) (1) (2) (1) (2) (1) (2) k7171 + k6161 k7172 + k6162 k7173 + k6163 k7181 + k6151 k7182 + k6152 k7183 + k6153  (1) (2) (1) (2) (1) (2) (1) (2) (1) (2)  k7272 + k6262 k7273 + k6263 k7182 + k6152 k7282 + k6252 k7283 + k6253  (1) (2) (1) (2) (1) (2) (1) (2)  k7373 + k6363 k7381 + k6351 k7382 + k6352 k7383 + k6353  (1) (2) (1) (2) (1) (2)  k8181 + k5151 k8182 + k5152 k8183 + k5153   (1) (2) (1) (2) k8282 + k5252 k8283 + k5253  (2) (1) k8383 + k5353

Observando la simetría del problema respecto al plano x1 x3 podemos concluir que los desplazamientos en ambos nudos en la dirección x2 son nulos. De ahí que el sistema se reduzca a cuatro ecuaciones con cuatro incógnitas: 

(1)

(2)

(1)

(2)

(1)

(2)

(1)

(2)

k7171 + k6161 k7173 + k6163 k7181 + k6151 k7183 + k6153  (1) (2) (2) (1) (2) (1) k7373 + k6363 k7381 + k6351 k7383 + k6353   (1) (2) (1) (2)  k8181 + k5151 k8183 + k5153 (1) (2) k8383 + k5353 Cada componente de la matriz de rigidez se calcula: Z (e) (e) (e) kncmd = Cijkl Bijnc Bklmd dV



  9 u91 f   9   19   u3   f3   12  =  12 f1  u1 12 u3 f312

V(e)

siendo la matriz C la relación entre tensión y deformación σ  λ + 2µ λ λ 0 0 0  λ λ + 2µ λ 0 0 0   λ λ λ + 2 µ 0 0 0 C=  0 0 0 2µ 0 0   0 0 0 0 2µ 0 0 0 0 0 0 2µ donde λ = 62

νE (1+ν)(1−2 ν)

yµ=

E 2 (1+ν)

= Cǫ:        

son las constantes de Lamé.

   

         

G. Rus, E. Puertas La matriz B de interpolación del campo de deformaciones ǫ = B un es:   Hn,1 0 0  0 Hn,2 0      0 0 H n,3  B=  1 Hn,2 1 Hn,1 0  2   12 1   Hn,3 0 H n,1 2 2 1 1 Hn,3 2 Hn,2 0 2

Presentamos el cálculo de uno de los elementos de la matriz de rigidez: Z 1Z 1Z 1 (1) (1) (1) k7171 = Cijkl Bij71 Bkl71 |J| dx′1 dx′2 dx′3 −1

−1

−1

¢ R1 R1 R1 ¡ (1) 1 1 H k = 7,1 0 0 2 H7,2 2 H7,3 0 7171 −1 −1 −1    λ + 2µ λ λ 0 0 0 H7,1    λ λ + 2µ λ 0 0 0    0     λ λ λ + 2µ 0 0 0     1 0  |J| dx′1 dx′2 dx′3    0 0 0 2µ 0 0   2 H7,2     0 0 0 0 2 µ 0   12 H7,3  0 0 0 0 0 0 2µ ¶ µ Z 1Z 1Z 1 1 1 (1) (λ + 2 µ) H7,1 H7,1 + µ H7,2 H7,2 + µ H7,3 H7,3 |J| dx′1 dx′2 dx′3 k7171 = 2 2 −1 −1 −1 Vector de cargas Únicamente existe fuerza aplicada sobre el nodo 12 y en la primera dirección coordenada, por lo que será necesario el cálculo de la carga en este nodo: (1)

(2)

f112 = f81 + f51

El primero de estos sumandos se obtiene: ¶ Z Z 1Z 1 µ x′3 dx2 dx3 ′ ′ (1) V p 1− f81 = f1 H8 dV = H8 |x′ =0 dx dx 1 3 dx′2 dx′3 2 3 −1 −1 V

63

G. Rus, E. Puertas

Problema 14 Se define un muro de contención de hormigón (E=20 GPa, ν=0.3) de 30 metros de altura, 12 de base y 4 de coronación, cuya geometría y condiciones de contorno se describen en el dibujo. El muro se carga con una presión hidrostática procedente del peso del agua que baña la cara izquierda, de modo que en la coronación tiene valor 0 y en la base valor 2p=0.3 MPa.

Considerando sección en tensión plana de espesor unitario, y una discretización de un elemento cuadrado lineal, y un origen de coordenadas globales en la esquina inferior izquierda, se pide: 1. Indicar gráficamente los grados de libertad de la estructura, distinguiendo cuáles son de desplazamientos y cuáles de fuerzas. 2. Obtener las derivadas de las funciones de forma. 3. Obtener la componente k1111 de la matriz de rigidez del elemento que lo compone, donde el nodo 1 es el superior derecho, y la dirección 1 es la horizontal. (Examen Final de febrero de 2007. 1h 30min.) 64

G. Rus, E. Puertas

Solución 14 Grados de libertad La discretización de la estructura se toma tal y como indica en el enunciado del problema (figura 9)

(a)

(b)

Figure 12: Discretización de la estructura: (a) Nodos globales; (b) Nodos locales La estructura es bidimensional, por lo que hay que considerar dos grados de libertad por nodo, que en total serán 8 grados de libertad. Teniendo en cuenta las condiciones de contorno y la simetría de la estructura, únicamente se definen grados de libertad en desplazamientos en los nodos superiores, correspondientes con los nodos globales 1 y 2. Para los nodos 3 y 4 los grados de libertad vienen definidos en tensiones. Derivadas de las funciones de forma El elemento más sencillo de clase C 0 es el de cuatro nodos que se muestra en la figura 12. Las coordenadas locales de cada uno de los nodos se definen en la siguiente tabla: 65

G. Rus, E. Puertas Nodo 1 2 3 4

x′1 1 -1 -1 1

x′2 1 1 -1 -1

Las funciones de forma de un nodo se obtienen, como producto de las dos funciones de una sola variable correspondientes a cada una de las dos direcciones x′1 , x′2 en ese nodo. H1 H2 H3 H4

= = = =

1 (1 4 1 (1 4 1 (1 4 1 (1 4

+ x′1 )(1 + x′2 ) − x′1 )(1 + x′2 ) − x′1 )(1 − x′2 ) + x′1 )(1 − x′2 )

Las derivadas de las funciones de forma respecto a las coordenadas locales son: Ã ∂H ! µ ¶ 1 1 ′ (1 + x ) ∂x′1 2 4 = 1 ∂H1 ′ (1 + x ′ 1) 4 ∂x2 ! Ã ∂H ¶ µ 2 − 41 (1 + x′2 ) ∂x′1 = 1 ∂H2 (1 − x′1 ) 4 ∂x′2 ! Ã ∂H ¶ µ 1 3 − 4 (1 − x′2 ) ∂x′1 = ∂H3 − 14 (1 − x′1 ) ∂x′2 ! Ã ∂H µ 1 ¶ 4 ′ (1 − x ) ∂x′1 2 4 = ∂H4 − 14 (1 + x′1 ) ∂x′ 2

El jacobano de la transformación se halla teniendo en cuenta la relación entre las coordenadas globales y locales: Ã ∂x ∂x ! J=

1

2

∂x′1 ∂x1 ∂x′2

∂x′1 ∂x2 ∂x′2

A modo de ejemplo calculamos un elemento de la matriz: ∂x1 ∂H1 1 ∂H2 2 ∂H3 3 ∂H4 4 = x + x + x + x = 2(2 − x′2 ) ′ ∂x1 ∂x′1 1 ∂x′1 1 ∂x′1 1 ∂x′1 1 Obteniendo que el jacobiano es: ¶ µ 2(2 − x′2 ) 0 J= −2(1 − x′1 ) 15 66

G. Rus, E. Puertas Su inverso: J −1 =

Ã

1 2(2−x′2 ) 1+x′1 15(2−x′2 )

0 1 15

!

y su determinante es |J| = 30(2 − x′2 ) Una vez calculado la inversa del jacobiano de la transformación, las derivadas en coordenadas globles de las funciones de forma se obtienen multiplicando éste por las derivadas en coordenadas locales. Ã ! ¶ µ ∂Hn ′ Hn,1 ∂x1 = J −1 Hn,i = ∂Hn Hn,2 ∂x′ 2

H1,i =

µ

H2,1 H2,2

¶

H3,i =

µ

H3,1 H3,2

¶

H4,i =

µ

H4,1 H4,2

¶

H2,i =

µ

H1,1 H1,2

1+x′2 8(2−x′2 ) 1+x′1 20(2−x′2 )



¶

=

 

 1+x′ − 8 2−x2′ =  1−3 (x′ −22x)′  

1

60(2−x′2 )

2

 1−x′ − 8 2−x2′ =  3−x(′1 −2 x2′2)  60(2−x′2 )   1−x′ 

2

=

8(2−x′2 ) 1+x′ − 60 2−x1 ′ ( 2)



Elemento k1111 de la matriz de rigidez Cada componente de la matriz de rigidez se calcula: Z (e) (e) (e) kncmd = Cijkl Bijnc Bklmd dV V(e)



1 ν Ez  ν 1 C= 1 − ν2 0 0

 0 0 

1−ν 2

Y el tensor B que relaciona tensiones con deformaciones según εij = Bijnc unc es:       µ n¶ H1n,1 0 ε11 B11n1 B11n2 u H2n,2  1n Bijnc = B22n1 B22n2  según ε22  =  0 u2 1 1 ε12 B12n1 B12n2 H H 2 1n,2 2 2n,1

67

G. Rus, E. Puertas Presentamos el cálculo de uno de los elementos de la matriz de rigidez: Z 1Z 1 (1) (1) Bij11 Cijkl Bkl11 |J| dx′1 dx′2 = k1111 = −1

−1



  0 H1,1 Ez 0   0  |J| dx′1 dx′2 = = −1 −1 H1,1 0 21 H1,2 1−ν 2 1−ν 1 H 2 2 1,2 ¶ µ Z 1Z 1 1 Ez H1,1 H1,1 + (1 − ν) H1,2 H1,2 |J| dx′1 dx′2 = = 2 8 −1 −1 1 − ν R1 R1 ¡

¢

1 ν  ν 1 0 0

= 3.93136 1010 N/m

68

G. Rus, E. Puertas

Problema 15 Se considera un muro de hormigón (E = 20GP a, ν = 0.3) sometido a una carga p = d7 104 P a , con las dimensiones representas en la figura. Empleando una discretización de dos elementos de forma cúbica y con funciones de forma lineales.

Se pide: 1. Obtener el desplazamiento en los nodos centrales del muro, A y B explicando y estructurando cada paso seguido en la resolución. 2. Resolver el problema anterior empleando el programa de elementos finitos FEAP. 3. Resolver el problema doblando el número de elementos en cada una de las direcciones x, y, z utilizando FEAP. 4. Comparar los resultados obtenidos. Nota: Los parámetros di coinciden con las cifras de su DNI.

69

G. Rus, E. Puertas

Problema 16 Resuelva por el Método de los Elementos Finitos la presa de la figura (cotas en metros), cuyo módulo elástico es E = 3.10d3d4 1010 P a, el coeficiente de Poisson ν = 0.5 y la densidad ρ = 26d5d6 kg/m3 . La estructura está empotrada en su base, estando el lado izquierdo sometido a una presión hidrostática, tal y como indica la figura. Utilícese para su resolución elementos cuadrados de cuatro nodos con funciones de forma lineales y tenga en cuenta que se considera en tensión plana.

Se pide: 1. Calcular los desplazamientos de los nodos considerando que la presa está únicamente sometida a presión hidrostática. Explique y estructure cada paso seguido en la resolución. 2. Resolver el problema anterior utilizando FEAP. 3. Resolver el problema considerando también el peso propio de la presa mediante FEAP. 4. Analizar los resultados obtenidos. Nota: Los parámetros di coinciden con las cifras de su DNI.

70

G. Rus, E. Puertas

Problema 17 Se considera una placa cuadrada (E = 206GP a, ν = 0.3) sometida a una carga puntual p = d7 106 N , con las dimensiones representadas en la figura. La placa se discretiza mediante cuatro elementos de cuatro nodos con funciones de forma lineal.

Se pide: 1. Calcular el desplazamiento en el centro de la placa considerando la simetría del problema, explicando y estructurando cada paso seguido en la resolución. 2. Resolver el problema anterior empleando el programa de elementos finitos FEAP. 3. Resolver el problema aumentando el número de elementos en cada una de las direcciones x, y utilizando FEAP, para los casos siguientes: (a) Doble número de elementos en la dirección x (b) Doble número de elementos en la dirección y (c) Doble número de elementos en ambas direcciones x, y 4. Comparar los resultados obtenidos. Nota: Los parámetros di coinciden con las cifras de su DNI.

71

G. Rus, E. Puertas

Problema 18 Preguntas de teoría. Se evaluará la capacidad de seleccionar y sintetizar la información relevante. 1. Describir y expresar matemáticamente la diferencia entre un problema de tensión plana y uno de deformación plana. Poner un ejemplo de cada problema. 2. Dada una barra de un elemento lineal (dos nudos) de 3 unidades de longitud, calcular el vector de cargas f debido a una carga distribuida triangularmente, con valor nulo a la izquierda y valor 5 unidades a la derecha. 3. Coméntense ventajas y limitaciones de la idea de hacer los elementos isoparamétricos. 4. Calcular el número de nudos que tendrá un elemento sólido 2D de forma cuadrangular tal que en la dirección x1 sea cuadrático y en la dirección x2 sea lineal. Escribir las funciones de forma asociadas a cada nodo. 5. Describir el sentido geométrico del jacobiano, y dónde se usa. 6. Describir en qué ocasiones es mejor elegir una discretización 3D con elementos cúbicos y en qué ocasiones con elementos tetraédricos. 7. Calcular el número de nodos que tiene un elemento cúbico 3D tal que tenga todas las funciones de forma polinómicas de hasta orden cúbico. 8. Sea un elemento cuadrado 2D cuadrático de 9 nodos, cuyas funciones de forma incluyen los términos (1, x, x2 , y, xy, x2 y, y2 , xy2 , x2 y2 ), al que se le elimina un nodo de modo que de sus funciones de forma el término (x), quedando 8. Explíquese si es posible, y el motivo, representar una deformación ε constante según la dirección x. 9. Enumerar las ventajas e inconvenientes de elementos triangulares frente a cuadrados en 2D. 10. Calcular el área de un elemento cuadrado cuyos nodos tienen coordenadas (4,3) (2,3) (1,1) (5,2) mediante integración numérica, usando un total de cuatro puntos de Gauss. 72

G. Rus, E. Puertas 11. Calcular el momento de inercia respecto al eje y=0 y respecto al eje y=centro de gravedad del elemento, de un elemento cuadrado cuyos nodos tienen coordenadas (4,3) (2,3) (1,1) (5,2) mediante integración numérica, usando un total de cuatro puntos de Gauss. 12. Calcular el área de un elemento cúbico cuyos nodos tienen coordenadas (4,3,1) (2,3,1) (1,1,1) (5,2,1) (4,3,4) (2,3,4) (1,1,5) (5,2,5) mediante integración numérica, usando un total de ocho puntos de Gauss. 13. Explíquese porqué en el problema 7, aunque todas las fuerzas aplicadas verticales son nulas, existen desplazamientos verticales. 14. Relacionese el resultado del problema 12 con el fenómeno de los baches que se producen al entrar en un puente al poco tiempo de construirlos. 15. Indique qué tipología de elementos utilizaría para calcular una presa de contrafuertes y describa sucintamente éstos. 16. Describir cómo se calculan las tensiones mediante el método de los elementos finitos.

73