Cálculo de un depósito de hidrógeno
Simulación realizada con FEMCOM 1.0. Modelos con elementos hexaedros
M. Coma, A. Férriz, O. Fruitos, X. Martinez
Centre Internacional de Metodes Numerics a l'Enginyeria - CIMNE, Barcelona, Spain
1 INTRODUCCIÓN
En el presente informe se exponen los resultados de cuatro modelos obtenidos con la versión de FEMCOM 1.0. En los dos casos primeros se han utilizando elementos hexaédricos. En el primer caso se estudia la geometría proporcionada por ACCIONA. El segundo caso de estudio es una variación del modelo anterior, también con elementos hexaédricos, añadiendo un refuerzo en el domo. En los dos últimos casos son simulaciones con elemento finito de lámina con espesor según la ley de ACCIONA y la segunda simulación con espesores incrementados.
Se ha incorporado el liner de aleación de aluminio de 4.2 mm de espesor, que no se había incluido en las simulaciones con láminas remitidas a ACCIONA en el informe [3].
Se ha pensado en utilizar elementos hexaédricos para obtener los resultados con un modelo que considere con mayor precisión la geometría en el domo, dada su complejidad.
2 OBJETIVOS
- Obtener los resultados con un modelo de elementos hexaédricos de la geometría propuesta en [1] y su comparación con elementos de lámina..
- Estudiar posibles mejoras en el diseño del tanque de hidrógeno.
3 SISTEMA DE COORDENADAS Y VARIABLES DEL PROBLEMA
Para aclarar la nomenclatura utilizada en todo el documento, a continuación se describen las variables del modelo. El centro del sistema de coordenadas se sitúa en el centro de la esfera inscrita en el domo. El eje X apunta al polo, centro del agujero del domo. El eje r está contenido en el plano perpendicular a X. El ángulo de elevación es el ángulo entre el vector que va del origen al punto genérico con su proyección sobre el plano perpendicular a X.
Las variables del problema son:
x Coordenada axial.
α Ángulo de bobinado (dirección de fibra) respecto al meridiano.
r Distancia entre la proyección del punto de la esfera sobre el plano de unión cilindro-domo y el centro de la esfera.
αc Ángulo de bobinado (dirección de fibra) respecto del meridiano en la unión cilindro-domo (ángulo de elevación = 0º)
thick Espesor del domo, medido en dirección radial.
R Diámetro exterior del liner de aluminio que se corresponde con el diámetro interior del material compuesto.
4 MODELO GEOMÉTRICO PROPORCIONADO POR ACCIONA
4.1 Definición de la geometría
La Figura 2 muestra las dimensiones básicas del tanque estudiado.
En [1] se especifica un diámetro del agujero de 40 mm. El modelo realizado se ha definido con un diámetro de 101,2 mm en el interior del conducto de aluminio, de tal forma que el diámetro del agujero del material compuesto en esa zona es de 110 mm.
A continuación se detallan las distintas zonas del depósito.
En la simulación se ha definido un octavo de la geometría aprovechando la simetría tanto del modelo como de la carga aplicada.
4.1.1 Domo
Se ha utilizado la distribución de espesor resultante de la expresión (1), ver [1]:
La expresión del thick en función de r tiene una asíntota vertical en el agujero. Para evitar un espesor infinito se ha modificado el criterio enunciado en [1] según el cual:
“en la abertura del domo, este espesor es corregido para ser constante desde el 98% de la longitud meridional desde la unión cilindro-domo”
No se ha aplicado este criterio porqué la geometría resultante era poco creíble. Se ha optado por mantener el espesor constante en el último 10% de la longitud del arco comprendido entre la unión domo-cilindro hasta el inicio del agujero del domo, que correspondería a un 90% y no el 98% indicado en [1].
Calculando el espesor con la expresión (1) para un radio interior R=314.2 mm, y asumiendo que el espesor es constante a partir de r=97.54 mm, correspondiente al 90% indicado en el párrafo anterior, se obtiene la distribución de espesor mostrada en el Gráfico 1.
El contorno exterior del domo se consigue transformando la curva del Gráfico 1 a coordenadas esféricas y por revolución de esta geometría se genera el volumen.
4.1.2 Cilindro
El cilindro no presenta variaciones de diseño respecto de la geometría descrita en [1].
La geometría completa del modelo, definida y visualizada en FEMCOM, se muestra en la Figura 3. El cilindro en el modelo de la simulación tiene una longitud de 1500 mm y un ángulo circunferencial de 90º debido a la simetría.
En la Figura 3 el cilindro se muestra dividido en dos volúmenes dado que son de distintos materiales compuestos (capa circunferencial y capa helicoidal).
4.2 Materiales
4.2.1 Material base
Las propiedades asignadas a la fibra de carbono del modelo numérico calculado con el software FEMCOM son las siguientes:
Los valores a los que se asignan valores de 0 y 0,0, se calculan de forma automática dentro de FEMCOM (Characteristic lenght) o bien sus valores en cero son los adecuados en este tipo de simulación.
Las propiedades de la matriz de epoxi se detallan en la siguiente figura:
En cuanto a las propiedades asignadas a la aleación de aluminio, componente que no se tuvo en cuenta en análisis anteriores a los expuestos en este informe pero que se ha considerado que puede tener relevancia, son las descritas en la Figura 6.
4.2.2 Compuestos
Una vez definidos los materiales base, se definen los materiales compuestos explicitando el porcentaje de cada uno en el compuesto. La fracción de fibra respecto del total es del 60%.
En lo referente a las orientaciones de fibra que se han aplicado a cada parte del tanque de hidrógeno, se describen a continuación.
En la parte exterior del cilindro, se han dispuesto fibras a +55º y -55º (esta área corresponde a la capa circunferencial mostrada en la Figura 2). En el interior del cilindro, las fibras estan a +10º y -10º. Éstas últimas fibras se corresponden con la capa helicoidal mostrada en la Figura 2. En este caso, la capa helicoidal tiene continuidad en el domo pero su orientación cambia con la posición.
En el Gráfico 2 se observa la variación de α en función de r, (ver expresión en [1]):
Se observa que en la coordenada r=55 mm el ángulo de bobinado respecto al meridiano es de 90º, dado que tiene que ser tangente al agujero del domo. Con un coeficiente de fricción de λ=0.355 se obtiene αc=-10º (ángulo de entrada al domo desde el cilindro r=314.2mm) como se observa en el Gráfico 2.
4.3 Cargas aplicadas y otras consideraciones
Las cargas aplicadas sobre la estructura son las siguientes:
:4.3.1 Presión interna=
La presión interna se ha informado sobre todo el interior del depósito, siguiendo una curva temporal según el Gráfico 3.
La presión de servicio es de 22MPa, ésta se mantiene entre los 0,0025s y los 0,005s. Luego se eleva la presión hasta 52,8MPa.
4.3.2 Efecto de la válvula en el agujero
Cómo la presión interna, solo actúa sobre la estructura definida en el modelo, se ha procedido a restringir el movimiento de los nodos del contorno del agujero en las direcciones radiales, solo se permite el movimiento en la dirección axial del depósito. De este modo, se asume una rigidez infinita, que no permite que el agujero se abra. Para añadir la presión sobre el agujero, ejercida por la válvula u otros elementos de cierre, se ha añadido la fuerza en el perímetro del agujero en la dirección X, aplicada en los nodos. Esta fuerza ha sido calculada a partir de la presión interna y del área del agujero.
4.4 Malla – Consideraciones sobre la discretización
Se utiliza una malla de elementos hexaédricos de 8 nodos estructurada. El número de divisiones (particiones) en la líneas utilizado se muestra en la Figura 8.
El número de divisiones en la dirección de los paralelos son 30. En la Figura 9 se muestran las distintas secciones que se han utilizado para informar la variación del ángulo de bobinado de las fibras en el domo en función del ángulo del elevación.
4.5 Resultados del modelo proporcionado por ACCIONA
En este apartado se muestran los resultados obtenidos en la simulación del depósito sin refuerzo. Tal y como se ha comentado en apartados anteriores, este modelo se ha informado con orientación y espesor variable.
La Figura 10 muestra el daño del compuesto para la presión de servicio:
En la Figura 11 se muestra el daño en la matriz polimérica y en la fibra de carbono. Tal y como se observa la fibra no está dañada.
 |
 |
Tal y como se puede observar en la Figura 12, el daño se inicia cuando la presión alcanza los 19,15 MPa (87,06% de la presión nominal). En este caso el daño se inicia en la parte del agujero.
En la Figura 13 se muestra el daño en el compuesto para 25,6 MPa:
En la Figura 14 se muestra el daño en la matriz polimérica y en la fibra de carbono. Tal y como se observa la fibra no está dañada.
 |
 |
La tensión a tracción para la presión de servicio se muestra en la Figura 15:
La tensión principal en la dirección 2 se muestra en la Figura 16:
La tensión principal en la dirección 3 se muestra en la Figura 17:
En la Figura 18 se muestra el daño medio de fibra y matriz para la presión de explosión:
5 PROPUESTA DE REFUERZO
Debido a que el depósito elaborado con elementos sólidos mostrado anteriormente daña antes de llegar a la presión de servicio, se ha decido informarle un refuerzo por tal de evitar este daño prematuro.
5.1 Modificación de la geometría
El refuerzo añadido consiste en una prolongación de la capa circunferencial del cilindro hacia el domo. Ésta se prolonga hasta 36º de elevación con espesor constante de 10,5 mm,
En la zona del refuerzo, midiendo el ángulo de la misma forma que en la capa helicoidal en el domo, el ángulo de bobinado respecto a los meridianos varía entre 55º en la unión cilindro-domo hasta 90º.
5.2 Materiales añadidos
Para este modelo no se ha añadido ningún material base debido a que se ha usado el mismo tipo de matriz y fibras que las ya informadas en el resto del modelo (tanto domo como cilindro).
5.3 Malla de la parte añadida
La malla de la parte añadida debe tener el mismo número de divisiones que la malla colindante del domo ya definido en el modelo anterior. Para el refuerzo se ha decidido asignar un elemento en el espesor del refuerzo. En este caso también ha sido definida mediante hexaedros.
5.4 Resultados del modelo con refuerzo
Los resultados para este caso se muestran a continuación.
El N Damage para la presión de servicio (22 MPa) se muestra en la Figura 21.
El inicio del daño se retrasa hasta la presión de 25,6 MPa (ver Figura 22). Para esta presión se inicial el daño en la zona de unión entre cilindro y domo pero en las capas que no se prolongan en el domo. Se observa además una fuerte retracción longitudinal de las capas exteriores del cilindro. En el detalle de la izquierda de la Figura 22 se observa la ausencia de daño en la cara interior del domo.
Las tensiones principales 1, 2 y 3 en el depósito para la presión de servicio se muestran desde la Figura 24 hasta la Figura 26.
La simulación se ha dejado continuar hasta la presión de explosión sin incrementarse el daño más allá del 20% de su valor máximo (ver Figura 27) y manteniéndose estable hasta la finalización del cálculo (0.01 segundos). En ningún momento se observa daño en el liner de aluminio o fibras de carbono (ver Figura 28).
Debido a que el daño en las fibras es de cero, todo el daño observado en el compuesto (ver Figura 27) reside en la matriz.
5.5
Discusión de los resultados
5.5.1 Modelo de sólidos sin refuerzo
Tal y como se puede ver en las imágenes anteriores el daño empieza antes de que se alcance la presión de servicio (22 MPa), aunque solo en la parte del agujero. Para la presión de servicio el daño está por encima del 8%, en cambio, la tensión a tracción es inferior a la máxima soportada por las fibras (378 MPa en frente de los 3800 MPa que soportan las fibras de carbono).
Finalmente, para la presión de explosión se puede observar que prácticamente todo el depósito está dañado al 100%, y que las fibras, debido a errores numéricos, están incluso por encima del valor límite real (mayor que el 200%).
5.5.2 Modelo de sólidos con refuerzo
Se observa que el daño es nulo. La presencia de refuerzo en la zona media del domo (donde la orientación de fibra es casi paralela a los meridianos) y hasta la unión con el cilindro ha contribuido a la desaparición del daño observado en los modelos sin refuerzo.
La retracción longitudinal de las capas exteriores del cilindro es producida por un efecto similar al efecto Poisson derivado de la expansión del cilindro.
Se puede concluir que el modelo con refuerzo es mucho más resistente que el que no dispone de refuerzo. Es por este motivo que, en base a los resultados simulados obtenidos, se recomienda la adición del refuerzo a la geometría original del depósito.
:6 MODELO DE LÁMINAS
6.1 Modelo hid15r
El modelo presentado es un modelo en el que se ha informado la variación de espesor según el Gráfico 1 .
La geometría considerada es la misma que en el modelo de sólidos. Ver apartado 4.1.
Los materiales base son los mismos que en los modelos de sólidos. Estos se pueden ver en el apartado 4.2.1.
Las diferencias de este modelo respecto al modelo de sólidos, residen principalmente en:
- Discretización del espesor en el domo. Ver Gráfico 1.
- Discretización de la orientación de fibras. Ver Gráfico 2.
- Modelado de la unión cilindro-domo.
6.1.1 Malla – Consideraciones sobre la discretización
Debido a que existe un cambio de espesor entre las dos estructuras involucradas (domo y cilindro) y a que en FEMCOM es necesario definir la superficie de la lámina en el plano medio del espesor. Entonces, los planos medios del domo y cilindro no son coincidentes y no transmiten el momento. Por lo tanto, el modelado de la unión entre el cilindro y el domo se ha realizado mediante una superficie auxiliar por tal de que se transmitan los momentos del cilindro al domo y viceversa.
Por otro lado, la discretización de la orientación de fibras en el domo es la siguiente:
Tabla 1. Distribución de espesores
| Capas en modelo | Radio (mm) | Thickness (mm) |
| D7 y D8 | 1.38*-82 | 119 |
| D6 | 82-123 | 76 |
| D5 | 123-164 | 67 |
| D4 | 164-205 | 62 |
| D3 | 205-246 | 57 |
| D2 | 246-287 | 52,35 |
| D1 | 287-314.2 | 35 |
| Tapa | 314.2 | 352** |
| Cilindro | 314.2 | 127 |
| *El origen de coordenadas está situado en el polo del domo.
**Para aumentar la rigidez debida al conducto conectado al domo. | ||
| 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 |
Se utiliza una malla no estructurada de elementos triangulares de 3 nodos. El tamaño medio asignado a la malla del domo es de 15 mm. Por otro lado, para el cilindro se ha seleccionado un tamaño de 30 mm.
En la Figura 30 se muestran las secciones en las que ha quedado dividido el domo para informar la variación de bobinado de las fibras.
6.1.2 Resultados del modelo hid15r
En este apartado se muestran los resultados obtenidos en la simulación del depósito mallado con elementos de lámina.
La Figura 31 muestra el daño medio de matriz y fibras para la presión de servicio:
Tal y como se puede observar en la Figura 32, el daño se inicia para cuando la presión alcanza los 31,86 MPa (60,35% de la presión máxima). En este caso el daño se inicia en la parte del agujero.
La tensión de Von Mises para la presión de servicio se muestra en la Figura 33:
En la Figura 34 se muestra el daño medio de fibra y matriz para la presión de explosión:
El daño en la fibra para la presión de explosión se muestra en la Figura 35:
6.2 Modelo hid15q
La geometría utilizada es una modificación, ampliando el espesor del compuesto interior del depósito en 32 [mm], en domo y cilindro.
6.2.1 Malla – Consideraciones sobre la discretización
En la figura 36 se muestran las mallas usadas para las simulaciones:
| |
 |
Figura 36. Malla del domo con los anillos para informar el espesor variable.
En este apartado se muestran los resultados obtenidos por simulación para una propuesta de distribución de espesor y dirección de fibras variables en el domo calculada a partir de los criterios propuestos por ACCIONA en la referencia [2]. 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 1. 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 en domo y cilindro respectivamente, por lo que se asume provisionalmente un error mayor que en los casos de espesor constante en el domo.
| 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 | 152 |
| D6 | 82-123 | 130 |
| D5 | 123-164 | 101 |
| D4 | 164-205 | 86 |
| D3 | 205-246 | 76 |
| D2 | 246-287 | 69.5 |
| D1 | 287-314.2 | 65.5 |
| Tapa | 314.2 | 352** |
| Cilindro | 314.2 | 127 |
| *El origen de coordenadas está situado en el polo del domo.
**Para aumentar la rigidez debida al conducto conectado al domo. | ||
6.2.2 Resultados del modelo hid15q
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.
El daño aparece antes en la matriz que en la fibra. Éste se inicia, para una carga de 47.87 MPa (muy próxima a la de explosión), en la zona del domo y después aparece también en la zona de la unión. Se alcanza la presión de explosión y el daño no evoluciona al mantener la carga. El daño para la presión de explosión se muestra en la figura siguiente:
La tensión de Von Mises en la capa 1, se muestra en la siguiente figura. Efectivamente la tensión es mayor en la zona donde se inicia el daño y en la unión. La zona menos solicitada es el cilindro.
La capa exterior está menos solicitada en la zona que inicia el daño y más en la unión. Se muestra en la siguiente figura:
7 DISCUSIÓN DE LOS RESULTADOS
La simulación de láminas “hid15r” no detecta daño en la presión de servicio 22MPa ver Fig.31, no obstante la simulación con elementos hexaédricos tiene un leve daño en la transición domo y cilindro (más localizado en el agujero), ver Fig 10 por eso se adopta a reforzar domo y cilindro con elementos hexaédricos.
La simulación de láminas que no pondera los espesores del modelo “hid15r” da un contour del daño en matriz parecido (0.19max..) al hexaédrico para 52.8Mpa presión de explosión donde es máximo el daño en la transición cilindro y domo ver figuras Fig27; Fig.34.
En la simulación de láminas con espesores ponderados “hid15q” el daño se concentra en la corona intermedia del domo con tensiones máximas entre el domo y cilindro. Ver Fig37 y Fig.39
8 REFERENCIAS
[1] Diseño preliminar tanque para almacenaje de hidrógeno, Autor: no mencionado, Documento: “Tanque_t3.pdf”, Fecha última modificación: 23 de noviembre 2009.
[2] Informe proceso de simulación, autor: No mencionaldo, Documento: “informe simulación031008-2.pdf”, Fecha última modificación: 22 de noviembre de 2008.
[3] Estudio de un depósito de hidrógeno por el Método de los Elementos Finitos y validación analítica, Autor: no mencionado, 08 Abril, 2010.
[4] Jae-Sung Park, Chan-Sun Hong, Chun-Gon Kim, Cheol-Ung Kim, Análisis of filament wound composite structures considering the change of winding angles through the thickness direction, Composite Structures 55, pp 63-71, Ed. Elsevier, 2002.
[5] Cheol-Ung Kim, Chang-Sun Hong, Chun-Gon Kim, Jung-Yub Kim, Optimal design of filament wound type 3 tanks under internal pressure using a modified genetic algorithm, Composite Structures 71, pp 16-25, Ed. Elsevier, 2005.
[6] D.T. Jones, I.A. Jones and V Middleton, Improving composite lay-up for non-spherical filament-wound pressure vessels, Composites, Part A 27 A, pp 311-317, Ed. Elsevier, 1996.
Document information
Published on 01/01/2010
Licence: CC BY-NC-SA license
Share this document
Keywords
claim authorship
Are you one of the authors of this document?