Study of a hydrogen tank by the Finite Element Method and analytical validation

Introducción

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.

Geometría de los modelos estudiados

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:

• Caso A: Semiesfera con condiciones de simetría impuestas sobre su sección ecuatorial (modelo se03.gid, cálculos se03-*.dat).

• Caso B: Geometría del depósito completo (domo semiesférico + cilindro). Cilindro y domo se consideran solidariamente unidos. Se contempla simetría transversal en la sección situada en la mitad de su longitud (modelo hid-15b.gid, cálculos hid-15*.dat).

Materiales y otras condiciones

Otras condiciones de los modelos se describen a continuación:

• No se ha informado la capa interior de aleación de aluminio.

• Orientación de fibras del domo a +10º y -10º constantes para todos los modelos de espesor constante. En el caso del domo de espesor variable la dirección de fibras se ha definido siguiendo los criterios definidos por ACCIONA en el documento tanque_t3.pdf [4].

• Orientación de fibras para las 10 capas más superficiales del cilindro de ±55º y para las 10 capas interiores de ±10º.

• Diámetro de la tapa central del domo de 110mm.

• Apilado en la tapa de 8 capas de fibras orientadas para lograr comportamiento isótropo (0, ±45, 90).

• Tapa rígida con un espesor 10 veces superior al del domo.

• Espesores constantes en domo y cilindro, salvo para el caso de espesor variable, y con valores indicados en cada caso de estudio mostrados en el apartado “Descripción de las simulaciones”.

En las figuras 1 y 2 se muestran las mallas usadas para la simulación:

Figura 1. Malla del domo para el caso A (modelo se03.gid) con designación de los anillos destinados a la información del espesor variable.

Figura 2. Malla para las simulaciones del caso B (modelo hid-15b.gid)

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:

Tabla 1. Propiedades mecánicas de los materiales componentes del laminado

	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

Localización	Matriz	Fibra	% de fibra	Ángulo de fibras	Espesor del layup [mm]
Domo	Matriz polimérica	Fibras de carbono	60	±10º	10/32
Capas exteriores en cilindro	Matriz polimérica	Fibras de carbono	60	±10º	10/32
Capas interiores en cilindro	Matriz polimérica	Fibras de carbono	60	±55º	20/63

Descripción de las simulaciones

Para 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

Caso	Modelo geométrico	Simulación	Espesor cilindro (mm)	Espesor domo (mm)	Espesor tapa central (mm)
A*	se03.gid	se03-f2	-	10	100
A*	se03.gid	se03-k	-	32	320
B**	hid-15b.gid	hid-15k2	30	10	100
B**	hid-15b.gid	hid-15m	95	32	320
B**	hid-15b.gid	hid-15q	95	Variable entre 35 y 120	320
*Semiesfera con simetría **Depósito completo con simetría transversal en cilindro y unión solidaria entre domo y cilindro

Condiciones de carga

La 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

Consideraciones generales

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

(1)

(2)

(3)

(4)

donde: R: es el radio interior del domo o del cilindro
P_int: 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 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

Modelo hid-15k2 (Caso B, espesor domo 10mm y cilindro 30mm)

Este modelo tiene tanto la parte del domo como la parte del cilindro y presenta la distribución del daño en su fase preliminar que se muestra en las figuras 29 y 30. También en esta ocasión el daño en las fibras es nulo. El instante de aparición del daño corresponde a 6,16 MPa y se localiza inicialmente en la parte de los anillos centrales del domo. Por lo tanto la distribución de daño inicial es muy diferente a la que se producía en los modelos sin parte cilíndrica.

Figura 29. Daño en matriz para 6,16 MPa, 28% de la presión nominal

Figura 30. Daño en fibra para 6,16 MPa, 28% de la presión nominal

A continuación se muestran las deformaciones principales a lo largo de toda la geometría (figuras 31 y 32). En este caso también se observa que la rigidez mayor de la tapa produce menores deformaciones en esta zona. En el cilindro, las deformaciones son bastante menores a las del domo. Todos los valores son positivos.

Figura 31. Deformación E1 para 6,16 MPa, 28% de la presión nominal

En el caso de las E2, se producen obtienen valores positivos en toda la geometría excepto en la parte del cilindro que se une con el domo, donde los valores son negativos.

Figura 32. Deformación E2 para 6,16 MPa, 28% de la presión nominal

A continuación se muestra la deformada para la presión nominal, de 22 MPa.

Figura 33. Deformada para la presión nominal 22 MPa

Las siguientes figuras muestran la evolución del daño respecto de la presión (figura 34) y del tiempo (figura 35) en los tres nodos del anillo central de la figura 4.

Figura 34. Evolución del daño en función de la presión de los nodos indicados en figura 4

Figura 35. Evolución del daño de la matriz en función del tiempo en los nodos seleccionados de la zona marcada en la figura 4

En las figuras 36 a 43 se describen los esfuerzos de sección Nxx y Nyy para domo y cilindro. En todos los casos, los esfuerzos se refieren a presión y a tiempo de simulación. Tal y como se ha comentado anteriormente, en el caso analítico de referencia existe un solo Nii que se compara con los Nyy y a Nxx simulados. Los valores, al igual que en todas las simulaciones, se han obtenido en 3 puntos que siguen una línea meridional en el domo y en el cilindro.

Al igual que todas las simulaciones, los valores de Nyy y Nxx simulados se asemejan de forma importante en el caso de la parte del cilindro. En el caso del domo, el Nxx siempre queda por debajo del analítico y el Nyy siempre por encima debido a que el 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. Las dos curvas de referencia en e cilindro corresponden a los valores obtenidos en la parte exterior y en la parte interior ya que el cilindro tiene un espesor importante y los valores difieren entre ellos.

Figura 36. Gráfica de Nxx vs. presión y valor de referencia para el domo según expresión (1)

Figura 37. Gráfica de Nxx vs. tiempo y valor de referencia para el domo según expresión (1)

Figura 38. Gráfica de Nxx vs. presión y valor de referencia para el cilindro según expresión (3)

Figura 39. Gráfica de Nxx vs. tiempo y valor de referencia para el cilindro según expresión (3)

Figura 40. Gráfica de Nyy vs. presión y valor de referencia para el domo según expresión (2)

Figura 41. Gráfica de Nyy vs. tiempo y valor de referencia para el domo según expresión (2)

Figura 42. Gráfica de Nyy vs. presión y valor de referencia para el cilindro según expresión (4)

Figura 43. Gráfica de Nyy vs. tiempo y valor de referencia para el cilindro según expresión (4)

Los desplazamientos que se producen en la tapa fictícia a lo largo del eje z son también positivas y de 45 mm de valor (los valores positivos indican que la tapa va hacia adentro).

Figura 44. Evolución del desplazamiento axial de la tapa con la presión

Modelo hid-15m (Caso B, espesor domo 32 mm y cilindro 95 mm)

Este modelo tiene la geometría del depósito completo (domo semiesférico + cilindro) con una simetría transversal en la sección situada en la mitad de su longitud. Cilindro y domo se consideran solidariamente unidos. La tapa tiene un espesor diez veces mayor que el resto del domo.

Bajo estas condiciones aparece daño a los 18.48 MPa de presión, que corresponden al 84% de la presión nominal. El daño se localiza inicialmente en zonas intermedias del domo. No aparece daño en el cilindro inicialmente y tampoco en las fibras.

Figura 45. Daño en la matriz para 18,48 MPa, 84% de la presión nominal

Figura 46. Daño en la fibra para 18,48 MPa, 84% de la presión nominal

Las siguentes figuras muestran las deformaciones E1 y E2. Las E1 son todas positivas y con valores mayores en las zonas dañadas. Las E2 tienen valores sustancialmente inferiores y de valores tanto positivos como negativos.

Figura 47. Deformación E1 para 18,48 MPa, 84% de la presión nominal

Figura 48. Deformación E2 para 18,48 MPa, 84% de la presión nominal

Las siguientes figuras muestran la evolución del daño respecto de la presión (figura 49) y del tiempo (figura 50) en los tres nodos del anillo central de la figura 4. Se puede observar que aparece daño cerca de la presión nominal (alrededor de los 19 MPa) de una forma brusca y luego, a medida que se incrementa la presión, el daño se incrementa con una pendiente menor.

Figura 49. Evolución del daño en función de la presión de los nodos indicados en figura 4

Figura 50. Evolución del daño de la matriz en función del tiempo en nodos seleccionados de la zona marcada en la figura 4

En las figuras 51 a 58 se describen los esfuerzos de sección Nxx y Nyy para domo y cilindro referidas o a presión o a tiempo de simulación. Los valores de Nyy y Nxx simulados se asemejan de forma importante en el caso de la parte del cilindro. Recordar que en el caso del domo, el Nxx simulado queda por debajo del analítico y el Nyy por encima debido a que el 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.

Figura 51. Gráfica de Nxx vs. presión y valor de referencia para el domo según expresión (1)

Figura 52. Gráfica de Nxx vs. tiempo y valor de referencia para el domo según expresión (1)

Figura 53. Gráfica de Nxx vs. presión y valor de referencia para el cilindro según expresión (3)

Figura 54. Gráfica de Nxx vs. tiempo y valor de referencia para el cilindro según expresión (3)

Figura 55. Gráfica de Nyy vs. presión y valor de referencia para el domo según expresión (2)

Figura 56. Gráfica de Nyy vs. tiempo y valor de referencia para el domo según expresión (2)

Figura 57. Gráfica de Nyy vs. presión y valor de referencia para el cilindro según expresión (4)

Figura 58. Gráfica de Nyy vs. tiempo y valor de referencia para el cilindro según expresión (4)

Finalmente, en cuanto al desplazamiento producido en la tapa ficticia, observamos que los valores son significativamente menores que en el anterior modelo (de espesor inferior).

Figura 59. Evolución del desplazamiento axial de la tapa con la presión

Comparación con resultados analíticos

En este apartado se comparan algunos de los resultados obtenidos en las simulaciones con los resultados obtenidos mediante los cálculos analíticos mostrados en [3]. De este modo se pretende validar, en la medida de lo posible, los resultados por simulación mediante la comparación con los resultados analíticos. Las magnitudes que se comparan son:

• La sobrepresión a la que está sometido el cilindro en la zona de unión con el domo (caso B).

• Los alargamientos anulares del domo en su línea ecuatorial para los casos A y B.

Para una mejor comprensión del aspecto analítico en este documento también se incluirán los cálculos de [3] relacionados, en cada caso, con el parámetro a estudiar.

Cálculo analítico

En este apartado se realiza el cálculo del estado de deformaciones para el domo en estado de servicio. El estado de deformación se obtiene a partir de las tensiones calculadas membranalmente (ver 3.1.1 de [3]) y del tensor constitutivo del material compuesto para el domo (ver Anexo 4.3 de [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.

Cálculo por simulación

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:

Figura 60. Diámetro del domo al 32% de la carga nominal (los colores representan cambios de espesor, véase figura 1)

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:

Cálculo analítico

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.

Cálculo por simulación

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:

Cálculo analítico

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

.

Cálculo por simulación de los esfuerzos de sección y del incremento de presión en la unión

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, λ^C_sim y λ^T_sim 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:

y

Figura 61. Nyy en la unión para domo y cilindro

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:

y

Con lo que el incremento de presión es

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:

Cálculo analítico

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.

Cálculo por simulación

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:

Figura 62. Diámetro en la unión par 6,16 MPa, 28% de la carga nominal (los colores representan cambios de espesor, véase figura 1)

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:

Estudio del caso real con fibras y 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 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.

Tabla 4. Características principales del modelo estudiado

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

Tabla 5. Á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 6. Distribución de espesores

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

Figura 63. Distribución de espesores

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.

Figura 64. Daño en matriz para un 56% de la presión de explosión (39,25MPa)

Figura 65. Daño en fibra para un 56% de la presión de explosión (39,25MPa)

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.

Figura 66. Deformada para 22 MPa, presión nominal

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.

Figura 67. Nxx para un 56% (39,25MPa) de la presión de explosión

Figura 68. Nyy para un 56% (39,25MPa) de la presión de explosión

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.

Figura 69. E11 para un 56% (39,25MPa) de la presión de explosión

Figura 70. E22 para un 56% (39,25MPa) de la presión de explosión

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:

Figura 71. Zona de toma de datos para el modelo hid-15q

Figura 72. Evolución del daño en la matriz con la presión

Figura 73. Gráfica de Nxx vs. presión para el domo según (1)

Figura 74. Gráfica de Nxx vs. tiempo para el domo según expresión (1)

Figura 75. Gráfica de Nxx vs. presión y valor de referencia para el cilindro según expresión (3)

Figura 76. Gráfica de Nxx vs. tiempo y valor de referencia para el cilindro según expresión (3)

Figura 77. Gráfica de Nyy vs. presión para el domo según expresión (2)

Figura 78. Gráfica de Nyy vs. tiempo para el domo según expresión (2)

Figura 79. Gráfica de Nyy vs. presión y valor de referencia para el cilindro según expresión (4)

Figura 80. Gráfica de Nyy vs. tiempo y valor de referencia para el cilindro según expresión (4)

Figura 81. Evolución del desplazamiento axial de la tapa con el tiempo

Conclusiones provisionales

Del estudio de la simulación MEF se puede concluir que:

• En los casos A (semiesfera con simetría) el daño se propaga desde la simetría (línea ecuatorial) hacia la tapa, mientras que en los casos B (domo + cilindro con simetría transversal) la zona de inicio del daño depende de la distribución del espesor:

• • Si el espesor es constante el daño comienza en la zona media y a partir de esta se propaga hacia el resto del domo.

• • Si el espesor es variable, según la distribución de la tabla 4 y figura 35, el daño se inicia en la zona de unión con el cilindro (con espesor mínimo) y se propaga hacia la zona media.

• Como era de esperar el caso de mayores espesores (32 mm en domo y 95 en el cilindro) soporta mayor presión pero sin alcanzar la presión nominal (22 MPa) por la aparición previa de daño en la matriz del domo (para el 84% de P_nom). Analíticamente se había establecido que este fallo se produciría para un 37% de la P_nom (ver [3] cuadro de la página 20).

• En ningún caso se detecta daño en las fibras.

• Se ha observado (figura 16) que una vez que la matriz ya está dañada, este daño puede aumentar aún a presión constante.

• Las fuerzas específicas de referencia en el domo y el cilindro se han calculado asumiendo la hipótesis de paredes delgadas y material elástico lineal e isótropo. Por este motivo deben ser tomadas como una referencia aproximada.

• Los valores de los esfuerzos de sección (Nxx y Nyy) obtenidos en el compuesto se atenúan o mantienen constante después de iniciado el daño en la matriz por lo que no pueden alcanzar los valores de referencia que suponen comportamiento lineal y elástico.

• Se observa que los esfuerzos de sección en la dirección longitudinal (siguiendo los meridianos del domo, Nyy) se asemejan mucho a los valores de referencia mientras que las tensiones en la dirección circunferencial (siguiendo los paralelos, Nxx) son sustancialmente menores. Este efecto puede ser debido a la menor rigidez del domo en la dirección circunferencial (puesto que las fibras están dispuestas a 10º resp. la dirección longitudinal).

• Se piensa que se obtendrían mejores resultados utilizando una resina con menor módulo elástico, de tal forma que al deformarse más transfiera la carga a la fibra dejando así que ésta desarrolle su capacidad resistente y a la vez se protegería la matriz al no haber llegado al estado de tensión máxima por ser más flexible.

• En el apartado Estudio del caso real con fibras y espesor variable se aprecia que con la distribución de espesor y fibra variables indicados en las tablas 3 y 4 se obtienen resultados próximos a los requeridos en las especificaciones. Cabe destacar de nuevo, que el error en la definición del plano medio puede afectar a los resultados, y por tanto se debería estudiar este caso con la geometría exacta y una técnica desarrollada add hoc para la aplicar una presión equivalente.

Referencias

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