Avances en la formulación de membranas

Avances en la formulación de membranas

- 3 -

CIMNE Reporte Interno. Mayo 2004

- 1 -

1 Avances en la Formulación de Membranas

1.2 CIMNE

Centro Internacional de Métodos Numéricos en Ingeniería

Edificio C1, Jordi Girona 1-3

08034 Barcelona, España

Resumen

Avances en el análisis mediante elementos finitos de membranas curvas con wrinkling sujetas a grandes desplazamientos son presentados. Este reporte trata el estudio de membranas por el método de los elementos finitos mediante un análisis con no linealidad geométrica para materiales hiperelásticos. Una formulación sencilla en términos del tensor gradiente de deformación es usada para elementos finitos triangulares lineales. También se hace el análisis de membranas presforzadas, como se encuentran generalmente en la vida cotidiana.

1. Introducción

Las membranas estructurales tienen muchas aplicaciones en nuestros días y están caracterizadas porque solamente poseen fuerzas axiales (de membrana). Este trabajo es una aportación a la teoría del análisis de membranas estructurales propuesto por Taylor [1], en donde ahora es posible efectuar el análisis de membranas con propiedades ortótropas. El objetivo principal en el análisis ortótropo de membranas es la correcta asignación de las direcciones del refuerzo principal a cada elemento finito de la estructura. Para comprender bien el problema de usar la teoría de Taylor para analizar membranas con propiedades ortótropas se hace el siguiente desarrollo.

Consideremos un sistema global cartesiano donde las coordenadas de la configuración de referencia se representan por ${\textstyle {X}}$ y las coordenadas de la configuración deformada por ${\textstyle {x}}$. Las coordenadas para cada elemento finito triangular serán entonces ${\textstyle X^{i}}$ y ${\textstyle x^{i}}$ para ${\textstyle i=1,2,3}$. Además un sistema local de coordenadas para cada elemento es usado, donde la letra ${\textstyle {Y}}$ representa la configuración de referencia local cuyas coordenadas son ${\textstyle {Y}^{j}}$, y la letra ${\textstyle {y}}$ es la configuración actual local con coordenadas ${\textstyle {y}^{j}}$ para ${\textstyle j=1,2,3}$. El origen del sistema local de coordenadas se encuentra en la coordenada ${\textstyle X^{1}}$ y ${\textstyle x^{1}}$ para la configuración de referencia y la actual respectivamente como se puede ver en la Fig. 1. Los vectores unitarios para el sistema coordenado local se construyen alineando el primer vector unitario con el lado de conectividades 1-2 del elemento. Entonces el primer vector unitario se expresa por

$v_{1}={\frac {x_{2}-x_{1}}{\left\|x_{2}-x_{1}\right\|}}$

(1)

donde

$\left\|x^{2}-x^{1}\right\|=\left[\left(x^{2}-x^{1}\right)^{T}\left(x^{2}-x^{1}\right)\right]^{\frac {1}{2}}$

(2)

Un vector normal al triángulo es construido con

$v_{3}=\left(x^{2}-x^{1}\right)\times \left(x^{3}-x^{1}\right)$

(3)

$X_{1}X_{2}X_{3}v_{3}nx^{1}v_{2}v_{1}x^{2},y^{2}x^{3},y^{3}{y}_{1}{y}_{2}X_{1}X_{2}X_{3}v_{3}nx^{1}v_{2}v_{1}x^{2},y^{2}x^{3},y^{3}{y}_{1}{y}_{2}$

Fig.1.- Coordenadas globales y locales para un elemento triangular

Después, el vector normal unitario se define por

$n={\frac {v_{3}}{\left\|v_{3}\right\|}}$

(4)

Finalmente el segundo vector unitario es

$v_{2}=n\times v_{1}$

(5)

Todos los desarrollos anteriores se han efectuado con elementos en configuración actual, y para obtener la configuración de referencia es suficiente con cambiar las letras minúsculas por mayúsculas.

Las coordenadas globales de cada elemento pasan a cada sistema de referencia local (dentro del plano de cada elemento finito triangular) usando la siguiente ecuación en función de los vectores unitarios

$y_{i}^{j}=\left(x^{j}-x^{1}\right)^{T}\cdot v_{i}$

(6)

donde ${\textstyle i=1,2}$ y ${\textstyle j=1,2,3}$. La interpolación isoparamétrica para el elemento es

${y}={\xi }_{j}y^{j}$

(7)

donde ${\textstyle {\xi }_{j}}$ son las coordenadas naturales (de área) que satisfacen

$\xi _{1}+\xi _{2}+\xi _{3}=1$

(8)

De la Ec. 6 se puede ver que para ${\textstyle j=1}$, ${\textstyle y_{i}^{1}=0}$. Entonces la Ec. 7 se puede reducir a

${y}=\xi _{2}y^{2}+\xi _{3}y^{3}$

(9)

2. Tensor Gradiente de Deformación

El tensor gradiente de deformación para cada elemento finito triangular que describe el movimiento desde la configuración de referencia a la configuración actual dentro del plano del elemento es

$F={\frac {\partial {y}}{\partial {Y}}}$

(10)

Usando la formulación isoparamétrica dada por el Ec. GrindEQ__9_, pero en configuración de referencia, el tensor gradiente de deformación se puede expresar por

$F{\frac {\partial {Y}}{\partial {\xi }}}={\frac {\partial {y}}{\partial {Y}}}{\frac {\partial {Y}}{\partial {\xi }}}={\frac {\partial {y}}{\partial {\xi }}}$

(11)

y el Jacobiano de transformación está dado por

$J={\frac {\partial {Y}}{\partial {\xi }}}{\;\;}{and\;\;}j={\frac {\partial {y}}{\partial {\xi }}}$

(12)

donde ${\textstyle J}$ es el Jacobiano en configuración de referencia y ${\textstyle j}$ es el Jacobiano para la configuración deformada, y entonces

$F=j{\;}J^{-1}$

(13)

Para poder desarrollar el cálculo del Jacobiano, de la Ec. GrindEQ__6_ se deduce

$y^{2}=\left[{\begin{array}{c}{\left(x^{2}-x^{1}\right)^{T}\cdot v_{1}}\\{\left(x^{2}-x^{1}\right)^{T}\cdot v_{2}}\end{array}}\right]{\;\;\;}{and\;\;\;}y^{3}=\left[{\begin{array}{c}{\left(x^{3}-x^{1}\right)^{T}\cdot v_{1}}\\{\left(x^{3}-x^{1}\right)^{T}\cdot v_{2}}\end{array}}\right]$

(14)

y el Jacobiano en configuración deformada se calcula a partir de la Ec. GrindEQ__9_ y Ec. GrindEQ__12_

$j={\frac {\partial y}{\partial {\xi }}}=\left[{\begin{array}{cc}{\frac {\partial y_{2}^{2}}{\partial \xi _{2}}}&{\frac {\partial y_{2}^{3}}{\partial \xi _{3}}}\\{}&{}\\{\frac {\partial y_{3}^{2}}{\partial \xi _{2}}}&{\frac {\partial y_{3}^{3}}{\partial \xi _{3}}}\end{array}}\right]=\left[{\begin{array}{cc}{\left(x^{2}-x^{1}\right)^{T}\cdot v_{1}}&{\left(x^{3}-x^{1}\right)^{T}\cdot v_{1}}\\{\left(x^{2}-x^{1}\right)^{T}\cdot v_{2}}&{\left(x^{3}-x^{1}\right)^{T}\cdot v_{2}}\end{array}}\right]$

(15)

Similarmente para la configuración de referencia obtenemos

$J=\left[{\begin{array}{cc}{\left(X^{2}-X^{1}\right)^{T}\cdot V_{1}}&{\left(X^{3}-X^{1}\right)^{T}\cdot V_{1}}\\{\left(X^{2}-X^{1}\right)^{T}\cdot V_{2}}&{\left(X^{3}-X^{1}\right)^{T}\cdot V_{2}}\end{array}}\right]$

(16)

Notando que ${\textstyle x^{2}-x^{1}}$ es ortogonal a ${\textstyle v_{2}}$ en la configuración actual y de manera similar para la configuración de referencia, las Ecs. GrindEQ__15_ y GrindEQ__16_ se pueden simplificar a

$j=\left[{\begin{array}{ccc}{\left\|x^{2}-x^{1}\right\|}&{}&{\frac {\left(x^{3}-x^{1}\right)^{T}\left(x^{2}-x^{1}\right)}{\left\|x^{2}-x^{1}\right\|}}\\{}&{}&{}\\{0}&{}&{\frac {\left\|v^{3}\right\|}{\left\|x^{2}-x^{1}\right\|}}\end{array}}\right]$

(17)

y

$J=\left[{\begin{array}{ccc}{\left\|X^{2}-X^{1}\right\|}&{}&{\frac {\left(X^{3}-X^{1}\right)^{T}\left(X^{2}-X^{1}\right)}{\left\|X^{2}-X^{1}\right\|}}\\{}&{}&{}\\{0}&{}&{\frac {\left\|V^{3}\right\|}{\left\|X^{2}-X^{1}\right\|}}\end{array}}\right]$

(18)

Entonces el tensor derecho de deformación de Cauchy-Green se puede expresar por

$C=F^{T}F$

(19)

y substituyendo la Ec. 13 en la Ec. 19 obtenemos

$C=J^{-T}{\;}j^{T}{\;}j{\;}J^{-1}=G^{T}g{\;}G$

(20)

donde ${\textstyle G}$ es usado para denotar la inversa de ${\textstyle J}$. Expandiendo la Ec. 20 encontramos

$C={\frac {1}{J_{11}^{2}J_{22}^{2}}}\left[{\begin{array}{cc}{J_{22}}&{0}\\{-J_{12}}&{J_{11}}\end{array}}\right]\left[{\begin{array}{cc}{g_{11}^{}}&{g_{12}}\\{g_{12}}&{g_{22}^{}}\end{array}}\right]\left[{\begin{array}{cc}{J_{22}}&{-J_{12}}\\{0}&{J_{11}}\end{array}}\right]$

(21)

donde

$g=\left[{\begin{array}{ccc}{\begin{array}{c}{\left\|x^{2}-x^{1}\right\|^{2}}\\{}\end{array}}&{}&{\begin{array}{c}{\left(x^{3}-x^{1}\right)^{T}\left(x^{2}-x^{1}\right)}\\{}\end{array}}\\{}&{}&{}\\{\left(x^{3}-x^{1}\right)^{T}\left(x^{2}-x^{1}\right)}&{}&{\frac {\left[\left(x^{3}-x^{1}\right)^{T}\left(x^{2}-x^{1}\right)\right]^{2}+\left\|v^{3}\right\|^{2}}{\left\|x^{2}-x^{1}\right\|^{2}}}\end{array}}\right]$

(22)

Simplificando la Ec. GrindEQ__22_ se obtiene

${\begin{array}{l}{g_{11}=j_{11}^{2}=\left(x^{2}-x^{1}\right)^{T}\cdot \left(x^{2}-x^{1}\right)}\\{g_{12}=j_{12}\cdot j_{11}=\left(x^{3}-x^{1}\right)^{T}\cdot \left(x^{2}-x^{1}\right)}\\{g_{22}=j_{12}^{2}+j_{22}^{2}=\left(x^{3}-x^{1}\right)^{T}\cdot \left(x^{3}-x^{1}\right)}\end{array}}$

(23)

3. Ecuación Constitutiva Isótropa y Ortótropa

La relación entre las tensiones y las deformaciones usadas en este trabajo viene dada por

$S=DE$

(24)

donde ${\textstyle S}$ es el segundo tensor de tensiones de Piola-Kirchhoff, ${\textstyle E}$ es el tensor de deformación de Green-Lagrange que se calcula mediante

$E={\frac {1}{2}}\left(C-I\right)$

(25)

Ya que las membranas poseen solamente fuerzas axiales, la ecuación constitutiva ${\textstyle D}$ es expresada para un material isótropo por la ecuación de tensión plana

$D={\frac {E}{1-v^{2}}}\left[{\begin{array}{ccc}{1}&{v}&{0}\\{v}&{1}&{0}\\{0}&{0}&{\left(1-v\right){\mathord {\left/{\left(1-v\right)2}\right.}}2}\end{array}}\right]$

(26)

en la cual ${\textstyle E}$ es el modulo elástico y ${\textstyle v}$ es el coeficiente de Poisson (Zienkiewicz y Taylor [2]). Sin embargo para un material ortótropo la ecuación de tensión plana es

$D={\frac {1}{1-v_{yx}v_{xy}}}\left[{\begin{array}{ccc}{E_{x}}&{v_{yx}E_{x}}&{0}\\{v_{xy}E_{y}}&{E_{y}}&{0}\\{0}&{0}&{\left(1-v_{yx}v_{xy}\right)G_{xy}}\end{array}}\right]$

(27)

Nuestro análisis de membranas está enfocado al análisis para velas, donde los materiales ortótropos son usados para simular un comportamiento más realista que con materiales isótropos. El material para una vela se puede ver en la Fig. 2.

Fig. 2.- Material ortótropo para una vela

Suponemos que una parte de una vela es mallada con elementos finitos triangulares para ser analizada usando nuestra teoría de membranas ortótropas. Aunque trabajemos con mallas estructuradas, el sistema local de coordenadas para cada elemento puede tener diferentes orientaciones como se puede ver en la Fig. 3. Hay que tener en cuenta que los vectores unitarios ${\textstyle v^{1}}$ y ${\textstyle v^{2}}$ se encuentran en el plano de cada elemento triangular. Si suponemos en la línea punteada roja indica la dirección de las fibras de refuerzo principal del material, nuestro problema consiste en encontrar el ángulo ${\textstyle \theta }$ que necesitamos rotar nuestro sistema local para poder aplicar nuestra teoría.

Fig. 3.-Dirección de fibras de refuerzo principal y vectores unitarios locales

Supongamos que ${\textstyle v^{d}}$ es el vector que define la dirección de la fibra principal de refuerzo para un elemento y que es coplanar con el mismo elemento como en la Fig. 4. Más adelante veremos como obtener ${\textstyle v^{d}}$ para cada elemento finito triangular.

$\theta v^{d}\theta v^{d}$

Fig. 4.-Vector de dirección de la fibra principal

El análisis ortótropo puede ser efectuado calculando primero la ecuación constitutiva local ${\textstyle D^{L}}$ y después rotándola para obtener la ecuación constitutiva global ${\textstyle D^{G}}$ usando

$D^{G}=T^{T}D^{L}T$

(28)

donde ${\textstyle T}$ es la matriz de rotación definida por Oñate [3] como

$T=\left[{\begin{array}{ccc}{\cos ^{2}\theta }&{\sin ^{2}\theta }&{-2\sin \theta \cos \theta }\\{\sin ^{2}\theta }&{\cos ^{2}\theta }&{2\sin \theta \cos \theta }\\{\sin \theta \cos \theta }&{-\sin \theta \cos \theta }&{\cos ^{2}\theta -\sin ^{2}\theta }\end{array}}\right]$

(29)

en la cual

$\sin \theta =-v^{d}\cdot v^{2}{\;\;}{and\;\;}\cos \theta =v^{d}\cdot v^{1}$

(30)

4. Asignación de la Dirección del Refuerzo Principal

Cuando una membrana pertenece a un plano, la asignación del vector de dirección de refuerzo principal ${\textstyle v^{d}}$ no tiene ningún problema. Simplemente asignamos ese vector a cada elemento finito. Pero cuando la membrana no se encuentra dentro de un plano, como por ejemplo las velas, la asignación de la dirección del refuerzo no resulta tan sencilla. Para abordar este problema hemos propuesto la metodología siguiente.

Primero hay que escoger cualquier elemento finito de la malla y asignarle la dirección ${\textstyle v^{d}}$ que será la dirección del refuerzo principal (elemento de referencia), como se muestra en la Fig. 5. El vector ${\textstyle v^{d}}$ en ese elemento se copiará a todos los elementos con que tengan un lado en común (elementos de destino), con su correspondiente dirección dentro del plano del elemento de destino (elementos en gris de la Fig. 5).

Para copiar el vector ${\textstyle v^{d}}$ del elemento de referencia a uno de los elementos de destino se necesita identificar cual es el lado en común entre ambos elementos. De la Fig. 6(a) podemos ver que el lado ${\textstyle i-j}$ es el lado en común. El punto medio de ese lado se definirá como el punto ${\textstyle o}$, y a continuación construiremos el siguiente sistema de coordenadas local coplanar con el elemento de referencia. Como la normal de cada elemento es conocida, pasamos la normal del elemento de referencia ${\textstyle n_{1}}$ al punto ${\textstyle o}$. Después el eje ${\textstyle x_{1}}$ se construye desde el punto ${\textstyle o}$ al nodo ${\textstyle j}$. Finalmente el eje ${\textstyle y_{1}}$ se construye mediante un producto cruz entre ${\textstyle n_{1}{\;y\;}x_{1}}$.

Fig. 5.-Elemento con su vector de refuerzo principal

${\textstyle _{n_{1}}}$${\textstyle _{n_{2}}}$${\textstyle _{i}}$${\textstyle _{j}}$${\textstyle _{n_{2}}}$${\textstyle _{x_{2}}}$${\textstyle _{y_{2}}}$${\textstyle _{_{a}}}$${\textstyle _{_{b}}}$${\textstyle _{_{a'}}}$${\textstyle _{_{b'}}}$${\textstyle _{o}}$${\textstyle _{o}}$${\textstyle _{x}}$${\textstyle _{_{1}}}$${\textstyle _{n}}$${\textstyle _{_{1}}}$${\textstyle _{y}}$${\textstyle _{_{1}}}$${\textstyle _{n_{1}}}$${\textstyle _{n_{2}}}$${\textstyle _{i}}$${\textstyle _{j}}$${\textstyle _{n_{2}}}$${\textstyle _{x_{2}}}$${\textstyle _{y_{2}}}$${\textstyle _{_{a}}}$${\textstyle _{_{b}}}$${\textstyle _{_{a'}}}$${\textstyle _{_{b'}}}$${\textstyle _{o}}$${\textstyle _{o}}$${\textstyle _{x}}$${\textstyle _{_{1}}}$${\textstyle _{n}}$${\textstyle _{_{1}}}$${\textstyle _{y}}$${\textstyle _{_{1}}}$El sistema local del elemento de destino dentro de su plano es construido pasando su normal ${\textstyle n_{2}}$ al punto ${\textstyle o}$. El eje ${\textstyle x_{2}}$ se define como en vector que va desde el punto ${\textstyle o}$ al nodo ${\textstyle i}$. Finalmente el eje ${\textstyle y_{2}}$ se encuentra con el producto cruz de ${\textstyle n_{2}{\;y\;}x_{2}}$.

Fig. 6.-Copia de la dirección de refuerzo

Una vez que ambos sistemas están construidos, el siguiente paso es escoger un punto arbitrario ${\textstyle a}$ en el eje ${\textstyle x_{1}}$ como se ilustra en la Fig. 6(b). Haciendo pasar el vector ${\textstyle v^{d}}$ de refuerzo principal del elemento de referencia por el punto ${\textstyle a}$, encontramos sobre el eje ${\textstyle y_{1}}$ el punto ${\textstyle b}$. Así el vector que va desde el punto ${\textstyle b}$ al punto ${\textstyle a}$ tendrá la misma dirección que el vector ${\textstyle v^{d}}$.

Finalmente para copiar ese vector al elemento de destino, la distancia del punto ${\textstyle o}$ al punto ${\textstyle a}$ en el eje ${\textstyle x_{1}}$ debe ser la misma que del punto ${\textstyle o}$ al punto ${\textstyle a'}$ en el eje ${\textstyle x_{2}}$, y la distancia del punto ${\textstyle b}$ sobre ${\textstyle y_{1}}$ debe ser la misma que desde ${\textstyle o}$ a ${\textstyle b'}$ sobre ${\textstyle y_{2}}$. Con los punto ${\textstyle a',{\;}b'}$ definidos dentro del plano del elemento de destino, construimos un vector ${\textstyle {\overline {a'b'}}}$ que define el vector de refuerzo principal en el elemento de destino. El proceso se repite hasta que todos los elementos de la malla tienen un vector que define su dirección de refuerzo principal.

5. Forma Débil de las Ecuaciones de Movimiento

La forma débil para una membranas se escribe con el trabajo virtual como

$\int _{\Omega }\delta x^{i}\rho _{o}h{\ddot {x}}^{i}d\Omega +\int _{\Omega }\delta x^{i}c_{o}h{\dot {x}}^{i}d\Omega +\int _{\Omega }\delta E_{IJ}S_{IJ}hd\Omega =\int _{\omega }\delta x^{i}b_{i}d\omega +\int _{\gamma }\delta x^{i}t_{i}d\gamma$

(31)

en donde ${\textstyle \rho _{o}}$ es la densidad en la configuración de referencia, ${\textstyle c_{o}}$ es un amortiguamiento lineal en la configuración de referencia, ${\textstyle h}$ es el espesor de la membrana, ${\textstyle S_{IJ}}$ son las componentes del segundo tensor de tensiones de Piola-Kirchhoff, ${\textstyle b_{i}}$ son las componentes de las fuerzas másicas en coordenadas globales y ${\textstyle t_{i}}$ son las componentes de las fuerzas de superficie de las membranas por unidad de longitud. Los términos en la parte izquierda de la Ec. GrindEQ__31_ están escritos en la configuración de referencia mientras que los de la parte derecha lo están en configuración actual. La forma residual de la Ec. 31 para cada elemento se escribe por

$\left[{\begin{array}{c}{R^{1}}\\{R^{2}}\\{R^{3}}\end{array}}\right]=\left[{\begin{array}{c}{f^{1}}\\{f^{2}}\\{f^{3}}\end{array}}\right]-M\left[{\begin{array}{c}{{\ddot {x}}^{1}}\\{{\ddot {x}}^{2}}\\{{\ddot {x}}^{3}}\end{array}}\right]-C\left[{\begin{array}{c}{{\dot {x}}^{1}}\\{{\dot {x}}^{2}}\\{{\dot {x}}^{3}}\end{array}}\right]-hAB^{T}\left[{\begin{array}{c}{S_{11}}\\{S_{22}}\\{S_{12}}\end{array}}\right]$

(32)

donde

$B^{T}=b^{T}Q^{T}$

(33)

en la cual

$b^{T}=\left[{\begin{array}{ccccc}{\left(x^{1}-x^{2}\right)}&{}&{\left(x^{1}-x^{3}\right)}&{}&{\left(x^{1}-x^{2}\right)+\left(x^{1}-x^{3}\right)}\\{\left(x^{2}-x^{1}\right)}&{}&{0}&{}&{\left(x^{3}-x^{1}\right)}\\{0}&{}&{\left(x^{3}-x^{1}\right)}&{}&{\left(x^{2}-x^{1}\right)}\end{array}}\right]$

(34)

y

$Q^{T}=\left[{\begin{array}{ccc}{G_{11}^{2}}&{G_{12}^{2}}&{2G_{11}G_{12}}\\{0}&{G_{22}^{2}}&{0}\\{0}&{G_{12}G_{22}}&{G_{11}G_{22}}\end{array}}\right]$

(35)

Para resolver la Ec. GrindEQ__32_ para todo el sistema necesitamos hacer una integración temporal y otra integración espacial. Para la integración temporal podemos usar el método de Newmark y para la integración espacial de elementos finitos usamos el método completo de Newton-Raphson. Una completa explicación de la solución numérica puede encontrarse en Taylor [1] y ejemplos en Valdés [4].

6. Wrinkling

En el uso de membranas estructurales en la construcción se presta especial atención a los fenómenos de wrinkling (arrugas) y folding (doblado). Desde el punto de vista estético, siempre se ha preferido que los fenómenos mencionados anteriormente no estén presentes en la estructura, pero desde un punto de vista ingenieril, lo que importa es la adecuada medición de las tensiones a las que estará sujeta la estructura para que así se pueda llevar a cabo un adecuado diseño.

Para saber cuando una membrana está sujeta a tensiones en todas sus direcciones o no, se sigue el siguiente criterio basado en tensiones y deformaciones principales de la forma

6.1 Wrinkling basado en la eliminación de tensiones

Esta forma de calcular las tensiones reales sobre una estructura consiste en una vez que se identifica el wrinkling en la estructura mediante la Ec. GrindEQ__36_ se procede a eliminar directamente la segunda tensión principal. Para ello el tensor de deformaciones se descompone en dos partes, una parte elástica que contiene las deformaciones positivas y una parte de wrinkling, que en teoría debe contener solamente las tensiones negativas. Esto se expresa mediante la siguiente ecuación

Para el caso de que exista wrinkling en la estructura se precede de la siguiente manera

donde ${\textstyle {S}_{I}=S_{I}{v}_{I}\otimes {v}_{I}}$ y la deformación buscada se encuentra mediante

Esta forma de calcular el wrinkling da resultados aceptables, pero tiene el inconveniente que presenta una convergencia demasiado mala al no tomar en cuenta para nada las segundas tensiones principales, cosa que sí hace el siguiente apartado.

6.2 Wrinkling basado en la eliminación de tensiones combinado con deformaciones

Esta forma de calcular el wrinkling es

donde ${\textstyle \alpha }$es un porcentaje de ${\textstyle E_{II}}$. El problema consiste nuevamente en hacer que la segunda tensión principal de la Ec. GrindEQ__40_ valga cero. A diferencia del método anterior, este método al tomar en cuenta ambas direcciones principales presenta una convergencia mucho mejor a la anterior, donde al principio del cálculo parece ser igual de lento que el método anterior, pero al paso de los distintos pasos de tiempo, las iteraciones comienzan a ser menores hasta cuando prácticamente llega a efectuar el mismo numero de iteraciones que en el caso en que no se considere el wrinkling (convergencia cuadrática).

La forma en que se encuentra el valor de ${\textstyle \alpha }$ es muy sencilla y se puede hacer por muchos métodos para encontrar ceros de polinomios, como lo son el método de bisección, de la regla falsa, de Newton-Raphson, de la secante, etc. Para mayor claridad sobre ambos métodos se pueden ver las referencias [5] y [6].

7. Pretensado

Las membranas se pueden aplicar en muchas situaciones cotidianas, pero según la aplicación estás tendrán condiciones iniciales diferentes. En el caso de estructuras hinchables, es decir estructuras cerradas sometidas a presión interna que es quien define su forma, basta con la formulación de la membrana más el algoritmo del wrinkling.

Sin embargo, en situaciones donde las membranas son usadas como en las alas de los aviones o como recubrimientos, las condiciones iniciales consisten en imponer unas ciertas tensiones de pretensado a la estructura para darle forma inicial aceptable. Este tipo de tensiones se pueden introducir directamente en la formulación de manera que las fuerzas internas siempre tomen en cuenta estos valores así como la matriz de rigidez geométrica para lograr una convergencia cuadrática además de encontrar la simulación y resultados adecuados al problema. Una forma más sencilla de tratar este problema es convertir las tensiones de pretensado directamente en unas deformaciones iniciales. De esta manera se evita el problema de introducir el efecto de pretensado en las fuerzas internas y en la matriz de rigidez geométrica, el cual se toma en cuenta indirectamente a partir de las deformaciones.

8. Implementación con Elementos Finitos

La formulación de la membrana se implementa mediante elementos finitos tomando en cuenta que se trata de un problema con no linealidad geométrica y se tiene que resolver con métodos apropiados. Así mismo al ser un problema dinámico, se tiene que tener en cuenta la solución temporal de dichas ecuaciones.

El primer paso consiste en calcular correctamente las fuerzas residuales del problema, las cuales vienen dadas por las fuerzas internas y las externas (de la presión de la estructura).

Asimismo se debe tomar en cuenta la matriz de rigidez del problema, que a su vez dependerá de una rigidez material, una rigidez geométrica y una rigidez dada por las fuerzas de presión (que siempre son normales a la superficie).

Para completar el planteamiento mediante elementos finitos del problema se tiene que hacer la parte dinámica, la cual se puede resolver mediante el método de Newmark, para el siguiente caso

Hay que tener en cuenta que el único inconveniente de la Ec. GrindEQ__41_ es encontrar debidamente el valor de la matriz de amortiguamiento ${\textstyle {C}}$. Para lograrlo hay que darle diferentes valores (dependiendo de las dimensiones y unidades del problema a resolver) hasta encontrar el valor que nos proporcione una matriz de amortiguamiento crítico, pero teniendo cuidado de no ir a pasarnos de este valor, ya que de esta manera aparentaremos haber llegado a la solución, siendo en realidad que al haberle dado un amortiguamiento mayor al crítico, estamos sobrerigidizando la estructura, impidiendo que llegue a sus tensiones y deformaciones reales. Otra forma de resolver el problema es mediante la ecuación

En este caso el problema se puede resolver mediante el método de Newmark (en el cual se desperdicia parte del método al no tener en cuenta la aceleración) o de una mejor manera mediante un método de backward Euler por ejemplo. Para ello se encontró que usando una matriz de amortiguamiento igual a la matriz de masa el problema daba buenas soluciones. La única diferencia entre esta forma de resolver el problema y la anterior es que la segunda maneja time steps de un valor mucho menor al de la primera. Además se tiene el inconveniente de que el tiempo manejado en la Ec. GrindEQ__42_ no se puede considerar un tiempo real a la hora de hacer estimaciones de tiempo. Sin embargo el número de time steps e iteraciones es similar en ambos casos.

Una vez que se tiene el problema linealizado y teniendo el cuenta la parte dinámica ya sea por la Ec. GrindEQ__41_ o la Ec. GrindEQ__42_ se procede a resolver el problema mediante le método completo de Newton-Raphson.

9. Postproceso

El postproceso de membranas estructurales curvas ha sido durante mucho tiempo un problema sin resolver debido a que los refuerzos en las membranas (fibras principales) siguen una dirección, la cual generalmente está fuera de un plano. Lo más problemático llega cuando una membrana estructural que está formada por ensamblajes de varias partes del mismo o diferentes materiales, hace cambiar constantemente la dirección de la membrana. De esta forma el postproceso se vuelve un problema casi imposible de resolver.

Lo que generalmente se venía haciendo era postprocesar la información pero solo con valores y direcciones principales. Esta idea es bastante aceptable desde el punto de vista práctico, pero si lo que se desea es saber el valor que está resistiendo la fibra de la membrana tal como ha sido construida, la idea anterior no sirve de mucho.

Para solucionar este problema se hace uso de la asignación de una dirección de fibras a cada elemento como se explicó en el apartado 4, y a partir de ahí se puede hacer un postproceso adecuado, no importa cuantos cambios de dirección llegue a sufrir la membrana. Solo hay que tener cuidado es asignar las direcciones principales a los elementos apropiados. Esta es una tarea que se hace desde el preproceso, o en su defecto, como una parte inicial del cálculo.

Esta idea de postprocesar es exactamente la misma que se sigue para poder hacer el cálculo ortótropo de membranas. Incluso la orientación de las fibras principales se debe hacer aunque se trate de un material isótropo.

10. Ejemplos

Se tratan varios ejemplos con el fin de validar la formulación de la membrana sin ninguna modificación, la membrana junto con los dos modelos de wrinkling y la membrana pretensada. Además se presentas ejemplos con el fin de que se pueda ver una de las aplicaciones de la formulación.

10.1 Airbag

Para comenzar con la validación de la membrana se realiza el análisis de un airbag cuadrado que tiene por lado ${\textstyle l=1.0\times 1.0{\;ft}}$ y su espesor es ${\textstyle h=0.0001{\;ft}}$ Las propiedades del material son

Densidad ${\textstyle \rho =1500{\;lb}{/ft}^{3}}$

Módulos de Young ${\textstyle E=4.32\times 10^{6}{\;\;lb/ft}^{2}}$

Coeficientes de Poisson ${\textstyle \upsilon =0.3}$

Por condiciones de simetría se analiza solamente un octavo de la geometría con una malla estructurada de 10x10 (200 elementos) usando triángulos lineales. Los resultados son comparados con el estudio presentado en la referencia [7]. La deformada de la estructura se muestra en la Fig. 7.

El desplazamiento máximo transversal es de 0.097 in, mientras que la presentada por Accorsi es de 0.100 in. El máximo desplazamiento en el plano es de 0.066 in y el obtenido en la referencia es de 0.077 in.

Hay que tomar en cuenta que la diferencia entre los resultados es debido a que en la referencia se utilizan elementos cuadriláteros de 9 nodos. A pesar de esta gran diferencia en tiempo de cálculo y de memoria requerida, se puede ver que ambas soluciones son bastante parecidas.

10.2 Airbag con wrinkling. Validación

Se resuelve ahora el mismo problema que en el apartado anterior, pero en la formulación se toma en cuenta el wrinkling, con lo que las tensiones de compresión desaparecen de la estructura. La nueva deformada se presenta en la Fig. 8.

Fig. 8.-Original y deformada

El desplazamiento máximo transversal es de 0.250 in, mientras que la presentada por Accorsi es de 0.252 in. El máximo desplazamiento en el plano es de 0.142 in y el obtenido en la referencia es de 0.144 in.

Tomando en cuenta de que la estructura analizada en la referencia tiene más grados de libertad, se hace un segundo análisis igualando dicho número y se obtienen los siguientes resultados. El desplazamiento máximo transversal es de 0.251 in y el máximo desplazamiento en el plano es de 0.1442 in. Con esto se demuestra el buen comportamiento del elemento con wrinkling.

10.3 Membrana pretensada. Validación

Cuando una membrana es utilizada como una estructura que deberá tener una configuración inicial determinada, lo más común es que dicha forma se adquiera gracias a unas tensiones de pretensado. Para hacer una validación de nuestra formulación incluyendo el pretensado se analiza la estructura de la referencia [8] y [9] de Levy y Gil respectivamente. La estructura tiene por lado ${\textstyle l=240\times 240{\;in}}$ y su espesor es ${\textstyle h=0.004167{\;in}}$ Las propiedades del material son

Densidad ${\textstyle \rho =500{\;lb}{/in}^{3}}$

Módulos de Young ${\textstyle E=30.0\times 10^{6}{\;\;psi}}$

Coeficientes de Poisson ${\textstyle \upsilon =0.3}$

La estructura es sometida completamente a unas tensiones de presfuerzo de 80,000 psi en x e y, sin tensiones cortantes de presfuerzo. La carga transversal es de –10,000 lb en la parte central de la malla. La estructura tiene la malla mostrada en la Fig. 9.

Fig. 9.-Malla y puntos de validación

Los resultados para comparar son

${\begin{array}{cccccccccc}{}&{}&{Levy}&{}&{}&{Gil}&{}&{}&{Valdes}&{}\\{Nodo}&{u}&{v}&{w}&{u}&{v}&{w}&{u}&{v}&{w}\\{4}&{0.015}&{-0.015}&{-1.431}&{0.014}&{-0.014}&{-1.423}&{0.014}&{-0.014}&{-1.429}\\{8}&{0.000}&{-0.017}&{-2.605}&{0.000}&{-0.017}&{-2.600}&{0.000}&{-0.017}&{-2.600}\\{9}&{0.000}&{0.000}&{-6.642}&{0.000}&{0.000}&{-6.626}&{0.000}&{0.000}&{-6.626}\end{array}}$

Gráficamente los resultados anteriores se muestran en la Fig. 10.

Fig. 10.-Desplazamientos x, y, z

Como se puede ver, los resultados de los anteriores trabajos son prácticamente los mismo, por lo que se da como superada nuestra validación. La convergencia en este tipo de casos no presentó ningún problema conservándose siempre cuadrática.

10.4 Paracaídas

Por último se hace el análisis de un paracaídas con la formulación sin wrinkling y también tomando en cuenta el wrinkling. Los datos se toman directamente de la referencia [7], donde las propiedades del material son

Densidad ${\textstyle \rho =0.6{\;lb}{/ft}^{3}}$

Módulos de Young ${\textstyle E=4.32\times 10^{6}{\;\;lb/ft}^{2}}$

Coeficientes de Poisson ${\textstyle \upsilon =0.3}$

El espesor de la membrana es 0.0001 ft y se ha tenido una presión interna de ${\textstyle 0.5{\;lb/ft}^{2}}$. La malla utilizada para el análisis se muestra en la Fig. 11. Ya que en la referencia no se da la geometría del paracaídas, simplemente se trató de hacer un análisis con una geometría similar y ver cuales eran los resultados.

Fig. 11.-Malla utilizada en el análisis

El primer análisis se llevó a cabo sin utilizar el wrinkling en la estructura, lo que indica que existían tensiones de compresión en la estructura. Los resultados de este análisis se muestran en la Fig.12.

Fig. 12.-Malla deformada en distintos pasos de tiempo

Como podemos ver en la figura anterior, el análisis se está llevando a cabo tal y como se esperaba, sin ninguna complicación, excepto que en el análisis existen tensiones de compresión en el material, lo cual no ocurre en la realidad.

Debido a que en el centro de la estructura cada vez aumentan más las tensiones de compresión, aparecen los problemas que se muestran en la Fig. 13 que hacen que el análisis tome un rumbo inesperado y nos aleje demasiado de la solución buscada.

Fig. 13.-Malla deformada en distintos pasos de tiempo

Si este paracaídas se hubiera construido con los materiales que tiene nuestro análisis y tomaran en cuenta las tensiones de compresión, las personas encargadas de probar su funcionamiento simplemente hubieran pasado a mejor vida.

La solución en este análisis tenía un paso de tiempo adecuado para su análisis y llevaba las iteraciones necesarias de un cálculo normal, sin complicación alguna hasta que apareció el efecto negativo que obtuvimos.

También se realizó un segundo análisis, pero esta vez adicionando a la formulación el wrinkling, lo cual produjo los resultados de la Fig. 14.

Fig. 14.-Malla deformada en distintos pasos de tiempo

Como podemos ver este resultado es exactamente lo que buscábamos simular. El inconveniente en este análisis es la gran cantidad de iteraciones que se necesitan al principio en cada paso de tiempo para converger, lo cual tiende a desaparecer al final del análisis.

El segundo inconveniente es que se necesitan muchísimos pasos de tiempo para llegar a la solución. Se necesitarán hacer más pruebas y compararlas con algunos resultados numéricos para saber si ese es el tiempo que se requiere para proceder a realizar dicho análisis antes de dar por satisfactorio este resultado, que sin duda es muy bueno.

En la Fig. 15 se muestra el resultado de los desplazamientos máximos alcanzados por el paracaídas.

11. Referencias

[1] Taylor R. L. Finite element analysis of membrane structures. CIMNE internal report. (2001).

[2] Zienkiewicz O. C. and Taylor R. L. The finite element method. McGraw-Hill 4${\textstyle {}^{th}}$ edition. Volumes 1 and 2. (2002).

[3] Oñate E. Cálculo de estructuras por el método de los elementos finitos. CIMNE 1${\textstyle {}^{st}}$ edición. (1992).

[4] Valdés, G. et al, Assignment of principal fibers direction to analyze orthotropic membranes. GiD 2004 2nd Conference on Advances and Applications of GiD. J. Mora et al Eds., pags. 105-109. CIMNE 2004.