Cálculo de un depósito de hidrógeno

Latest revision as of 15:21, 28 February 2018

Simulación realizada con FEMCOM 1.0. Modelos con elementos hexaedros

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.

Figura 1. Sistema de coordenadas y algunas variables del problema

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.

Figura 2. Dimensiones del tanque de hidrógeno en milímetros. Extraído de [1]

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]:

(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.

Gráfico 1. Variación del espesor en el domo respecto a r

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).

Figura 3: Geometría del depósito simulada

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:

Figura 4. Propiedades asignadas a las fibras de carbono

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:

Figura 5: Propiedades asignadas a la matriz de epoxi

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.

Figura 6. Propiedades asignadas al liner de aluminio

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%.

Figura 7. Porcentaje de fibra en el compuesto

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.

Gráfico 2: Variación del ángulo de bobinado en el domo respecto a r para las fibras con ángulo de entrada +10º y -10º

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.

Gráfico 3. Presión interna aplicada a lo largo del tiempo

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.

Figura 8. Número de divisiones de las líneas

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.

Figura 9. Secciones para informar la variación del ángulo de bobinado en el domo respecto al ángulo de 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:

Figura 10. N Damage a la presión de servicio (22 MPa)

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.

Figura 11. Daño en la matriz (izquierda) y en la fibra (derecha), sin daño, para 22 MPa

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.

Figura 12. Instante de inicio del daño en la unión

En la Figura 13 se muestra el daño en el compuesto para 25,6 MPa:

Figura 13. Daño promedio de matriz y fibras para presión de 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.

Figura 14. Daño en la matriz (izquierda) y en la fibra (derecha), sin daño, para 25,6 MPa

La tensión a tracción para la presión de servicio se muestra en la Figura 15:

Figura 15. Tensión principal 1 (tracción) a la presión de servicio (22 MPa)

La tensión principal en la dirección 2 se muestra en la Figura 16:

Figura 16. Tensión principal 2 a la presión de servicio (22 MPa)

La tensión principal en la dirección 3 se muestra en la Figura 17:

Figura 17. Tensión principal 3 a la presión de servicio (22 MPa)

En la Figura 18 se muestra el daño medio de fibra y matriz para la presión de explosión:

Figura 18. N Damage para una presión de 52,8 MPa

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º.

Figura 19. Geometría resultante del modelo con refuerzo

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.

Figura 20. Malla de la geometría con refuerzo

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.

Figura 21. N Damage a la presión de servicio (22 MPa)

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.

Figura 22. Daño en el compuesto para presión intermedia (25,6 MPa)

Figura 23. Daño en la matriz (para las fibras es de 0) para presión intermedia (25,6 MPa)

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.

Figura 24. Tensión principal 1 (tracción) a la presión de servicio (22 MPa)

Figura 25. Tensión principal 2 a la presión de servicio (22 MPa)

Figura 26. Tensión principal 3 a la presión de servicio (22 MPa)

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).

Figura 27. N Damage para una presión de 52,8 MPa

Figura 28 Daño en la fibra para una presión de 52,8 MPa

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:

Figura 29. División de superfícies en el domo

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.

Tabla 2. Ángulo de bobinado 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

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.

Figura 30. Secciones para informar variación del ángulo de bobinado en el domo

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:

Figura 31. N Damage a la presión de servicio (22 MPa)

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.

Figura 32. Instante de inicio del daño

La tensión de Von Mises para la presión de servicio se muestra en la Figura 33:

Figura 33. Tensión de Von Mises a la presión de servicio (22 MPa)

En la Figura 34 se muestra el daño medio de fibra y matriz para la presión de explosión:

Figura 34. N Damage para una presión de 52,8 MPa

El daño en la fibra para la presión de explosión se muestra en la Figura 35:

Figura 35. Daño en la fibra para una presión de 52,8 MPa

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.

Tabla 3. Ángulo de apilado variable

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

Tabla 4. Distribución de espesores

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:

Figura 37. Daño en matriz a 52.8 MPa (carga de explosión).

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.

Figura 38. Tensión de Von Mises cara interior a 52.8 MPa (carga de explosión).

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:

Figura 39. Tensión de Von Mises cara exterior a 52.8 MPa (carga de explosión).

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.