Desarrollo de una herramienta numérica de análisis para losas de hormigón armado sometidas a aceleraciones verticales sísmicas

Development of a numerical tool for the analysis of reinforced concrete slabs under seismic vertical accelerations

Deni Herrera

Constructora El Sauce, Av. General Borgoño 934, oficina 301, Antofagasta, Chile, dherrera@elsauce.cl

Durante el terremoto en el centro sur de Chile del 27 de febrero del 2010 muchas estructuras fueron solicitadas fuertemente debido a aceleraciones horizontales considerables, sin embargo, la componente vertical de aceleración también fue importante. Esta característica ha motivado este trabajo, donde se estudia el efecto que produce la componente vertical de aceleración sobre losas de hormigón armado, teniendo en cuenta que la aceleración se aplica sobre el borde de la losa. La metodología utilizada contempló la elaboración de un programa computacional escrito en lenguaje MATLAB. Este programa utiliza el MEF y permite calcular tanto las tensiones como deformaciones en la losa. Inicialmente se realizó un análisis estático consistente en obtener las distribuciones de tensiones que genera una carga estática uniformemente distribuida de acuerdo a la norma NCh 1537.OJ86. Luego, se realizó un análisis sísmico consistente en obtener las envolventes de tensiones que genera la aplicación del registro vertical de aceleraciones de San Pedro de la Paz del terremoto del 27/02/2010. Se presenta el contexto teórico, la formulación del MEF y la validación numérica de los resultados. Posteriormente se analizan 2 losas de distinta geometría con el propósito de obtener las distribuciones de tensiones para el caso estático y envolventes de tensiones para el caso sísmico; al igual que las frecuencias y periodos naturales de cada losa. Los resultados revelan que las tensiones inducidas por las cargas estáticas son mayores a las inducidas por las cargas sísmicas para la losa con un periodo fundamental de 0.031 s. Mientras que para la losa con un periodo fundamental de 0.104 s resultó lo contrario.

During the Chile Centre-South 27^th February 2010 earthquake, many structures were strongly shaken due to considerable horizontal accelerations; however, the vertical acceleration component was important too. This feature has motivated this work, where the vertical acceleration component effect on reinforced concrete slabs is studied, when this acceleration acts on the slab edges. The methodology contemplated a computer software written in MATLAB code. This program uses the FEM and allows the calculation of stresses and strains in the slab. Initially a static analysis was carried out, which consisted of calculating the stress distributions generated by a uniformly distributed static load according to the norm NCh 1537Of.86. Subsequently, the seismic analysis consisted of obtaining the stress envelopes caused by the vertical acceleration recording of San Pedro de la Paz for the 27/02/2010 earthquake. The theoretical framework, the FEM formulation and the numerical validation of the results are presented. Then, two slabs with different geometry are analysed with the purpose of estimating the stress distributions for the static case and the stress envelopes for the seismic case; as well as the natural frequencies and periods for each slab. The results reveal that stresses induced by the static loads are higher than the induced by the seismic loading for the slab with a fundamental period of 0.031 s, whilst the contrary resulted for the slab with a fundamental period of 0.104 s.

Keywords: vertical accelerations, MZC and DKT finite elements, reinforced concrete slabs, plates, discretization

Introducción

Debido a la magnitud momento de 8.8 del terremoto que sacudió la zona centro sur de Chile el 27 de febrero del año 2010, muchas estructuras de hormigón armado sufrieron deformaciones que sobrepasaron el límite elástico de varios de sus elementos estructurales. Deformaciones irreversibles se evidenciaron por medio de la aparición de fisuras, grietas e incluso desprendimiento del hormigón y falla del acero de refuerzo. Esto condujo a la falla de importantes elementos estructurales e incluso el colapso de algunas estructuras (Betanzo 2010). En general las fallas más importantes suceden en muros, vigas y columnas de hormigón armado, en cambio en las losas los daños son menores. Tal vez debido a esto las losas de hormigón armado son menos estudiadas. Sin embargo, las fotos de la Figura 1 muestran ejemplos de daños en losas ocurridas en el terremoto del 27/02/2010. Por lo tanto, también se hace necesario estudiar el comportamiento sísmico de las losas y en particular el efecto de la componente vertical del movimiento, el cual es aún menos estudiado.

El terremoto del 27/02/2010 se caracterizó por presentar altas aceleraciones verticales, alcanzando en algunos registros valores muy similares a las aceleraciones horizontales (Barrientos 2010; Boroschek et al. 2010). Sin embargo, la componente vertical de la aceleración es raramente considerada en el diseño sismorresistente de edificios y menos aún en el diseño sismorresistente de losas de hormigón armado.

Ecuación diferencial de placas

La expresión que permite estimar la deformación vertical de una placa y que puede ser usada para calcular la flecha de una losa de hormigón a partir de las propiedades geométricas y resistentes de este material y las cargas aplicadas viene dada por (Timoshenko, 1937; Chakraverty 2009):

La expresión (1) es una ecuación diferencial de cuarto orden que relaciona la flecha w con la carga repartida q y las propiedades del material. Donde w es el desplazamiento vertical o flecha, es la masa por unidad de superficie y q es la carga distribuida sobre la losa. El parámetro representa la rigidez flexural de la placa y su expresión está dada por:

donde E es el módulo de Young o de elasticidad del material, ν es la relación de Poisson y t es el espesor de la losa.

Discretizacion de la losa por medio del MEF

Para resolver (1) se ha utilizado el Método de los Elementos Finitos MEF. Para ello se han escogido dos elementos en particular: el elemento rectangular MZC (Melosh 1963, Zienkievicz y Cheung 1967) y el elemento triangular DKT (Discrete Kirchhoff Triangular, ver Batoz et al. 1980). La diferencia que existe entre estos dos elementos, es que el elemento triangular proporciona más versatilidad geométrica que el elemento rectangular, por ejemplo se pueden realizar mallados con geometrías de losas circulares.

Al discretizar la losa y al aplicar el MEF sobre (1), se obtienen matrices de rigidez y masa, tanto para el elemento MZC como para el elemento DKT. Una vez que han sido calculadas las matrices de rigidez y masa tanto para el elemento MZC como para el elemento DKT, la ecuación que permite encontrar los desplazamientos de la losa es,

La expresión (3) es una ecuación de movimiento para sistemas de múltiples grados de libertad. Esta ecuación matricial de movimiento posee la característica de involucrar las matrices de masa [M], amortiguación [C] y rigidez [K] de la losa. Cabe destacar que los coeficientes que componen la matriz de amortiguación no pueden ser obtenidos explícitamente, debido a esto se incorporan en las ecuaciones de movimiento solo a nivel modal (Paz 1992). También se debe mencionar que (3) es un sistema de ecuaciones acophdo, por lo que para resolverlo se hace necesario desacoplarloa través del Método de Superposición Monal (Paz 1992). Al aplicar el Método de Superposición Modal sobre (3) se obtiene la siguiente ecuación modal de movimiento,

La expresión (4) ha sido resuelta por el método Newmark (1959), donde el lado derecho corresponde al registro verticañ de aceleraciones del terremoto 27/02/2010 obtenido en el Colegio Concepción de San Pedro de la Paz (Barrientos 2010).

Matriz de rigidez y masa del elemento MZC

La Figura 2 muestra el elemento rectangular de cuatro nodos denominado MZC, el cual se basa en la siguiente aproximación,

La aproximación (5) ha sido resuelta por Melosh (1963) y Zienkiewicz y Cheung (1967), entregando la siguiente expresión para la flecha w,

donde N son las matrices de funciones de forma y a es el vector de movimiento del elemento e de un nodo i, respectivamente. Las expresiones analíticas de las funciones de forma en coordenadas naturales son,

La matriz de deformación generalizada de flexión se obtiene a partir del vector de deformaciones generalizadas de flexión de la siguiente manera,

con

y

donde

es la matriz constitutiva de flexión para un material elástico isótropo. La matriz de masa elemental es,

La expresión del vector de fuerzas nodales equivalentes para una carga q uniformemente distribuida sobre el elemento e es:

Las integrales (13) y (15) se calculan por el método de integración numérica de Gauss para dos dimensiones (Oñate 1992), obteniendo así la matriz de rigidez y de masa elemental. Las matrices en coordenadas globales de rigidez y de masa, al igual que el vector en coordenadas globales de fuerzas nodales, se obtienen por medio de un proceso de ensamblaje a través de los vectores de coordenadas y conectividades de cada elemento que compone la malla de elementos rectangulares. Los esfuerzos en puntos nodales se calculan con la siguiente expresión,

Matriz de rigidez y masa del elemento DKT

La formulación del elemento DKT se basa en la teoría discreta de Kirchhoff para placas delgadas, la que da origen al elemento triangular discreto de Kirchhoff (Figura 3). El elemento DKT ha demostrado tener un excelente comportamiento con respecto a su convergencia (Batoz et al. 1980).

Para placas delgadas la deformación transversal por corte y, por lo tanto su energía de deformación U_s, es despreciable comparado con la energía de deformación por flexión U. Teniendo en cuenta lo anterior, la formulación de la matriz de rigidez del elemento DKT se basa en la siguiente expresión,

donde U es la energía de deformación por flexión, A es (17) el área del plano medio del elemento, ε_f es el vector de deformaciones por flexión y es la matriz constitutiva. No se considera la deformación por corte ya que para placas delgadas el aporte energético de la deformación por corte transversal es despreciable. Además (18) contiene solamente las primeras derivadas de θ_x y θ_y por lo tanto es posible establecer funciones de interpolación que satisfagan los requerimientos de compatibilidad, relacionando la rotación de la normal (θ_x y θ_y) a la superficie media con los desplazamientos transversales w. Este objetivo se alcanza con las siguientes consideraciones:

a.	El elemento triangular debe tener solo 9 grados de libertad, los cuales son el desplazamiento w y los giros θx y θy y de los tres nodos esquineros.
b.	Los giros nodales deben ser , de modo que se satisfagan las condiciones de borde del elemento.
c.	Ya que los modelos de elementos de placas delgadas son gobernados por la teoría de Kirchhoff, ésta se puede aplicar sobre cualquier punto discreto del elemento.
d.	La compatibilidad de los giros θx y θy no debe perderse.

El punto de partida para obtener la matriz de rigidez del elemento de 3 nodos DKT, es el elemento triangular de placa de Reissner-Mindlin de 6 nodos mostrado en la Figura 4, sometido a varias condiciones y/o consideraciones (Batoz et al. 1980).

Utilizando las consideraciones anteriores se puede establecer una interpolación del desplazamiento w (flecha) y los giros nodales θ_x y θ_y a través de,

donde u es el vector de movimientos y a es el vector de movimientos de cada grado de libertad del elemento,

N es la matriz de funciones de forma del elemento DKT que aproximan o interpolan los desplazamientos y giros nodales. La matriz de deformación generalizada de flexión se obtiene a partir del vector ε_f y (21) como,

donde B_f es la matriz de deformación generalizada del elemento definida por:

y para la matriz de masa se tiene de acuerdo a Luo y Hutton (2002) que,

Finalmente, la expresión del vector de fuerzas nodales equivalente para una carga q uniformemente distribuida sobre el elemento es,

donde q es la carga repartida sobre el elemento triangular y A es el área del elemento triangular. Al igual que para el elemento rectangular, las ecuaciones (25) y (26) se calculan por el método de integración numérica de Gauss para dos dimensiones (Oñate 1992), obteniendo así la matriz de rigidez y de masa elemental. Las matrices en coordenadas globales de rigidez y de masa, al igual que el vector en coordenadas globales de fuerzas nodales, se obtienen con el proceso de ensamblaje a través de la matriz de transformación de coordenadas, y también, por los vectores de coordenadas y conectividades de cada elemento que compone la malla de elementos triangulares. Los esfuerzos en los puntos nodales se calculan con la siguiente expresión,

Nuevamente se tiene que las tensiones calculadas en los nodos de cada elemento puede causar que en nodos comunes sus valores varíen. Esto se puede solucionar usando alisado de tensiones. Como también se dijo anteriormente, una vez que se ha resuelto (1) y se han utilizado (3) y (4) se pueden obtener las tensiones sobre la losa a través de (28). Cabe destacar que (28) se ha utilizado para el cálculo de tensiones en el caso estático y sísmico del elemento triangular.

Validación numérica

Validación del caso estático

En este caso se estudia la convergencia de la flecha central de una losa cuadrada con todos sus bordes empotrados. La solución analítica del problema, según la teoría de Kirchhoff para placas delgadas, está dada por la siguiente expresión (Novoa 2005),

donde w es la flecha central, q es la carga repartida, L es el largo, E es el módulo de elasticidad estático, ν es la razón de Poisson y t es el espesor de la losa. Las propiedades de la losa con las que se ha realizado el análisis y la comparación de los resultados se muestran la Tabla 1 y en la Figura 5.

De la Figura 5 se puede deducir que tanto el elemento MZC como DKT entregan aproximaciones razonables de la solución analítica de la flecha central de la losa. En las Figuras 6 y 7 se muestran los resultados de las deformadas de la losa con el elemento MZC ν DKT.

Validación de caso dinámico

La validación del caso dinámico se ha llevado acabo a través de la comparación de las 6 primeras frecuencias naturales obtenidas de la teoría de vibración libre de placas y las obtenidas a través del programa. Cabe destacar que la losa que se ha utilizado para la validación es una losa cuadrada simplemente apoyada. Existe una expresión teórica para las frecuencias naturales válida para estos tipos de losas dada por (Timoshenko, 1937; Chakraverty 2009):

donde es la rigidez flexural de una placa, ρ es la densidad del material de la placa, t es el espesor de la placa, a y b son las longitudes de los lados de la placa, y finalmente, m indica la numeración de los modos de vibrar de la placa. La Tabla 1 presenta los valores de las propiedades de la losa usados en los cálculos numéricos.

Las Tablas 2, 3 y 4 presentan los resultados obtenidos de las frecuencias naturales, en donde la letra m indica el número de frecuencia natural calculada.

]]>

De las Tablas 2, 3 y 4 es posible apreciar que los valores de las frecuencias naturales de las placas calculadas con el elemento MZC y DKT son aproximaciones razonables del valor teórico otorgado por (30).

Validación de la integración en el tiempo del registro de aceleraciones

Se ha analizado la respuesta de un oscilador de 1 grado de libertad por medio de la solución exacta y posteriormente, mediante la implementación del método numérico de Newmark (1959). El oscilador que se ha considerado posee las propiedades y valores mostrados en la Tabla 1.

La solución analítica de este oscilador está dada por la siguiente expresión:

donde ω₀ es la frecuencia de la fuerza externa. Cabe destacar que esta expresión corresponde a la respuesta de un oscilador de 1 grado de libertad sometido a una carga sinusoidal monotónica, la cual no se puede extender a un registro aleatorio como lo es un registro sísmico. Finalmente, en la Figura 8 se muestra la comparación de ambas respuestas, en donde se aprecia que el resultado numérico es prácticamente igual a la solución exacta.

Análisis comparativo

Para la comparación de ambos casos en el análisis sísmico de la losa, se ha considerado el punto 5.5.1 de la norma NCh 433Of96, la cual establece que para el cálculo de masas se deben considerar las cargas permanentes más un porcentaje de la sobrecarga de uso, que no podrá ser inferior a 25 % en construcciones destinadas a la habitación o al uso público. La Tabla 6 resume las características del material y condiciones de apoyo de la losa. Los valores adoptados permanecen constantes durante el análisis. Cabe destacar que se consideró en el análisis el módulo de elasticidad estático del hormigón. Las características geométricas de las losas que se han analizado se muestran en las Figuras 9 y 10. La Tabla 7 entrega las cargas estáticas uniformemente repartida sobre la losa, sobrecarga, peso propio y carga total.

Se debe tener en cuenta que en este trabajo, la masa de la losa corresponde a su masa como tal más un 100 % de la sobrecarga de uso (según NCh433Of96), la cual es de 2 kPa. En este caso se ha utilizado el 100 % del valor de la sobrecarga debido a que es el caso más desfavorable.

Las Figuras 12, 13 y 14 muestran las distribuciones de tensiones a partir del elemento MZC y DKT para la losa 1 bajo carga estática.

Las Figuras 15, 16 y 17 muestran las envolventes de tensiones a partir del elemento MZC y DKT para la losa 1 bajo carga sísmica. Notar que la escala en la dimensión horizontal χ ha sido reducida y que se han separado las tensiones positivas y negativas.

Las Figuras 18, 19 y 20 muestran las distribuciones de tensiones a partir del elemento MZC y DKT para la losa 2 bajo carga estática.

Las Figuras 21, 22 y 23 muestran las envolventes de tensiones a partir del elemento MZC y DKT para la losa 2 bajo carga sísmica. Notar que nuevamente se ha cambiado la escala y se han incluído las tensiones negativas.

En este caso tanto las tensiones máximas positivas como negativas que produce la aplicación del registro vertical del terremoto son siempre mayores a las que produce la aplicación de una carga estática vertical. Excepto en el caso de σ_x positivo con elementos triangulares.

Los resultados con tensiones negativas de -4 MPa podrían explicar las grietas y desprendimientos de hormigón especialmente en las uniones de las losas con los muros o vigas tal como se mostró en la Figura 1.

Las periodos y frecuencias fundamentales calculados para la losa 1 y 2 son 0.031 s (198.9 rad/s) y 0.104 s (60.4 rad/s) respectivamente. De esta manera el periodo fundamental de la losa 2 conlleva a aceleraciones espectrales mayores que la losa 1. Esto es porque la losa 2 es más flexible, ya que posee un periodo fundamental mayor. Por lo tanto, este aumento en la aceleración espectral se debe a que este terremoto en particular (en ese rango de periodos la aceleración espectral es ascendente), afecta mayormente a la losa 2 que a la losa 1. Esto explica que las tensiones provocadas por el registro vertical de aceleraciones aplicadas sobre la losa 2 superen a las que origina la carga estática aplicada sobre la misma losa, y para la losa 1 resulten menores.

Se han utilizado las teorías de placas basadas en la elasticidad lineal y la del movimiento de estructuras de 1 GDL para caracterizar el comportamiento sísmico de una losa de hormigón armado. Para ello se ha implementado un programa computacional en código MATLAB, el cual permite obtener desplazamientos verticales y tensiones ±σ_x, ±σy y ±τxy, producto de la aplicación de cargas estáticas y sísmicas. Se utilizaron los elementos finitos MZC y DKT, encontrándose no mayor diferencia en los resultados obtenidos al usar uno u otro elemento.

Las tensiones máximas en losas cuadradas y rectangulares empotradas surgen en los bordes de éstas, tanto para el caso estático como para el sísmico. La forma general de las distribuciones de tensiones en losas cuadradas y rectangulares son muy similares entre si, tanto para el caso estático como para el sísmico. Por otro lado, las fibras de una losa que habitualmente se encuentran en compresión (o tracción) bajo carga estática vertical, cambian temporalmente a tracción y compresión (compresión y tracción) durante la aplicación del registro de aceleraciones vertical del terremoto.

Los registros que posean una componente de aceleración vertical importante, como el del 27/02/2010 (San Pedro de la Paz), pueden provocar tensiones superiores a las que provocan cargas estáticas, como por ejemplo las cargas de diseño de la norma NCh 1537Of.86. Es por ello que se debería considerar la componente vertical en los análisis sísmicos, ya que sismos con aceleraciones verticales mayores a la del terremoto del 27/02/2010 podrían, eventualmente, provocar tensiones superiores a las que se consideran en los análisis para el diseño, tal como ocurre en las losas analizadas en este trabajo. Se recomienda en futuras investigaciones incluir a los muros en el análisis de aceleración vertical para así evaluar el efecto de interacción losa muro.

Referencias

Barrientos, S. (2010). Terremoto de Cauquenes 27 de Febrero del 2010. Informe técnico actualizado 27 de Mayo del 2010. Servicio Sismológico, Universidad de Chile.         [  ]

Batoz, J.L., Bathe, K.J. and Ho, L. W (1980). A study of three-node triangular plate bending elements. International Journal for Numerical Methods in Engineering 15, 1771-1812        [  ]

Betanzo, R. (2010). Daños estructurales y lecciones del terremoto del 27/f en el Gran Concepción. Obras y Proyectos 8, 59-75        [  ]

Boroschek, R., Soto, P. and León, R. (2010). Maule Region Earthquake. February 27, 2010, M_w = 8.8. Renadic report 10/08 Rev. 2. Universidad de Chile        [  ]

Chakraverty, S. (2009). Vibration of plates. Crc Press Taylor & Francis Group. USA        [  ]

IDIEM (2010). Informe Inspección Post Sismo del 27 de Febrero del 2010, Edificio Don Manuel. Informe ejecutivo N° 595408        [  ]

Luo, Z. and Hutton, S.G. (2002). Formulation of a three-node traveling triangular plate element subjected to gyroscopy and in-plane forces. Computers & Structures 80, 1935-1944.         [  ]

Melosh, R.J. (1963). Basis for derivation of matrices for the direct stiffness method. Journal of the American Institute of Aeronautics and Astronautic AIAA 1, no. 7, 1631 - 1637        [  ]

NCh 433Of.96. Diseño sísmico de edificios. Instituto Nacional de Normalización, Chile        [  ]

NCh 1537Of.86. Diseño estructural de edificios, cargas permanentes y sobrecargas de uso. Instituto Nacional de Normalización, Chile        [  ]

Newmark, N.M. (1959). A method of computation of structural dynamics. Proceedings American Society of Civil Engineers 85, EM3, 67-94        [  ]

Novoa, N. (2001). Análisis de cargas críticas en estructuras de cáscara delgada. Memoria para optar al título de Ingeniero Civil. Universidad de Concepción        [  ]

Oñate, E. (1992). Cálculo de Estructuras por el Método de Elementos Finitos. Centro Internacional de Métodos Numéricos en Ingeniería, Primera Edición. España        [  ]

Paz, M. (1992). Dinámica estructural, teoría y cálculo. Editorial Reverté. Tercera edición, España        [  ]

Timoshenko, S. (1937). Vibration problems in engineering. Vin Nostrand Company, Inc., USA        [  ]

Zienkiewicz, 0. C. and Cheung, Y K. (1967). The finite element method in structural and continuum mechanics. McGraw-Hill        [  ]

Fecha de entrega: 25 de agosto 2011 . Fecha de aceptación: 22 de noviembre 2011