EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

EL METODO DE LOS ELEMENTOS FINITOS EN EL ANALISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. [email protected]

Centro de Investigación Científica Escuela Politécnica del Ejército

RESUMEN Se presenta una introducción a la utilización de los elementos finitos en el análisis estructural de placas, para lo que se revisan las ecuaciones básicas de la Teoría de la Elasticidad, y se da una interpretación física al Método de Ensamblaje Directo, empleado en el Análisis Matricial de Estructuras; esa misma interpretación se puede usar también en elementos finitos. Además se establece una metodología genérica y directa tipo Rayleigh-Ritz, basada en la minimización de la energía potencial y en campos de desplazamientos predefinidos, para la deducción de las matrices de rigideces de los elementos finitos, tomando como referencia a un cuadrilátero plano para placas. Se adjunta un programa de computación para el análisis de placas planas que incluyen elementos finitos cuadriláteros, y se comparan resultados numéricos con paquetes comerciales de análisis estructural. Se presenta una gran variedad de ejercicios de aplicación, incluyendo condiciones de borde especiales.

ABSTRACT An introduction to the use of finite elements in structural analysis of plates is presented. Basic equations of Elasticity Theory are reviewed, and a physical interpretation of the Direct Assembly Method is introduced; such interpretation can be adapted to finite elements. A generic methodology of a Rayleigh-Ritz type, based on potential energy minimization and predefined displacements is established to deduce stiffness matrices for finite elements, using as a reference a flat plate quadrilateral. A computer program for flat plate analysis is provided, which includes quadrilateral finite elements. Numerical results of the program are compared to results in commercial structural analysis software. A large number of examples are presented, including special border conditions.

1.

INTRODUCCIÓN:

El método de los elementos finitos es un método genérico para obtener soluciones numéricas, con una precisión aceptable, a muchos problemas complejos de ingeniería, constituidos o modelados mediante continuos. A través del método de los elementos finitos se ha conseguido abordar, con eficiencia, problemas tan disímiles como el análisis estructural, la transferencia de calor, el flujo de fluidos, los campos eléctricos, etc.

1

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Si se quisieran determinar los desplazamientos en la placa de la figura, los métodos clásicos nos conducirían al planteamiento de ecuaciones diferenciales parciales sin solución numérica específica, debido a que la estructura y el estado de carga son demasiado complicados. Para utilizar el método de los elementos finitos, por otro parte, se requiere discretizar el continuo material en un número finito de sectores (elementos finitos con dimensiones finitas), con geometría más simple, interconectados entre sí a través de nudos.

En cierto modo, los elementos finitos son pequeños pedazos de la estructura real. El hecho de idealizar la interconexión entre los elementos finitos exclusivamente a través de sus nudos, podría determinar que solamente en tales nudos se cumplan obligatoriamente las condiciones de compatibilidad de deformación. El resultado que se obtendría es una flexibilización excesiva de la estructura, pues se permitirían traslapes, separaciones o quiebres entre caras de los elementos contiguos.

2

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Es evidente que éste no es el comportamiento de la estructura real, por lo que para un modelamiento más apropiado, los elementos finitos aparentemente deberían deformarse siguiendo elásticas que mantengan la continuidad entre elementos, consiguiéndose de este modo compatibilidad de deformaciones entre las caras adyacentes de los elementos (no siempre ese enfoque es el más conveniente, pero es un buen punto de partida).

Los triángulos y los cuadriláteros planos, constituyen los elementos finitos bidimensionales más utilizados en el análisis estructural, tanto por la facilidad con que se adaptan a casi cualquier configuración geométrica, como por la relativa simplicidad de determinación de sus matrices de rigideces.

Las barras lineales, que conforman las estructuras aporticadas y las celosías, constituyen los elementos finitos naturales. El estudio de las barras lineales ha sido extenso, y los tratados de Análisis Matricial de Estructuras detallan la manera de modelar su comportamiento. 3

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Cuando se presentan continuos tridimensionales (como la presa de la siguiente figura), es usual la utilización de elementos finitos poliédricos, como el hexaedro o el tetraedro.

En el caso de continuos superficiales curvos, se suelen utilizar cuadriláteros fuera de plano (los cuatro vértices del cuadrilátero no pertenecen a un mismo plano).

4

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

2.

ECUACIONES DE LA TEORÍA DE LA ELASTICIDAD:

La Teoría de la Elasticidad es un auxilio importante para comprender el Método de los Elementos Finitos. La siguiente figura representa un elemento diferencial plano de espesor constante “t” (no es un elemento finito pues tiene dimensiones infinitamente pequeñas en lugar de dimensiones finitas).

Las fuerzas por unidad de volumen “Fx” y “Fy”, que actúan sobre el cuerpo, pueden provenir de la acción de la aceleración de la gravedad, aceleraciones sísmicas, campos magnéticos, etc. a)

Ecuaciones Diferenciales de Equilibrio:

Planteando equilibrio de fuerzas, en el elemento diferencial bidimensional, en las direcciones “x” y “y”, se tiene: 5

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

∂τ xy ∂σ x ⋅ dx ⋅ dy ⋅ t + ⋅ dy ⋅ dx ⋅ t + Fx ⋅ dx ⋅ dy ⋅ t = 0 ∂x ∂y ∂τ yx ∂σ y ⋅ dx ⋅ dy ⋅ t + ⋅ dy ⋅ dx ⋅ t + Fy ⋅ dx ⋅ dy ⋅ t = 0 ∂x ∂y Simplificando: ∂σ x ∂τ xy + + Fx = 0 ∂x ∂y ∂τ yx ∂σ y + + Fy = 0 ∂x ∂y Donde: τ yx = τ xy

Lo que transforma las ecuaciones previas en: ∂σ x ∂τ xy + + Fx = 0 ∂x ∂y ∂τ xy ∂σ y + + Fy = 0 ∂x ∂y Por analogía, las ecuaciones diferenciales de equilibrio en un elemento diferencial tridimensional son: ∂σ x ∂τ xy ∂τ xz + + + Fx = 0 ∂x ∂y ∂z ∂τ xy ∂σ y ∂τ yz + + + Fy = 0 ∂x ∂y ∂z ∂τ xz ∂τ yz ∂σ z + + + Fz = 0 ∂x ∂y ∂z b)

Compatibilidad de Deformaciones:

Cuando un cuerpo elástico se deforma, el campo de desplazamientos es continuo, sin que se produzcan aberturas, traslapes o quiebres de la elástica, lo que da lugar a las condiciones de compatibilidad de deformaciones. Al considerar la compatibilidad de deformaciones en el elemento diferencial plano, las tres deformaciones unitarias “ex ”, “ey”, “γ xy”, están interrelacionadas, y son función de dos campos de desplazamientos: u = u ( x , y) v = v( x , y)

De igual manera, al considerar la compatibilidad de deformaciones en el elemento diferencial tridimensional, las seis deformaciones unitarias “ex ”, “ey”, “ez”, “γ xy”, “γ yz”, “γ zx ”, son función de tres campos de desplazamientos: u = u ( x , y , z)

6

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

v = v ( x , y , z) w = w ( x , y , z)

c)

Relación entre Desplazamientos y Deformaciones Unitarias:

La relación existente entre los desplazamientos y las deformaciones unitarias es fundamental en la formulación de la matriz de rigideces de los elementos finitos.

Si se expresa matricialmente la relación entre desplazamientos y deformaciones unitarias para el elemento diferencial bidimensional, se tiene:

 ∂u   ∂     ε x   ∂x   ∂x {ε} =  ε y  =  ∂v  =  0 γ   ∂y    xy   ∂u ∂v  ∂ +   ∂y ∂x   ∂y   

 0  ∂  u  ⋅  ∂y   v ∂  ∂x 

Por analogía, la relación entre desplazamientos y deformaciones unitarias para el caso del elemento diferencial tridimensional es:

7

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

 ∂u   ∂x    ∂v    εx   ε   ∂y   y   ∂w   ε   {ε} =  z  =  ∂u ∂z ∂v  = γ xy    +  γ yz   ∂y ∂x     ∂v ∂w   γ zx   +   ∂z ∂y   ∂u + ∂w   ∂z ∂x  d)

∂  ∂x  0   0  ∂  ∂y  0  ∂   ∂z

0 ∂ ∂y 0 ∂ ∂x ∂ ∂z 0

0   0  ∂ u    ∂z    ⋅ v    0  w     ∂ ∂y  ∂  ∂x 

Relaciones Esfuerzo Unitario - Deformación Unitaria:

Para el caso de materiales ortotrópicos (materiales con características elásticas diferentes en cada una de las tres direcciones ortogonales principales), en continuos tridimensionales, se tienen las siguientes relaciones: εx = + εy = − εz = − γ xy = γ yz =

γ xz =

µ xy 1 µ ⋅ σx − ⋅ σ y − xz ⋅ σ z Ex Ey Ez µ xy

[1]

1 µ ⋅ σ y − xz ⋅ σ z Ey Ez

[2]

µ yz µ xz 1 ⋅ σx − ⋅ σy + ⋅ σz Ex Ey Ez

[3]

Ex

⋅ σx +

τ xy

[4]

G xy τ yz

[5]

G yz

τ xz G xz

[6]

Para el caso de materiales isotrópicos (materiales con características elásticas idénticas en todas las direcciones), en continuos tridimensionales, se tienen las siguientes relaciones simplificadas: 1 µ µ ⋅ σx − ⋅ σy − ⋅ σz E E E µ 1 µ ε y = − ⋅ σx + ⋅ σy − ⋅ σz E E E µ µ 1 ε z = − ⋅ σx − ⋅ σy + ⋅ σz E E E τ xy γ xy = G εx = +

[1’] [2’] [3’] [4’] 8

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

γ yz =

γ xz =

τ yz

[5’]

G τ xz

[6’]

G

Donde: G=

E 2(1 + µ)

Expresando matricialmente las relaciones correspondientes a elementos bidimensionales isotrópicos, bajo condición de esfuerzos planos (se eliminan las ecuaciones “3’ ”, “5’ ” y “6’ ”, y los esfuerzos “s z”, “t yz” y “t zx ”), se tiene:  εx  0   σx   +1 − µ   1    {ε} =  ε y  = − µ + 1 0  ⋅  σ y  = [C] ⋅ {σ} γ  E  0 0 2(1 + µ)  τ xy    xy 

La relación matricial inversa para esfuerzos planos es:

σx  {σ} =  σ y  = E 2 τ  1 − µ  xy 

  0   εx  1 µ   µ 1 0  ⋅  ε y  = [E] ⋅ {ε}    0 0 1 − µ  γ xy   2 

Evidentemente la matriz [E] es la matriz inversa de [C]. [E] = [C]-1 La matriz [C] recibe el nombre de matriz de deformabilidad del material, y la matriz [E] se denomina matriz de elasticidad del material. La relación matricial entre deformaciones unitarias y esfuerzos unitarios, para continuos tridimensionales, con materiales isotrópicos es:

0 0 0   σx   εx   1 −µ −µ ε  − µ 1 − µ 0 0 0   σ y   y  0 0 0   σz   ε z  1 − µ − µ 1  =   ⋅  γ E 0 0 0 2 ( 1 + µ ) 0 0 xy     τ xy   γ yz   0 0 0 0 2(1 + µ ) 0   τ yz         γ xz  0 0 0 0 2(1 + µ)   τ xz   0 La relación matricial inversa es:

9

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

µ µ   1 1− µ 1−µ  µ µ  1  σx  1 − µ 1−µ σ   µ µ  y 1   σz  E (1 − µ) 1 − µ 1 − µ  = τ xy  (1 + µ)(1 − 2µ)  0 0 0  τ yz      0  τ xz  0 0   0 0  0 

0

0

0

0

0

0

1 − 2µ 2(1 − µ)

0

0

1 − 2µ 2(1 − µ)

0

0

   0   εx      εy 0     ⋅  ε z  0  γ xy   γ yz  0  γ xz   1 − 2µ   2(1 − µ)  0

Para el caso de continuos tridimensionales, con materiales isotrópicos, bajo condiciones de deformaciones planas, se descartan las filas “3”, “5” y “6”, de la matriz [E] de 6x6, y se define “ez = 0”, “γ yz = 0” y “γ xz = 0”, obteniéndose:

  σx  µ 0  ε x  1 − µ E    µ 1− µ {σ} =  σ y  = 0   ε y  = [E]{ε}   τ  (1 + µ )(1 − 2µ )  1 − 2µ  γ  0 0  xy   xy   2  Para el caso de placas planas delgadas, cuyo comportamiento está gobernado por el efecto de flexión, las deformaciones de interés son las curvaturas de la superficie neutra de la estructura, y los esfuerzos requeridos son los momentos flectores por unidad de longitud.

  M xx  1 µ 0   æ xx  E ⋅ t3      µ 1 0  æ yy  = [E]{ε}  M yy  = 2   M  12(1 − µ )  1 − 2µ  æ  xy 0 0    xy   2  Donde: Mxx : Momento flector por unidad de longitud alrededor del eje “x”. Myy : Momento flector por unidad de longitud alrededor del eje “y”. Mxy : Momento torsor por unidad de longitud sobre el plano “xy” (alrededor del eje “z”). æ xx =

æ xx = æ xx

∂ 2w ∂x 2 ∂ 2w

∂y 2

:

curvatura de flexión alrededor del eje “x”

:

curvatura de flexión alrededor del eje “y”

∂ 2w = : ∂x.∂y

curvatura de torsión

10

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

3.

INTERPRETACIÓN FÍSICA DEL MÉTODO DE ENSAMBLAJE DIRECTO:

El Método de Ensamblaje Directo, empleado tradicionalmente en el análisis matricial de estructuras aporticadas y en celosía, es la utilización del principio de que, la solicitación (fue rza o momento) nodal que se requiere, para que varios elementos que convergen a un mismo nudo de la estructura tengan una misma magnitud de corrimiento nodal (desplazamiento o rotación compatible, usualmente unitario), es igual a la suma de las solicitaciones que se requieren para conseguir dicha deformación en cada uno de los elementos que convergen al nudo. Dado que los componentes de las matrices de rigideces de cada uno de los elementos de una estructura, son las fuerzas o momentos nodales que se necesitan para mantener una deformación unitaria en uno de los nudos de un elemento estructural, la formación de la matriz de rigideces global de la estructura puede reducirse a la suma selectiva de los componentes de las matrices de rigideces de todos los elementos de una estructura. Este proceso se conoce como Ensamblaje Directo. Los mismos criterios empleados para la utilización del Método de Ensamblaje Directo en el análisis de pórticos y celosías, pueden ser empleados para analizar continuos discretizados mediante elementos finitos.

Elásticas de Deformación Fundamentales de las Barras Planas:

11

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Elásticas de Deformación Correspondientes a los Corrimientos Unitarios de los Grados de Libertad de la Estructura:

Por cada desplazamiento nodal desconocido de la estructura se plantea una ecuación de equilibrio de fuerzas, y por cada rotación nodal desconocida se plantea una ecuación de equilibrio de momentos. Cada componente de la matriz de rigideces de la estructura (matriz de coeficientes del sistema de ecuaciones de equilibrio), se puede obtener directamente de las elásticas de deformación para corrimientos unitarios, sumando las solicitaciones de todas las barras que concurren al nudo donde se está especificando la condición de equilibrio. 12

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Ecuaciones de equilibrio del pórtico: δ x3

δ y3

320 +4 0

θz3

δ x4

δ y4

θz4

= Término Independiente 0 2 ΣFx3

0

600

-320

0

256

0

-1.024

256

-7.50 ΣFy3

600

600 +1.024 256

0

-256

85333

-625 ΣMz3

-320

0

180000 +17066 7 0

0

600

0 ΣFx4

0

-1.024

-256

320 +4 0

-256

-7.50 ΣFy4

0

256

85333

600

600 +1.024 -256

180000 +17066 7

625 ΣMz4

En análisis matricial de estructuras, en lugar de emplear las elásticas de deformación para corrimientos unitarios de los grados de libertad, se calculan matrices de rigideces, en coordenadas globales, para cada barra (cada componente de la matriz de rigideces se calcula en base a corrimientos unitarios en los extremos de barra) y, durante el ensamblaje de la matriz de rigideces de la estructura global se realiza la suma de componentes consistentes de las matrices de rigideces de diferentes elementos. Este proceso es numéricamente equivalente a utilizar las elásticas de deformación, y recibe el nombre de Ensamblaje Directo. Los términos independientes, son las solicitaciones nodales más las solicitaciones de barra transformadas a solicitaciones nodales. ( δx 1 )  4.000  0   - 600 [K1−3 ] =   - 4.000  0   - 600

( δy1 ) ( θz1 ) 0 - 600 600 0 0 180000 0 600 - 600 0 0 90000

(δx 3 ) - 4.000 0 600 4.000 0 600

13

( δy 3 ) (θz 3 ) 0 - 600  - 600 0  0 90000   0 600  600 0   0 180000

(Fx1 ) (Fy1 ) (Mz 1 ) (Fx 3 ) (Fy 3 ) ( Mz 3 )

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

( δx 2 )  4.000  0   - 600 [K 2−4 ] =   - 4.000  0   - 600     [K 2−3 ] =     

( δx 2 ) 320 0 0 − 320 0 0

(δy 2 ) (θz 2 ) ( δx 4 ) 0 - 600 - 4.000 600 0 0 0 180000 600 0 600 4.000 - 600 0 0 0 90000 600 (δy 2 ) (θ z 2 ) 0 0 1.024 256 256 170667 0 0 − 1.024 − 256 256 85333

(δy 4 ) (θ z 4 ) 0 - 600  - 600 0  0 90000   0 600  600 0   0 180000

(Fx 2 ) (Fy 2 ) ( Mz 2 ) (Fx 4 ) (Fy 4 ) ( Mz 4 )

(δx 3 ) ( δy 3 ) (θ z 3 ) − 320 0 0  0 − 1.024 256  0 − 256 85333   320 0 0  0 1.024 − 256   0 − 256 170667

(Fx 2 ) (Fy 2 ) (Mz 2 ) (Fx 3 ) (Fy 3 ) ( Mz 3 )

Si mediante algún proceso especial (luego se describirá tal proceso), se pudieran determinar las rigideces de los elementos finitos que conforman un continuo (por ejemplo un cuadrilátero plano en placas delgadas), no existiría ningún obstáculo para que se construyan elásticas de deformación correspondientes a corrimientos unitarios de nudo, que permitan visualizar físicamente los componentes de las diferentes ecuaciones de equilibrio que deberían plantearse. Como alternativa podrían utilizarse las matrices de rigideces de los elementos finitos, en conjunto con el método de ensamblaje directo, para conseguir el mismo objetivo.

Si se supone que los nudos de la placa solamente admiten un desplazamiento transversal y dos rotaciones como corrimientos, y que la estructura está empotrada en su perímetro, el número total de grados de libertad de la estructura sería de 9. Las elásticas correspondientes a los 3 corrimientos unitarios del nudo “7” de la estructura serían:

14

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

w7 = 1

θx 7 = 1

θy 7 = 1

4.

LA ENERGÍA POTENCIAL Y EL MÉTODO RAYLEIGH - RITZ:

a)

Energía Potencial:

La Energía Potencial de un sistema estructural se designa “pP”, y se puede expresar como función de los corrimientos. Cuando “pP” se minimiza con respecto a los corrimientos, da lugar a ecuaciones de 15

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

equilibrio de la forma

[K ]⋅ {D} = {R} Un sistema estructural es conservativo si, partiendo de una configuración inicial, sufre corrimientos arbitrarios y retorna a la configuración inicial sin efectuar trabajo físico alguno (realiza trabajo nulo). Una configuración o un corrimiento es admisible cuando no viola ni las condiciones internas de compatibilidad, ni las condiciones de borde esenciales.

b)

Energía Potencial en Sistemas con un Grado de Libertad:

Como ejemplo, se puede tomar un resorte suspendido, de longitud “L”, de rigidez axial “k”, en cuyo extremo libre se aplica una fuerza “P”, y se permite un desplazamiento vertical “D” en el lugar de aplicación de la fuerza.

La energía potencial (capacidad de realizar trabajo a futuro) de un sistema estructural tiene dos componentes: Ø Ø

Energía Potencial de Deformación de la estructura. Energía Potencial de las solicitaciones. 16

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Si se toma como nivel de referencia al extremo libre del resorte cuando no está solicitado por la fuerza, la energía potencial del sistema, después de aplicada la fuerza y producidos los corrimientos es: πP =

1 ⋅ k ⋅ D 2 − P.D 2

El signo negativo de la energía potencial de la fuerza obedece a que, una vez realizado el trabajo, la fuerza ha perdido capacidad de realizar trabajo a futuro. Derivando la energía potencial “pP” con respecto a “D”, e igualando a cero para obtener un mínimo, se tiene:

k ⋅D − P = 0 k⋅D = P Esta ecuación es exactamente la misma que se plantearía al imponer condiciones de equilibrio en el sistema. D eq =

P k

Si alternativamente se toma como nivel de referencia a un punto ubicado “H” unidades hacia abajo del extremo libre del resorte cuando no está solicitado por la fuerza, la energía potencial del sistema se describiría como: πP =

1 ⋅ k ⋅ D 2 − P ⋅ ( H − C) 2

Derivando la nueva ecuación de energía potencial con respecto a “D”, e igualando a cero para obtener un mínimo, se tiene:

k ⋅D − P = 0 k⋅D = P Nuevamente se obtiene que: P D eq = k El resultado obtenido es independiente de cualquier nivel de referencia que se escoja para definir la energía potencial del sistema estructural, por lo que resultaría conveniente escoger aquel que defina las expresiones más sencillas o las más convenientes para simplificación. Cualitativamente se puede decir que las solicitaciones pierden energía potencial cuando han realizado trabajo sobre una deformación en la misma dirección que la solicitación, mientras que los resortes almacenan energía potencial positiva sin importar la dirección de la deformación. La representación gráfica de las dos ecuaciones antes detalladas indica que la energía potencial ha sido minimizada, y que los mínimos son coincidentes:

17

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

De todas las configuraciones admisibles de la elástica de deformación de un sistema conservativo, aquellas que satisfacen las ecuaciones de equilibrio convierten a la energía potencial del sistema en estacionaria con respecto a pequeñas variaciones de los corrimientos. c)

Energía Potencial en Sistemas con Múltiples Grados de Libertad:

Se dice que un sistema tiene “n” grados de libertad, si se requieren “n” magnitudes independientes para definir su configuración. En este caso, la energía potencial del sistema será función de la magnitud de los diferentes grados de libertad. π P = F( D1 , D 2 , D3 , ..., D n ) Si se aplica la condición estacionaria de la energía potencial se tiene:

∂π P =0 ∂D i

(i = 1, 2, 3, ..., n )

El resultado es un sistema de “n” ecuaciones simultáneas con “n” incógnitas. Como ejemplo, se pueden tomar tres resortes en serie, de rigideces axiales “k1”, “k2” y “k3”, en cuyos extremos se aplican fuerzas “P1”, “P2” y “P3”, respectivamente, y se permiten desplazamientos absolutos “D1”, “D2” y “D3” en los sitios de aplicación de las fuerzas.

La energía potencial del sistema, después de las deformaciones es: πP =

1 1 1 k 1 ⋅ D 12 + k 2 ⋅ ( D 2 − D1 ) 2 + k 3 ⋅ ( D 3 − D 2 ) 2 − P1 ⋅ D 1 − P2 ⋅ D 2 − P3 ⋅ D 3 2 2 2

Derivando con respecto a cada corrimiento, e igualando a “0”, para minimizar la energía potencial del sistema, se tiene: 18

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

∂π P = k1 ⋅ D1 − k 2 ⋅ ( D 2 − D1 ) − P1 = 0 ∂D1 ∂π P = k 2 ⋅ ( D 2 − D1 ) − k 3 ⋅ ( D3 − D 2 ) − P2 = 0 ∂D 2 ∂π P = k 3 ⋅ ( D 3 − D 2 ) − P3 = 0 ∂D 3 Organizando matricialmente el sistema de ecuaciones se tiene:

k 1 + k 2  −k 2   0

− k2 k2 + k3 − k3

0   D1   P1      − k 3  ⋅ D 2  = P2  k 3   D3   P3 

La expresión matricial simplificada es: [K ] ⋅ {D} = {P} El sistema de ecuaciones es exactamente igual al que se obtendría planteando ecuaciones de equilibrio en los puntos de aplicación de las fuerzas. Una manera alternativa de plantear la ecuación de energía potencial, en términos matriciales, es: πP =

1 {D}T ⋅ [K] ⋅ {D} − {D}T .{P} 2

En el caso del ejemplo previo, la expresión desarrollada de la ecuación de energía potencial es:

1 π P = [D1 2 d)

D2

k 1 + k 2 D3 ] ⋅  − k 2  0

− k2 k 2 + k3 − k3

0   D1    − k 3  ⋅ D 2  − [D1 k 3  D 3 

D2

 P1    D3 ].P2  P   3

Expresiones para la Energía Potencial:

Se puede tomar como referencia el caso general de esfuerzos tridimensionales { s } y deformaciones tridimensionales {e}: {σ} = σ x σ y σ z τ xy τ yz τ zx T

[ {ε} = [ε x

]

εy

εz

γ xy

γ yz

γ zx

]T

La relación esfuerzo unitario - deformación unitaria en coordenadas rectangulares es:

{σ} = [E ] ⋅ {ε} Si se define como “Uo” a la energía potencial de deformación por unidad de volumen, tal magnitud representa el trabajo realizado por las fuerzas internas. En un cubo de dimensiones unitarias, el esfuerzo unitario es igual a la fuerza y la deformación unitaria es igual al desplazamiento sobre el que actúa el esfuerzo unitario. Con estas consideraciones, incluyendo todos los esfuerzos unitarios, las deformaciones unitarias infinitesimales producen un cambio en la energía de deformación interna de acuerdo a la siguiente expresión. dU o = {σ}⋅ {dε}T 19

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

dU o = σ x ⋅ dε x + σ y ⋅ dε y + σ z ⋅ dε z + τ xy ⋅ dγ xy + τ yz ⋅ dγ yz + τ zx ⋅ dγ zx

Derivando parcialmente con relación a cada deformación unitaria se tienen las siguientes expresiones.

dU o = σx ∂ε x dU o = σy ∂ε y dU o = σz ∂ε z dU o = τ xy ∂γ xy dU o = τ yz ∂γ yz

dU o = τ zx ∂γ zx Generalizando las expresiones anteriores, con ecuaciones matriciales:  dUo    = {σ} = [E] ⋅ {ε}  ∂ε  Integrando con respecto a las deformaciones unitarias se tiene: 1 U o = {ε}T ⋅ [E] ⋅ {ε} 2 Se definen a los desplazamientos de un punto arbitrario de coordenadas “x”, “y”, “z” con la siguiente expresión:

{f } = [u

v w ]T

Donde “u”, “v”, “w” son función de las coordenadas “x”, “y”, “z”. Las fuerzas por unidad de volumen pierden potencial cuando ocurren los desplazamientos en la misma dirección de las fuerzas. En un volumen unitario el cambio de energía potencial es: Cambio de Potencial = − Fx ⋅ u − Fy ⋅ v − Fz ⋅ w

Un cuerpo de volumen “V” tiene una energía potencial total: πP =

T T ∫V U o .dV − ∫V {f } ⋅ {F}⋅ dV − {D} ⋅ {P}

La primera expresión es la energía de deformación interna, la segunda es el cambio de potencial en las fuerzas volumétricas, y la tercera es el cambio de potencial en las fuerzas que actúan sobre los nudos. Previamente se estableció que, para el caso de sistemas estructurales sin fuerzas volumétricas, la energía potencial podía calcularse con la siguiente expresión: 20

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

πP =

1 {D}T ⋅ [K] ⋅ {D} − {D}T .{P} 2

Comparando las dos ecuaciones se deduce que: 1 {D}T ⋅ [K] ⋅ {D} = ∫V U o ⋅ dV 2

Reemplazando “Uo” en la expresión anterior: 1 {D}T ⋅ [K] ⋅ {D} = 1 ∫V {ε}T ⋅ [E ] ⋅ {ε}⋅ dV 2 2

Simplificando:

{D}T ⋅ [K ]⋅ {D} = ∫V {ε}T ⋅ [E] ⋅ {ε}⋅ dV Esta relación permite la determinación de la matriz de rigideces de un elemento finito, en función de sus deformaciones unitarias internas. e)

El Método Rayleigh - Ritz:

Las estructuras con miembros discretos, como los pórticos y celosías, tienen un número finito de grados de libertad, pero los sistemas continuos pueden tener grados de libertad en cada uno de sus puntos, y su comportamiento se describe mediante ecuaciones diferenciales parciales simultáneas. Se puede evitar resolver dichas ecuaciones (en la gran mayoría de los casos no tienen solución cerrada), empleando el Método Rayleigh - Ritz, que utiliza expresiones matemáticas de interpolación para expresar los corrimientos de cada punto, en función de un número finito de grados de libertad. El Método Rayleigh - Ritz se vuelve más exacto mientras mayor sea el número de grados de libertad que se utilice.

5.

EL ELEMENTO FINITO CUADRILÁTERO PLANO PARA MODELAR DEFORMACIONES FLEXIONANTES EN PLACAS DELGADAS:

A continuación se presenta una metodología genérica para formular la matriz de rigideces de un elemento finito cuadrilátero plano, empleado en el análisis estructural de placas delgadas. Procedimientos muy similares se emplean en la deducción de matrices de rigideces de otros tipos de elementos finitos, utilizados en otros tipos de problemas estructurales.

21

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

a)

El Cuadrilátero Plano en Coordenadas Globales:

El elemento finito cuadrilátero plano puede tener una geometría real arbitraria.

b)

El Cuadrilatero Plano en Coordenadas Naturales:

Para efectos de simplificar las operaciones se utiliza como referencia al elemento finito cuadrilátero plano en coordenadas normalizadas (coordenadas naturales).

22

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

A cada punto del cuadrilátero plano real le corresponde un punto del cuadrilátero normalizado con coordenadas naturales. Las ecuaciones de transformación entre los sistemas de coordenadas se discuten posteriormente. c)

Grados de Libertad del Cuadrilátero Plano:

Los grados de libertad (corrimientos) del cuadrilátero plano utilizado en el modelamiento de placas, tanto en coordenadas globales como en coordenadas naturales, son un desplazamiento “w” transversal al plano principal y dos rotaciones (“θx=∂w/∂y”, “θy=-∂w/∂x”) por cada nudo, lo que significa un total de 12 corrimientos referenciales para el elemento finito. d)

Funciones de Forma de los Desplazamientos Nodales en el Cuadrilátero Plano en Coordenadas Naturales:

Se definen las siguientes Funciones de Forma de los Corrimientos Nodales, cuya característica es la de ser funciones de dos variables (“r”, “s”) simples y manejables, que tienen valor unitario para uno de los grados de libertad de los nudos del elemento finito y valor nulo para los restantes once grados de libertad. Por facilidad de formulación se utiliza como base al elemento finito en coordenadas naturales, y posteriormente se realiza una transformación consistente de coordenadas para modelar los corromientos en el elemento finito en coordenadas reales. Los desplazamientos transversales en el elemento finito se describen mediante la variable “w”, la rotación de nudo alrededor del eje “x” es “θx=∂w/∂y”, y la rotación de nudo alrededor del eje “y” es “θy=-∂w/∂x”. Desplazamiento Unitario Perpendicular al Nudo “1” (w1=1): w=

1 3 ( r − 3r + 2)(s 3 − 3s + 2) 16

23

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Rotación Unitaria Alrededor del Eje “x” en el Nudo “1” (θx1=∂w1/∂s=1): w=

1 3 ( r − 3r + 2)(s 3 − s 2 − s + 1) 16

Rotación Unitaria Alrededor del Eje “y” en el Nudo “1” (θy1=-∂w1/∂r=1): w=

1 ( −r 3 + r 2 + r − 1)(s 3 − 3s + 2) 16

24

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Desplazamiento Unitario Perpendicular al Nudo “2” (w2=1): w=

1 ( −r 3 + 3r + 2)(s 3 − 3s + 2) 16

Rotación Unitaria Alrede dor del Eje “x” en el Nudo “2” (θx2=∂w2/∂s=1): w=

1 ( −r 3 + 3r + 2)(s 3 − s 2 − s + 1) 16

25

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Rotación Unitaria Alrede dor del Eje “y” en el Nudo “2” (θy2=-∂w2/∂r=1): w=

1 ( −r 3 − r 2 + r + 1)(s 3 − 3s + 2) 16

Desplazamiento Unitario Perpendicular en el Nudo “3” (w3=1): w=

1 ( −r 3 + 3r + 2)( −s 3 + 3s + 2) 16

26

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Rotación Unitaria Alrededor del Eje “x” en el Nudo “3” (θx3=∂w3/∂s=1): w=

1 ( −r 3 + 3r + 2)(s 3 + s 2 − s − 1) 16

Rotación Unitaria Alrededor del Eje “y” en el Nudo “3” (θy3=-∂w3/∂r=1): w=

1 ( −r 3 − r 2 + r + 1)( −s 3 + 3s + 2) 16

27

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Desplazamiento Unitario Perpendicular al Nudo “4” (w4=1): w=

1 3 ( r − 3r + 2)( −s 3 + 3s + 2) 16

Rotación Unitaria Alrede dor del Eje “x” en el Nudo “4” (θx4=∂w4/∂s=1): w=

1 3 ( r − 3r + 2)(s 3 + s 2 − s − 1) 16

28

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Rotación Unitaria Alrede dor del Eje “y” en el Nudo “4” (θy4=-∂w4/∂r=1): w=

1 ( −r 3 + r 2 + r − 1)( −s 3 + 3s + 2) 16

f)

Funciones de Transformación de Coordenadas:

Para transformar coordenadas entre el elemento finito en coordenadas naturales y el elemento finito en coordenadas globales se utilizan las siguientes funciones paramétricas cuyo valor es unitario para las coordenadas de uno de los nudos y nulo para los tres restantes: 1 3 ( r − 3r + 2)(s 3 − 3s + 2) 16 1 N II = (− r 3 + 3r + 2)(s 3 − 3s + 2) 16 1 N III = ( −r 3 + 3r + 2)( −s 3 + 3s + 2) 16 1 N IV = (r 3 − 3r + 2)( −s 3 + 3s + 2) 16 NI =

Las cuatro funciones paramétricas detalladas son coincidentes con las funciones que describen los cuatro desplazamientos transversales de nudo. Las expresiones de transformación entre los dos sistemas de coordenadas, que utilizan las funciones paramétricas detalladas anteriormente son: x = N I ⋅ x 1 + N II ⋅ x 2 + N III ⋅ x 3 + N IV ⋅ x 4 y = N I ⋅ y1 + N II ⋅ y 2 + N III ⋅ y 3 + N IV ⋅ y 4 Reemplazando se tiene:

29

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador 3 3 3 3  1 ( r − 3r + 2)(s − 3s + 2) x 1 + ( −r + 3r + 2)( s − 3s + 2) x 2 x=   16 + ( −r 3 + 3r + 2)( −s 3 + 3s + 2) x 3 + (r 3 − 3r + 2)( −s 3 + 3s + 2) x 4  3 3 3 3  1 ( r − 3r + 2)(s − 3s + 2) y1 + ( −r + 3r + 2)( s − 3s + 2) y 2 y=   16 + ( −r 3 + 3r + 2)( −s 3 + 3s + 2) y 3 + (r 3 − 3r + 2)( −s 3 + 3s + 2) y 4 

g)

Campo de Desplazamientos:

En el literal “d” se detallan doce funciones de forma para describir los campos de desplazamientos transversales del elemento finito, basadas en los desplazamientos transversales y rotaciones de nudo. 1 3 (r − 3r + 2)(s 3 − 3s + 2) 16 1 = ( r 3 − 3r + 2)(s 3 − s 2 − s + 1) 16 1 = ( −r 3 + r 2 + r − 1)( s 3 − 3s + 2) 16 1 = ( −r 3 + 3r + 2)( s 3 − 3s + 2) 16 1 = ( −r 3 + 3r + 2)( s 3 − s 2 − s + 1) 16 1 = (− r 3 − r 2 + r + 1)(s 3 − 3s + 2) 16 1 = ( −r 3 + 3r + 2)( −s 3 + 3s + 2) 16 1 = ( −r 3 + 3r + 2)( s 3 + s 2 − s − 1) 16 1 = (− r 3 − r 2 + r + 1)( −s 3 + 3s + 2) 16

N1 = N2 N3 N4 N5 N6 N7 N8 N9

N 10 =

1 3 (r − 3r + 2)( −s 3 + 3s + 2) 16

N 11 =

1 3 (r − 3r + 2)( s 3 + s 2 − s − 1) 16

N 12 =

1 (− r 3 + r 2 + r − 1)( −s 3 + 3s + 2) 16

El campo de desplazamientos transversales en el elemento finito normalizado es una función de los corrimientos nodales. La contribución de cada corrimiento nodal está definida por la respectiva función de forma. w = N1 ⋅ w1 + N 2 ⋅ θx 1 + N3 ⋅ θy1 + N 4 ⋅ w 2 + N5 ⋅ θx 2 + N 6 ⋅ θy 2 + N 7 ⋅ w 3 + N 8 ⋅ θ x 3 + N9 ⋅ θy 3 + N10 ⋅ w 4 + N11 ⋅ θx 4 + N12 ⋅ θy 4 Se podrían utilizar formulaciones más sencillas (menos elaboradas) que incluyan como únicos grados de libertad los desplazamientos transversales de nudo, en dicho caso existirán solamente cuatro 30

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

funciones de forma que describan los desplazamientos, y las mismas ecuaciones deben emplearse para las transformaciones de coordenadas. 1 (1 − r )(1 − s) 4 1 N II = (1 + r )(1 − s) 4 1 N III = (1 + r )(1 + s) 4 1 N IV = (1 − r )(1 + s) 4 NI =

Basado en las últimas funciones de forma, la relación entre coordenadas sería: 1 {(1 − r )(1 − s ) x1 + (1 + r )(1 − s) x 2 + (1 + r )(1 + s ) x 3 + (1 − r )(1 + s) x 4 } 4 1 y = {(1 − r )(1 − s) y1 + (1 + r )(1 − s ) y 2 + (1 + r )(1 + s ) y 3 + (1 − r )(1 + s) y 4 } 4 x=

h)

Relaciones Deformaciones Unitarias - Desplazamientos:

Por el principio de Kirchoff, que es una extensión de una hipótesis básica de flexión en vigas, cualquier línea recta perpendicular al plano principal de una placa, antes de la aplicación de solicitaciones, mantiene su condición de recta después de la deformación provocada por las solicitaciones, lo que determina las siguientes equivalencias: u = z ⋅ θy v = −z ⋅ θx

De las seis componentes de deformación unitaria en problemas tridimensionales (ex , ey, ez, γ xy, γ yz, γ zx ). En el presente caso ez se considera despreciable y se asume nula. εz =

∂w =0 ∂z

Las restantes deformaciones unitarias se describen con las siguientes igualdades: ∂u ∂x ∂v εy = ∂y ∂u ∂v γ xy = + ∂y ∂x ∂v ∂w γ yz = + ∂z ∂y ∂u ∂w γ xz = + ∂z ∂x εx =

Se requiere transformar las derivadas parciales respecto a las coordenadas x, y, z del elemento finito en coordenadas reales, tomando como base las derivadas respecto a r, s, t del elemento finito en coordenadas normalizadas, para lo que se utilizan las expresiones de la derivación en cadena para 31

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

funciones de funciones. ∂u ∂u ∂x ∂u ∂y ∂u ∂z = ⋅ + ⋅ + ⋅ ∂r ∂x ∂r ∂y ∂r ∂z ∂r ∂u ∂u ∂x ∂u ∂y ∂u ∂z = ⋅ + ⋅ + ⋅ ∂s ∂x ∂s ∂y ∂s ∂z ∂s ∂u ∂u ∂x ∂u ∂y ∂u ∂z = ⋅ + ⋅ + ⋅ ∂t ∂x ∂t ∂y ∂t ∂z ∂t Para el caso de placas la coordenada perpendicular al plano es independiente de las coordenadas coplanares, y es conveniente que sean coincidentes tanto en coordenadas reales como en coordenadas normalizadas. z=t

De donde: ∂z ∂r ∂z ∂s ∂z ∂t ∂x ∂t ∂y ∂t

=0 =0 =1 =0 =0

Reemplazando en las derivaciones en cadena, y simplificando, se tiene: ∂u ∂u ∂x ∂u ∂y = ⋅ + ⋅ ∂r ∂x ∂r ∂y ∂r ∂u ∂u ∂x ∂u ∂y = ⋅ + ⋅ ∂s ∂x ∂s ∂y ∂s ∂u =1 ∂t Se puede emplear un procedimiento similar para relacionar las derivadas parciales de v, w, respecto a r, s, t. Si se agrupan matricialmente las expresiones, se obtiene:

 ∂u   ∂s     ∂u   =  ∂y   ∂u   ∂x   

 ∂x  ∂r  ∂x   ∂s 0 

∂y ∂r ∂y ∂s 0

  ∂u  0   ∂x   ∂u  0  ⋅    ∂y  1   ∂u    ∂z 

Expresiones similares pueden obtenerse para las derivadas parciales de v, w respecto a r, s, t. 32

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

La matriz de transformación obtenida recibe el nombre de matriz Jacobiana.  ∂x  ∂r  [J ] =  ∂x  ∂s 0 

∂y ∂r ∂y ∂s 0

0   0  1  

La inversa de la matriz Jacobiana puede representarse:

[Γ] = [J ]

−1

 Γ11 Γ12 = Γ21 22  0 0

0  0  1 

Las relaciones inversas agrupadas se pueden expresar:

  u ′x    x ′r  u′  x′  y   s  u ′z    0     v′x      v′y  =   v′    z   w′x   w ′    y   w ′z   

y ′r x ′s 0

0 0 1

−1

 x ′r x′  s  0

y′r x ′s 0

0 0 1

−1

 x ′r x′  s  0

y′r x ′s 0

 u′   r    u s′     u ′t    v′r      ⋅  v′s    v′   t −1  0  w ′r     0  w s′  1  w ′t  

Las derivadas de u, v, w respecto a r, s, t (u’y es la derivada parcial de “u” respecto a “y” o “∂u/∂y”)se calculan con las funciones de forma de w y las expresiones que relacionan las rotaciones de nudo con los corrimientos u, v.

 u ′r   0  u′   0  s   u ′t   0     v′r  n  0    v′s  = ∑  0  v ′  i =1  0   t  − z ⋅ Ni′  w′r  r   w′  ′ − z ⋅ N i  s  s   w ′t   − Ni

− z ⋅ Ni ′r − z ⋅ Ni ′s − Ni 0 0 0 0 0 0

0  0  0   − z ⋅ Ni ′r   wi    − z ⋅ Ni ′s  ⋅ θ xi   − Ni  θyi  0   0  0 

Donde n es el número de funciones de forma (12 para la solución más exacta y 4 para la solución simplificada). Se debe recordar que: 33

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

ε x = u ′x ε y = v′y ε z = w ′z =0 γ xy = u ′y + v′x γ yz = v′z + w ′y

γ xz = u ′z + w ′x Condensando las últimas expresiones para obtener matrices de cinco filas se tiene:   εx     εy   n    γ xy  = Σ  i =1 γ    yz    γ xz  

0

0

0

0

0

0

bi

0

ai

− Ni

0     0   0 +   − Ni     0  

0

− z.a i

0

0

0

− z.b i

0

0

0

0

  − z.b i   w i    − z.a i  θ xi   0  θ yi   0  0

O en su forma simplificada: {ε} = [B] ⋅ {d} Donde: a i = Γ11 ⋅ Ni ′r + Γ12 ⋅ Ni ′s b i = Γ21 ⋅ Ni ′r + Γ22 ⋅ Ni ′s i)

Derivadas Parciales de las Funciones de Forma:

A continuación se presenta un resumen de las derivadas parciales que se requieren para determinar la matriz de rigideces de los elementos finitos cuadriláteros en placas. Las derivadas parciales que aparecen en la Matriz Jacobiana, y que son utilizadas con las 4 funciones de forma nodales, son: 2 3 2 3 2 3 ∂x 1 (3r − 3)( s − 3s + 2) x 1 + ( −3r + 3)(s − 3s + 2) x 2 + (−3r + 3)( −s + 3s + 2) x 3 + =   ∂r 16 + (3r 2 − 3)( −s 3 + 3s + 2) x 4  2 3 2 3 2 3 ∂y 1 (3r − 3)(s − 3s + 2) y1 + ( −3r + 3)( s − 3s + 2) y 2 + ( −3r + 3)( −s + 3s + 2) y 3 +  =   ∂r 16 + (3r 2 − 3)( −s 3 + 3s + 2) y 4  3 2 3 2 3 2 ∂x 1 ( r − 3r + 2)(3s − 3) x 1 + ( −r + 3r + 2)( 3s − 3) x 2 + (− r + 3r + 2)( −3s + 3) x 3 + =   ∂s 16 + ( r 3 − 3r + 2)( −3s 2 + 3) x 4  3 2 3 2 3 2 ∂y 1 (r − 3r + 2)( 3s − 3) y1 + ( −r + 3r + 2)(3s − 3) y 2 + ( −r + 3r + 2)( −3s + 3) y 3 +  =   ∂s 16 + ( r 3 − 3r + 2)( −3s 2 + 3) y 4  Las derivadas parciales de las funciones de forma respecto a r son: 1 N 1′ r = ( 3r 2 − 3)(s 3 − 3s + 2) 16 1 N 2′ r = (3r 2 − 3)(s 3 − s 2 − s + 1) 16

34

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

N 3′ r = N 4′ r N 5′ r N 6′ r N 7′ r N 8 ′r N 9′ r

1 ( −3r 2 + 2 r + 1)(s 3 − 3s + 2) 16 1 = (−3r 2 + 3)(s 3 − 3s + 2) 16 1 = ( −3r 2 + 3)(s 3 − s 2 − s + 1) 16 1 = ( −3r 2 − 2r + 1)( s 3 − 3s + 2) 16 1 = (−3r 2 + 3)( −s 3 + 3s + 2) 16 1 = ( −3r 2 + 3)( s 3 + s 2 − s − 1) 16 1 = ( −3r 2 − 2r + 1)( −s 3 + 3s + 2) 16

N 10′ r =

1 ( 3r 2 − 3)( −s 3 + 3s + 2) 16

1 N 11′ r = ( 3r 2 − 3)(s 3 + s 2 − s − 1) 16 N 12′ r =

1 ( −3r 2 + 2r + 1)( −s 3 + 3s + 2) 16

Las derivadas parciales de las funciones de forma respecto a s son: 1 N 1′ s = ( r 3 − 3r + 2)( 3s 2 − 3) 16 1 N 2′ s = (r 3 − 3r + 2)( 3s 2 − 2s − 1) 16 1 N 3′s = ( −r 3 + r 2 + r − 1)(3s 2 − 3) 16 1 N 4′ s = (− r 3 + 3r + 2)( 3s 2 − 3) 16 1 N 5′s = ( −r 3 + 3r + 2)( 3s 2 − 2s − 1) 16 1 N 6′ s = (− r 3 − r 2 + r + 1)(3s 2 − 3) 16 1 N 7′ s = (− r 3 + 3r + 2)( −3s 2 + 3) 16 1 N 8 ′s = ( −r 3 + 3r + 2)(3s 2 + 2s − 1) 16 1 N 9′ s = ( −r 3 − r 2 + r + 1)( −3s 2 + 3) 16 N 10′ s =

1 3 ( r − 3r + 2)( −3s 2 + 3) 16

35

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

1 N 11′ s = ( r 3 − 3r + 2)( 3s 2 + 2s − 1) 16 N 12′ s =

1 ( −r 3 + r 2 + r − 1)( −3s 2 + 3) 16

A partir de estas expresiones y de la inversa de la Matriz Jacobiana se obtienen las derivadas de las 4 funciones de forma nodales (N1, N4, N7, N10) respecto a las variables (x, y), quedando establecida la matriz [B]. Para las 8 funciones complementarias de flexión (N2, N3, N5, N6, N8, N9, N11, N12) se calcula la inversa de la Matriz Jacobiana, evaluando solamente las derivadas de las funciones antes detalladas para valores de (r=0; s=0). j)

Relaciones Esfuerzos Unitarios - Deformaciones Unitarias:

La matriz de elasticidad que incluye las deformaciones por corte es:

0 0  E ′x E ′′ 0  E ′′ E ′ 0 0 0  y  [E ] =  0 0 G 0 0    0 0 G yz 0  0  0 0 0 0 G xz  Donde:

E ′x = E′y =

E

1−µ2 µ⋅ E

E ′′ =

2(1 + µ 2 ) E G= 2(1 + µ) De la teoría de Elasticidad se conoce la siguiente relación para el caso de placas planas: {s } = [E].{e} {s } = [E].[B].{d} k)

Matriz de Rigideces del Elemento Finito:

La siguiente expresión relaciona la matriz de rigideces de un elemento finito genérico con sus deformaciones unitarias internas:

{d}T ⋅ [K] ⋅ {d} = ∫V {ε}T ⋅ [E] ⋅ {ε}.dV Donde:

{ε} = [B]⋅ {d} {ε}T = {d}T ⋅ [B]T Reemplazando {e} y {e}T se tiene: 36

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

{d}T ⋅ [K] ⋅ {d} = ∫V {d}T ⋅ [B]T ⋅ [E] ⋅ [B] ⋅ {d}⋅ dV Simplificando:

[K ] = ∫V [B]T ⋅ [E] ⋅ [B] ⋅ dV Efectuando la integración en coordenadas normalizadas (se utiliza el determinante de la matriz jacobiana), la matriz de rigideces del elemento finito cuadrilátero plano queda definida como:

[K ] = ∫−1 ∫−1 [B]T ⋅ [E] ⋅ [B] ⋅ Det [J ] ⋅ dV 1

1

Donde: dVn : diferencial volumétrico en el elemento finito normalizado dVn = espesor.dAn

[K ] = ∫−1 ∫−1 [B]T ⋅ [E] ⋅ [B] ⋅ Det [J ] ⋅ espesor ⋅ dAn 1

1

Para integrar numéricamente la expresión se pueden utilizar 4 puntos de integración (un punto de Gauss-Legendre por cada cuadrante), cuyas coordenadas son: (0.57735,0.57735), (0.57735,0.57735), (-0.57735,-0.57735) y (0.57735,-0.57735), o ( 3 / 3 , 3 / 3 )...

[K ] = ∑ [[Bi ]T ⋅ [Ei] ⋅ [Bi ] ⋅ Det [J ] ⋅ espesor ⋅ Ai ] 4

i =1

El área de influencia de cada punto de integración (Ai) es unitaria, por lo tanto:

[K ] = ∑ [[Bi ]T ⋅ [Ei] ⋅ [Bi] ⋅ Det[J ] ⋅ espesor ] 4

i =1

La matriz de rigideces obtenida es de 12x12 Para obtener una mayor precisión en la integración numérica se podría emplear un mayor número de puntos de integración, con coordenadas y áreas de influencia descritas por los coeficientes de Gauss para los polinomios de Legendre (Gauss-Legendre). Las coordenadas de los puntos de Gauss y sus respectivos pesos de integración (áreas de influencia), para diferentes números de puntos de integración son:

37

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

Número de Puntos de Integración 4

9

16

l)

Coordenadas

Área de Influencia de cada punto (Ai)

(-0.57735027,-0.57735027) (+0.57735027,-0.57735027) (+0.57735027,+0.57735027) (-0.57735027,+0.57735027) (-0.77459667,-0.77459667) (-0.77459667,0.00000000) (-0.77459667,+0.77459667) (0.00000000,-0.77459667) (0.00000000,0.00000000) (0.00000000,+0.77459667) (+0.77459667,-0.77459667) (+0.77459667,0.00000000) (+0.77459667,+0.77459667) (-0.86113631,-0.866113631) (-0.86113631,-0.33998104) (-0.86113631,+0.33998104) (-0.86113631,+0.866113631) (-0.33998104,-0.866113631) (-0.33998104,-0.33998104) (-0.33998104,+0.33998104) (-0.33998104,+0.866113631) (+0.33998104,-0.866113631) (+0.33998104,-0.33998104) (+0.33998104,+0.33998104) (+0.33998104,+0.866113631) (+0.86113631,-0.866113631) (+0.86113631,-0.33998104) (+0.86113631,+0.33998104) (+0.86113631,+0.866113631)

1.0000000000x1.0000000000=1.0000000000 1.0000000000x1.0000000000=1.0000000000 1.0000000000x1.0000000000=1.0000000000 1.0000000000x1.0000000000=1.0000000000 0.5555555556x0.5555555556=0.3086419753 0.5555555556x0.8888888889=0.4938271605 0.5555555556x0.5555555556=0.3086419753 0.8888888889x0.5555555556=0.4938271605 0.8888888889x0.8888888889=0.7901234568 0.8888888889x0.5555555556=0.4938271605 0.5555555556x0.5555555556=0.3086419753 0.5555555556x0.8888888889=0.4938271605 0.5555555556x0.5555555556=0.3086419753 0.3478548451x0.3478548451=0.1210029933 0.3478548451x0.6521451549=0.2268518518 0.3478548451x0.6521451549=0.2268518518 0.3478548451x0.3478548451=0.1210029933 0.6521451549x0.3478548451=0.2268518518 0.6521451549x0.6521451549=0.4252933031 0.6521451549x0.6521451549=0.4252933031 0.6521451549x0.3478548451=0.2268518518 0.6521451549x0.3478548451=0.2268518518 0.6521451549x0.6521451549=0.4252933031 0.6521451549x0.6521451549=0.4252933031 0.6521451549x0.3478548451=0.2268518518 0.3478548451x0.3478548451=0.1210029933 0.3478548451x0.6521451549=0.2268518518 0.3478548451x0.6521451549=0.2268518518 0.3478548451x0.3478548451=0.1210029933

Los Términos de Carga de las Ecuaciones de Equilibrio:

Los términos independientes de las ecuaciones de equilibrio, al igual que en el Análisis Matricial de Estructuras Aporticadas y en Celosía, son las solicitaciones nodales que actúan sobre la estructura, más las solicitaciones sobre las caras transformadas a solicitaciones nodales equivalentes. Para determinar las solicitaciones nodales equivalentes se puede igualar el trabajo virtual de las solicitaciones sobre las caras, al trabajo virtual de las solicitaciones nodales equivalentes, tomando como deformación virtual en ambos casos a la elástica de deformación genérica de la cara del elemento. 6.

PROGRAMA DE COMPUTACIÓN Y EJEMPLOS:

6.1

Programa de Análisis de Placas Planas con Elementos Finitos Cuadriláteros Planos:

A continuación se presenta un programa ilustrativo del uso de elementos finitos en el análisis de 38

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador

placas, escrito en lenguaje GWBASIC, cuyos algoritmos pueden ser fácilmente adaptados a cualquier lenguaje científico. 10 REM PROGRAMA DE ANALISIS DE PLACAS PLANAS CON ELEMENTOS FINITOS CUADRANGULARES 20 REM DESARROLLADO POR MARCELO ROMO EN NOVIEMBRE DE 1993 30 A1$=" ### ######.### ######.### ## ## ##" 40 A2$=" ### ### ### #######.### ######.### ########.###" 50 A3$=" ### #.######^^^^ #.######^^^^ #.######^^^^" 60 A4$=" ### #.######^^^^ #.######^^^^ #.######^^^^" 70 A5$=" ### #.###^^^^ #.###^^^^ #.###^^^^ #.###^^^^ #.###^^^^ #.###^^^^" 80 A6$=" ### ######.### ######.### ######.###" 90 A7$=" ### #####.#### #####.####" 100 A8$=" ### ### ### ### ### #######.## ##.### ####.###" 110 REM LEE NOMBRE DEL ARCHIVO DE ENTRADA DE DATOS 120 CLS:LOCATE 10,1:PRINT "DEME NOMBRE DEL ARCHIVO DE ENTRADA DE DATOS"; 130 INPUT ARCHIVO$ 140 REM ABRE ARCHIVOS DE ENTRADA DE DATOS Y SALIDA DE RESULTADOS 150 CLS:LOCATE 10,1:PRINT "ABRE ARCHIVOS DE ENTRADA DE DATOS Y SALIDA DE RESULTADOS" 160 OPEN ARCHIVO$ FOR INPUT AS#1 LEN=128 170 A$=ARCHIVO$+".RES" 180 OPEN A$ FOR OUTPUT AS#2 LEN=128 190 REM LEE E IMPRIME TITULO DEL PROBLEMA 200 INPUT#1,A$ 210 INPUT#1,TITULO$ 220 PRINT#2,TITULO$ 230 PRINT#2," " 240 PRINT#2," DATOS DE LA ESTRUCTURA:" 250 CLS:LOCATE 10,1:PRINT "LEE DATOS DE LA ESTRUCTURA" 260 INPUT#1,A$,A$,A$ 270 REM LEE CARACTERISTICAS BASICAS DE LA ESTRUCTURA Y DIMENSIONA LOS ARREGLOS: 280 REM LEE NUMERO DE NUDOS, NUMERO DE ELEMENTOS FINITOS, NUMERO DE ESTADOS DE CARGA 290 INPUT#1,NNUDOS,NFINITOS,NCARGAS 300 DIM X(NNUDOS),Y(NNUDOS),ORDEN(NNUDOS,3),P1#(NNUDOS,3),CORRIM#(NNUDOS,3) ,P#(3*NNUDOS),PUN(3*NNUDOS) 310 DIM NUDO1F(NFINITOS),NUDO2F(NFINITOS),NUDO3F(NFINITOS),NUDO4F(NFINITOS) ,EF#(NFINITOS),POISSON#(NFINITOS),ESPESOR#(NFINITOS) 320 DIM KM#(12,12),IT(12),F#(12),REACCION#(12),CORR#(12),F1#(12),RE#(12) 330 DIM ESFUERZOS#(3),JCB#(3,3),JI#(3,3),B#(5,12),ELAS#(5,5),PROD#(12,5),R( 16),S(16),P(16),EPSILON#(3),NUDO(4) 340 PRINT#2,"NUMERO DE NUDOS =";NNUDOS 350 PRINT#2,"NUMERO DE ELEMENTOS FINITOS =";NFINITOS 360 PRINT#2,"NUMERO DE ESTADOS DE CARGA =";NCARGAS 370 PRINT#2," " 380 PRINT#2," CARACTERISTICAS DE NUDO:" 390 PRINT#2," NUDO COORDENADAS RESTRICCIONES" 400 PRINT#2," X Y DES.Z ROT.X

39

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador ROT.Y" 410 REM LEE COORDENADAS Y RESTRICCIONES DE NUDO 420 REM # NUDO, COORDENADA X, COORDENADA Y, RESTRICCION DESPL. Z, RESTRICCION ROTAC. X, RESTRICCION ROTAC. Y 430 CLS:LOCATE 10,1:PRINT "LEE COORDENADAS Y RESTRICCIONES DE NUDO" 440 INPUT#1,A$,A$,A$,A$,A$,A$,A$ 450 FOR I=1 TO NNUDOS 460 INPUT#1,J,X(J),Y(J),ORDEN(J,1),ORDEN(J,2),ORDEN(J,3) 470 PRINT#2,USING A1$;J,X(J),Y(J),ORDEN(J,1),ORDEN(J,2),ORDEN(J,3) 480 NEXT I 490 PRINT #2," " 500 PRINT#2," PROPIEDADES DE LOS ELEMENTOS FINITOS:" 510 PRINT#2," ELEMENTO NUDO 1 NUDO 2 NUDO 3 NUDO 4 MODULO MODULO ESPESOR" 520 PRINT#2," ELASTICO POISSON" 530 INPUT#1,A$,A$,A$,A$,A$,A$,A$,A$,A$ 540 REM LEE PROPIEDADES DE LOS ELEMENTOS FINITOS 550 REM # ELEMENTO, NUDO 1, NUDO 2, NUDO 3, NUDO 4, MODULO ELASTICO, MODULO DE POISSON, ESPESOR 560 FOR I=1 TO NFINITOS 570 INPUT#1,J,NUDO1F(J),NUDO2F(J),NUDO3F(J),NUDO4F(J),EF#(J),POISSON#(J ),ESPESOR#(J) 580 PRINT#2,USING A8$;J,NUDO1F(J),NUDO2F(J),NUDO3F(J),NUDO4F(J),EF#(J),POISSON#(J),ES PESOR#(J) 590 NEXT I 600 REM CALCULA EL NUMERO DE GRADOS DE LIBERTAD Y ORDENA LAS ECUACIONES 610 NGRADOS=0 620 FOR I=1 TO NNUDOS 630 FOR J=1 TO 3 640 IF ORDEN(I,J)=0 GOTO 670 650 ORDEN(I,J)=0 660 GOTO 690 670 NGRADOS=NGRADOS+1 680 ORDEN(I,J)=NGRADOS 690 NEXT J 700 NEXT I 710 PRINT#2,"NUMERO DE GRADOS DE LIBERTAD =";NGRADOS 720 CLS:LOCATE 10,1:PRINT "NUMERO DE GRADOS DE LIBERTAD =";NGRADOS 730 REM CALCULA EL NUMERO DE ELEMENTOS POR COLUMNA MATRICIAL Y DETERMINA PUNTEROS DEL VECTOR SKYLINE 740 FOR I=1 TO NGRADOS 750 PUN(I)=1 760 NEXT I 770 REM DETERMINA VECTOR DE PUNTEROS EN FUNCION DE LOS ELEMENTOS FINITOS 780 FOR I=1 TO NFINITOS 790 FOR J=1 TO 3 800 IT(J)=ORDEN(NUDO1F(I),J) 810 IT(J+3)=ORDEN(NUDO2F(I),J) 820 IT(J+6)=ORDEN(NUDO3F(I),J) 830 IT(J+9)=ORDEN(NUDO4F(I),J) 840 NEXT J 850 REM ORDENA DE MENOR A MAYOR LOS GRADOS DE LIBERTAD DEL ELEMENTO FINITO 860 FOR J=2 TO 12

40

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador 870 FOR K=1 TO J-1 880 IF IT(J)>IT(K) GOTO 920 890 TEMP=IT(J) 900 IT(J)=IT(K) 910 IT(K)=TEMP 920 NEXT K 930 NEXT J 940 REM CALCULA LONGITUDES DE COLUMNAS MATRICIALES DEL ELEMENTO FINITO 950 FOR J=2 TO 12 960 IF IT(J)=0 GOTO 1020 970 FOR K=1 TO J-1 980 IF IT(K)=0 GOTO 1010 990 IF IT(J)-IT(K)+1IT(K) OR IT(J)=0 OR IT(K)=0 GOTO 1330 1270 TEMP=PUN(IT(K))+IT(J)-IT(K) 1280 GET #3,TEMP 1290 C1#=CVD(C$) 1300 C1#=C1#+KM#(J,K) 1310 LSET C$=MKD$(C1#) 1320 PUT #3,TEMP 1330 NEXT K 1340 NEXT J 1350 NEXT ELEMENTO 1360 REM OPERA LA MATRIZ DE COEFICIENTES CON LA TECNICA DEL SKYLINE 1370 CLS:LOCATE 10,1:PRINT "OPERA LA MATRIZ DE COEFICIENTES" 1380 EPSILON#=1E-10 1390 GET #3,1 1400 C1#=CVD(C$)

41

EL MÉTODO DE LOS ELEMENTOS FINITOS EN EL ANÁLISIS ESTRUCTURAL DE PLACAS Marcelo Romo Proaño, M.Sc. Escuela Politécnica del Ejército - Ecuador 1410 IF ABS(C1#)