A. Ferriz, O. Fruitos, M. Coma, S. Oller
Centre Internacional de Metodes Numerics a l'Enginyeria - CIMNE, Barcelona, Spain
En este informe se muestran distintas simulaciones realizadas sobre dos tipos de modelos: sólo el domo (caso A) y el conjunto domo y cilindro (caso B). Se pretenden determinar los niveles de presión interior a partir de los cuales aparece daño en la matriz y/o fibra en cada caso. Se evalúa no solo su localización inicial sino también su evolución al aplicar presión interior de forma progresiva hasta el valor nominal (22 MPa) y hasta una sobrecarga de 2,4 veces la presión nominal (52,8 MPa).
Se considerarán distintos espesores: 10 y 32 mm para el domo; 30 y 95 mm para el cilindro, ambos definidos por ACCIONA en los documentos: Tank.pdf [1] y informe simulación031008-2.pdf [2]. Adicionalmente se incluyen resultados preliminares para el modelo B con domo de espesor y dirección de fibras variables.
En ninguno de los casos se ha considerado el liner de aleación de aluminio.
Para realizar este estudio se ha escogido en todos los casos la geometría correspondiente al plano medio. Se ha incluido la tapa central del domo de diámetro 110 mm, y las condiciones de contorno estudiadas analíticamente en [3]. Los modelos geométricos estudiados son los siguientes:
Otras condiciones de los modelos se describen a continuación:
En las figuras 1 y 2 se muestran las mallas usadas para la simulación:
Estos modelos han sido realizados completamente con elementos de lámina BST, por lo que la geometría dibujada corresponde al plano medio de la geometría real.
Las propiedades de los materiales informados son las definidas en el documento [3] y que se muestran sintéticamente en la tabla 1:
Fibras de Carbono | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Módulo de Young | 242 | GPa | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Módulo de Poisson | 0.2 | - | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Tensión de rotura | 3800 | MPa | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Densidad | 1810 | Kg/m3 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Matriz Polimérica | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Módulo de Young | 4 | GPa | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Módulo de Poisson | 0.353 | - | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Tensión de rotura | 59 | MPa | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Densidad | 1496 | Kg/m3
|}
Tabla 2. Propiedades de los compuestos
Descripción de las simulacionesPara tener una visión global de las simulaciones que se describen en este documento se destacan sus principales diferencias en la tabla 2:
Tabla 3. Características principales de los modelos estudiados
Condiciones de cargaLa condición de carga informada en todas las simulaciones es una presión interior en función del tiempo según la curva de la figura siguiente:
Figura 3. Gráfico de la evolución de la presión interna aplicada en función del tiempo
Resultados obtenidos en las simulacionesConsideraciones generalesEn este apartado se muestran imágenes del daño en la matriz y/o fibra en el instante en que este se inicia así como su evolución en función del tiempo para las distintas simulaciones en nodos de la zona señalada en la figura 4. Se escoge esta zona de estudio por estar alejada de las singularidades (unión del cilindro y el domo en el caso B y la tapa en los casos A y B).
Figura 4. Zona para toma de datos de daño (marcada en blanco) y nodos donde se obtienen los esfuerzos de sección
Los valores de referencia de los esfuerzos de sección Nxx (según los paralelos del domo) y Nyy (según los meridianos del domo) son los calculados analíticamente para paredes delgadas de material homogéneo lineal. Las expresiones de cálculo para estos valores de referencia son las siguientes: Donde: R es el radio interior del domo o del cilindro Pint es la presión interior aplicada que varía en función del tiempo
Modelo se03-f2 (Caso A, espesor domo 10 mm)A continuación se muestran resultados del daño obtenido en la matriz polimérica (figura 5) y en las fibras de carbono (figura 6) para el instante de tiempo en que aparece daño en el anillo intermedio del domo (zona marrón de toma de resultados de figura 4). Dicho instante corresponde a una presión de 7,04 MPa, que es un 32% de la presión nominal (ver figura 3).
Figura 5. Daño en matriz para 7,04 MPa, 32% de la presión nominal
Tal y como se observa en la figura de daño en fibras (fig. 6), para este instante de tiempo el valor de daño es zero en todo el domo, lo que significa que no han empezado a dañarse. Cuando el valor de daño es 1 significa que está completamente dañado en toda su sección.
Figura 6. Daño en fibra para 7,04 MPa, 32% de la presión nominal
A continuación se muestran las deformaciones principales a lo largo de toda la geometría (figuras 7 y 8). Como a la tapa fictícia se se ha dado un espesor diez veces mayor, la deformación es mínima en esta zona.
Figura 7. Deformación E1 para 7,04 MPa, 32% de la presión nominal
Figura 8. Deformación E2 para 7,04 MPa, 32% de la presión nominal
En la siguiente figura se puede ver la deformada de los elementos BST (elementos shell o de lámina que se han empleado en las distintas simulaciones) para una presión muy próxima a la nominal (92%).
Figura 9. Deformada para 20,24 MPa, 92% de la presión nominal (22 MPa)
A continuación se puede ver como evoluciona el daño en la matriz en los puntos representativos indicados en la figura 4. La evolución del daño se ha realizado en función de la presión del cilindro (fig. 10) y también a lo largo del tiempo (fig. 11) en que se ha aplicado en la simulación, descrito en la figura 3.
Figura 10. Evolución del daño en función de la presión de los nodos indicados en la figura 4
Figura 11. Evolución del daño de la matriz en función del tiempo en nodos seleccionados de la zona marcada en la figura 4
Las figuras que vienen a continuación describen los esfuerzos de sección en las dos direcciones principales Nxx y Nyy (al tratarse de elementos de lámina la Nzz no se obtiene). Los valores se han obtenido en los 3 puntos que siguen una línea meridional en el domo y en el cilindro.
Figura 12. Gráfica de Nxx vs. presión y valor de referencia para el domo según la expresión (1)
Figura 13. Gráfica de Nxx vs. tiempo y valor de referencia para el domo según la expresión (1)
Figura 14. Gráfica de Nyy vs. presión y valor de referencia para el domo según la expresión (2)
Figura 15. Gráfica de Nyy vs. tiempo y valor de referencia para el domo según la expresión (2)
En la siguiente figura aparece el desplazamiento que se produce en la tapa en el eje z, correspondiente al eje longitudinal del depósito. Los valores son referidos a la presión a la que se producen. Se puede observar que los valores son positivos, cosa que indica que esta parte de la geometría se desplaza hacia dentro del depósito.
Figura 16. Evolución del desplazamiento axial de la tapa con la presión
Modelo se03-k (Caso A, espesor domo 32 mm)El segundo modelo simulado empieza a dañar en la matriz a una presión que duplica el valor de la observada en el modelo anterior. Esto se ve reflejado a continuación en las figuras 17 y 18 que contienen los outputs de daño en matriz y fibras respectivamente. La presión para la cual aparece daño en la zona central (criterio mencionado ya en la anterior simulación) corresponde a una presión de 12,36 MPa, que es un 56% de la presión nominal.
Figura 17. Daño en matriz para 12,36 MPa, 56% de la presión nominal
En este caso (fig. 18) tampoco se produce daño en las fibras en este instante de tiempo crítico.
Figura 18. Daño en fibra para 12,36 MPa, 56% de la presión nominal
A continuación se muestran las deformaciones principales a lo largo de toda la geometría (figuras 19 y 20) para la presión en que aparece el daño. Como a la tapa fictícia se se ha dado un espesor diez veces mayor, la deformación es mínima en esta zona. Las deformaciones mayores se producen en la zona donde en la realidad física se une el domo con el cilindro (al igual que en el caso del daño).
Figura 19. Deformación E1 para 12,36 MPa, 56% de la presión nominal
Figura 20. Deformación E2 para 12,36 MPa, 56% de la presión nominal
Figura 21. Módulo de la deformación para 22 MPa, presión nominal
Figura 22. Evolución del daño en función de la presión de los nodos indicados en figura 4
Figura 23. Evolución del daño de la matriz en función del tiempo en nodos seleccionados de la zona marcada en la figura 4
Las figuras que siguen describen los esfuerzos de sección en las dos direcciones principales Nxx y Nyy en los mismos 3 puntos de la anterior simulación. En las figuras 24 y 26 se refieren a la presión interior a la que está sometido el domo. En las figuras 25 y 27 se refieren al tiempo de simulación a la que se encuentra el domo (la correspondencia entre tiempo de simulación y presión se encuentra en la figura 3).
Figura 24. Gráfica de Nxx vs. presión y valor de referencia para el domo según la expresión (1)
En el caso del valor de referencia, no tiene en cuenta las orientaciones de fibra y, en el caso de una semiesfera, alcanza valores teóricos exactamente iguales para Nxx y Nyy. Es por este motivo que en el caso de la simulación mediante MEF (con orientaciones de fibra asignados), los valores Nxx quedan por debajo de los valores teóricos de referencia y en cambio las Nyy quedan por encima.
Figura 25. Gráfica de Nxx vs. tiempo y valor de referencia para el domo según expresión (1)
Figura 26. Gráfica de Nyy vs. presión y valor de referencia para el domo según la expresión (2)
Figura 27. Gráfica de Nyy vs. tiempo y valor de referencia para el domo según expresión (2)
En la siguiente figura aparece el desplazamiento que se produce en la tapa en el eje longitudinal del depósito. Los valores son también función de la presión. Se puede observar que los valores son positivos, cosa que indica que en este caso la tapa también se desplaza hacia dentro del depósito.
Figura 28. Evolución del desplazamiento axial de la tapa con la presión
|
Debido a que la simulación llega a daño antes de llegar al estado de carga nominal, se ha realizado un cálculo no definido en [3] para conocer cuál es el valor del incremento de radio, que en base a los cálculos analíticos, debería tener el domo.
Para el cálculo por simulación se han verificado los radios iniciales y finales de la geometría para el instante en que la carga llega a un 32% de la nominal (cuando la matriz empieza a dañar en la zona media del domo). Los resultados del incremento de diámetro de la simulación se muestran en la figura 35:
A partir del diámetro inicial y final obtenido en la simulación se obtiene el incremento de radio mediante:
Del mismo modo que en el cálculo analítico se obtiene el alargamiento anular en la simulación:
En este apartado se muestra el error relativo entre los cálculos analíticos y los de simulación para los alargamientos anulares en el caso A:
Para la predicción del daño en la matriz se compara la tensión de servicio en la matriz calculada analíticamente en [3] en el estado de servicio con el valor máximo que soporta la matriz (59MPa).
De este modo, considerando que el material se comporta de manera elástica y lineal para estados previos al daño, y sabiendo que la curva de presiones aplicada al depósito en la simulación es la que se muestra en la figura 3.
Se considera que el porcentaje de presión de servicio aplicable sin superar la tensión máxima admisible en la matriz sería:
Según esto se podría esperar un inicio del daño en la matriz para un 37,67% de la presión de servicio.
El cálculo por simulación se ha mostrado anteriormente en el modelo se03-f2, que es el que tiene las mismas condiciones que los cálculos analíticos mostrados en [3]. En este modelo, el daño se generaliza en la zona media de la tapa para un 32% de la presión nominal.
En este apartado se muestra el error relativo entre los cálculos analíticos y los de simulación para el porcentaje de carga nominal en el que se inicializa el fallo:
Se impone una condición de compatibilidad circunferencial entre el domo y el cilindro, teniendo en cuenta los estados de tensión que se desarrolla en un punto de la unión entre ambos. De esta condición de compatibilidad resultará una presión
que representa la interacción entre las dos estructuras, como consecuencia del contacto entre ambas (unión solidaria). Este cambio de presión actuará incrementando la presión en el cilindro (estructura más rígida) y disminuyendo la presión en el domo (estructura menos rígida). Esto es:
Siendo la condición de compatibilidad entre cilindro y domo, la siguiente:
Teniendo en cuenta que:
Y sustituyendo en la ec.(24) resulta,
Sustituyendo en esta ecuación
y
, resulta la siguiente magnitud para el incremento de presión
Sustituyendo en esta expresión las magnitudes geométricas y mecánicas allí indicadas (ver Anexo), resulta
.
Para la comparación de los resultados se ha encontrado el valor de la presión en la zona de unión entre cilindro y domo tanto para el cilindro como para el domo. Con este cálculo se pretende encontrar los valores de ∆p, λCsim y λTsim y compararlos con los valores analíticos.
Para el cálculo de la presión en cada una de las partes (domo y cilindro) se han graficado los esfuerzos de sección (ver figura 41) y se ha encontrado la presión interior en cada caso a partir de las expresiones:
A partir de los valores de Nyy que se muestran en la gráfica para la presión nominal (t=0,0025 segundos) se encuentran los valores de la presión interior del domo y del cilindro:
En este apartado se muestra el error relativo entre los cálculos analíticos y los de simulación para los factores de corrección:
Al igual que en el apartado anterior, aquí se estudia el estado de tensión en los materiales componentes del domo que están en contacto con el cilindro. El estado de deformación se obtiene a partir de las tensiones calculadas membranalmente (3.2 en [3]) y del tensor constitutivo del material compuesto para el cilindro (ver Anexo 4.3 en [3]).
Debido a que la simulación llega a daño antes de llegar al estado de carga nominal, se ha realizado un cálculo no definido en [3] para conocer cuál es el valor del incremento de radio, que en base a los cálculos analíticos, debería tener el domo.
Del mismo modo que en el caso anterior del domo se han realizado las verificaciones geométricas para la geometría con la unión. En este caso la deformada se ha tomado en un 28% de la carga nominal (cuando la matriz empieza a dañar en la zona media del domo). Los resultados del incremento de diámetro de la simulación se muestran en la figura 37:
A partir del diámetro inicial y final obtenido en la simulación se obtiene el incremento de radio mediante:
Del mismo modo que en el cálculo analítico se obtiene el alargamiento anular en la simulación:
En este apartado se muestra el error relativo entre los cálculos analíticos y los de simulación para los alargamientos anulares en el caso B:
En este apartado se muestran los resultados obtenidos por simulación para una propuesta de distribución de espesor y dirección de fibras variable en el domo calculada a partir de los criterios propuestos por ACCIONA en la referencia [4]. En las tablas 3 y 4 se muestran los ángulos de apilado y espesores aplicados en cada uno de los anillos indicados en la figura 4. Estos apilados y espesores están informados sobre el plano medio de la geometría correspondiente a los espesores 36,2 y 99,2 mm (ver figura 2), por lo que se asume provisionalmente un error mayor que en los casos de espesor constante en el domo.
Caso | Modelo geométrico | Simulación | Espesor cilindro (mm) | Espesor domo (mm) | Espesor tapa central (mm) |
B* | hid-15b.gid | hid-15q | 95 | Variable entre 35 y 120 | 320 |
*Depósito completo con simetría transversal en cilindro y unión solidaria entre domo y cilindro |
Coordenada axial aproximada (mm) | División en la geometría de simulación (ver figura 1) | Ángulo de apilado real aproximado (º) |
293 | D1 | ±10 |
250 | D2 | ±10 |
211 | D3 | ±12 |
171 | D4 | ±14 |
130 | D5 | ±16 |
90 | D6 | ±19 |
48 | D7 | ±31 |
25 | D8 | ±58 |
Capas en modelo | Radio (mm) | Thickness (mm) |
D7 y D8 | 1.38*-82 | 120 |
D6 | 82-123 | 98 |
D5 | 123-164 | 69.5 |
D4 | 164-205 | 54 |
D3 | 205-246 | 44 |
D2 | 246-287 | 37.5 |
D1 | 287-314.2 | 33.5 |
Tapa | 314.2 | 362** |
*El origen de coordenadas está situado en el polo del domo.
**Para aumentar la rigidez debida al conducto conectado al domo. |
En la figura 35 se muestran los espesores indicados en la tabla 5 y la curva de distribución correspondiente a la propuesta descrita en la referencia [4].
A continuación se muestran los resultados obtenidos en la simulación con la orientación de fibras y el espesor en cada anillo que se muestran en las tablas 3 y 4.
En este caso aparece daño a los 39,25 MPa de presión, que corresponden al 56% de la presión de explosión. El daño se localiza inicialmente en zona del domo que se une con el cilindro y en menor grado en las proximidades de la tapa. En toda la geometría no aparece daño en las fases iniciales de daño en la matriz. El cilindro no daña.
A continuación se muestra la deformada (factor 1) para 22 MPa de presión.
Las figuras 67 y 68 muestran la distribución de los esfuerzos de sección en las direcciones principales Nyy y Nxx. En el caso de los Nxx, los valores son positivos a lo largo de todo el modelo. En lo que se refiere a Nyy, los valores son positivos excepto en el domo, justo en la parte de unión con el cilindro.
Las figuras 69 y 70 muestran la distribución de las deformaciones principales E11 y E22. En el caso de los E11, los valores son positivos a lo largo de todo el modelo. En lo que se refiere a las E22, los valores son positivos excepto en el cilindro, justo en la parte de unión con el domo.
Debido a que en la zona media del domo la matriz no llega a dañar, se han tenido que escoger nodos distintos a los mostrados para el resto de simulaciones. En la figura 59 se muestran los nodos en los cuales se ha realizado la toma de datos para el daño en matriz:
Del estudio de la simulación MEF se puede concluir que:
[1] Documento Tank.pdf
[2] Informe proceso de simulación. ACCIONA (documento informe simulación031008-2.pdf
[3] Cálculo de un Depósito de Hidrógeno a Presión construido en Material Compuesto por el método “filament-winding” – Cálculo elástico analítico, Sergio Oller, Barcelona, 2009. (documento: Tanque Hidrogeno-T Mezclas - ANALITICO.pdf).
[4] Diseño preliminar tanque para almacenaje de hidrógeno. ACCIONA (documento Tanque_t3.pdf)
===
Anexos 1, 2 y 3===
Published on 01/01/2010
Licence: CC BY-NC-SA license
Are you one of the authors of this document?