Improving a strongly coupled method for FSI by a simple approximation of the pressure tangent matrix

Resumen

La interacción de los fluidos con las estructuras de su entorno es un desafío clásico para las técnicas numéricas. El objetivo de este trabajo es doble: en primer lugar se proporciona una explicación teórica simple de los principales problemas que se deben superar cuando se trata con un fluido incompresible. A continuación se introduce y justifica un nuevo procedimiento para la solución de problemas complejos de interacción fluido-estructura. Dicha estrategia se basa en la introducción de un «laplaciano de interfase» en el contorno común entre ambos medios. La idea es considerar la dependencia entre la presión del fluido y la velocidad de la estructura como un problema no lineal, que se va a resolver mediante un esquema cuasi-Newton. Se demuestra que el término de interfase resultante es una aproximación de la matriz tangente de dicho problema no lineal, usando exclusivamente álgebra lineal discreta. Finalmente, se verifica la validez de la técnica propuesta mediante su aplicación a algunos ejemplos.

Abstract

The interaction of fluids with surrounding structures constitutes a classical challenge for the different numerical techniques. The aim of current work is twofold: first we provide a simple theoretical explanation of the problems to be faced in incompressible FSI. Then we introduce and justify a new procedure for the solution of complex fluid-structure interaction problems. Such a new strategy is based on the introduction of an «interface Laplacian» at the coupling boundary. The idea is to consider the dependence between fluid pressure and structural velocity as a non linear problem for which a Quasi-Newton scheme is sought. The new interface term is then proved to be an approximation of the tangent matrix for such non-linear problem. In the derivation of this result we make use exclusively of discrete linear algebra. Finally, we prove the efficiency of the new approach showing its ability to tackle standard benchmark problems.

Palabras clave

Interacción fluido-estructura ; Acoplamiento fuerte ; Método particionado

Keywords

Fluid-structure interaction ; Strong coupling ; Partitioned method

1. Introducción

La simulación de la interacción entre estructuras flexibles y fluidos (FSI por la siglas en inglés de Fluid-Structure Interaction ) representa una de las áreas de investigación más activas en el campo de los métodos numéricos. A lo largo de los años, un gran número de autores, procedentes de múltiples campos de la ingeniería, se han enfrentado a este problema, como se puede ver, por ejemplo, en [1] , donde se hace una revisión de las principales técnicas de solución segregada para FSI. Debido a su importancia inherente para el sector aeroespacial, un número de técnicas FSI se diseñaron con el objetivo de simular problemas aeroelásticos, para los cuales se usaron modelos compresibles [2] . Se desarrollaron múltiples algoritmos débilmente acoplados y se estudió su estabilidad para su aplicación tanto en problemas simples de prueba como en estructuras de ingeniería complejas [3]  and [4] . La importancia de la llamada ley de conservación geométrica también fue destacada y puesta en relación con la precisión de la conservación de energía para el problema acoplado [5] , [6]  and [7] .

Múltiples autores intentaron aplicar tales métodos de acoplamiento a problemas que incluyen la interacción de estructuras con fluidos incompresibles, pero rápidamente se demostró que los algoritmos débilmente acoplados eran altamente inestables, con lo que se vio claramente la necesidad de utilizar esquemas de solución fuertemente acoplados. Los primeros resultados en esta dirección se obtuvieron por iteración simple entre los distintos dominios, mientras que la convergencia se garantizaba mediante un factor de subrelajación. A la vez, se desarrollaron estrategias de aceleración que eligen factores de subrelajación «óptimos». La más común de estas técnicas es el acelerador de Aitken, propuesto por primera vez en [8] y publicado posteriormente en [9] . En [10] se presentan técnicas alternativas para obtener el mismo resultado. Una de las características más importantes de este tipo de métodos acoplados es su gran modularidad, lo cual permite utilizar distintas estrategias de solución en cada campo como «cajas negras». Esquemas de proyección estándares como [11] , [12]  and [13] se usan frecuentemente para aprovechar su eficiencia computacional. Como ejemplo, la aplicación de un esquema de paso fraccionado para la solución de problemas acoplados se presenta en [14]  and [15] , tanto en forma débilmente acoplada como dentro de un esquema iterativo encaminado hacia una técnica fuertemente acoplada.

Solamente en tiempos recientes se empezó a ver el problema FSI como un problema acoplado en el cual la estructura y el fluido se presentan como partes de una única discretización. Este enfoque presenta claras ventajas sobre anteriores formulaciones ya que, en particular, permite realizar un análisis matemático del conjunto del proceso acoplado.

Este concepto dio lugar a la aparición tanto de técnicas monolíticas para la solución del problema FSI (véase por ejemplo [16] ) como de esquemas semi-explícitos basados en la aplicación del concepto de separación de la presión partiendo del problema monolítico. Tales esquemas semi-implícitos fuertemente acoplados se trataron por primera vez en [17] . Un posterior desarrollo de dicha idea se presenta en [18] y [19] . El papel que juega la estabilización en este campo se destaca de forma muy clara en [6] y [20] . En esta misma tesis se presenta además una interesante discusión sobre la importancia de la ley de conservación geométrica dependiendo de la elección de la forma conservativa o no conservativa en la implementación del algoritmo de solución de la ecuación de Navier-Stokes.

Un paradigma de acoplamiento alternativo se basa en la llamada condición Robin-Robin y se presenta en [21] , mientras que una formulación distinta basada en el uso de una proyección modificada en la interfase se presenta en [22] .

El objetivo de este artículo es presentar una forma de acoplamiento cuasi-implícita que sea adecuada para problemas de interés en ingeniería. Se pretende proporcionar una justificación tan simple como sea posible de las fuentes de los problemas encontrados con frecuencia en el enfoque desacoplado de FSI, proponiendo un remedio que elimine en la práctica el efecto de masa añadida.

Con la intención de hacer esta presentación más clara el método de acoplamiento se introduce como una iteración de Newton aproximada, lo cual permite justificar la inclusión de no linealidades a todos los niveles y conduce de forma natural a la obtención del esquema final.

El algoritmo que se presenta está enfocado a un procedimiento de solución que sea «menos acoplado» que las alternativas presentadas en trabajos relacionados y presenta como ventaja una completa modularidad de la solución de los dominios fluido y estructural. En esta forma se permite el uso de métodos de solución directos para el dominio estructural.

Para lograr este resultado se encontrará una aproximación simple para el operador tangente del problema no lineal y se justificará su uso en el marco de un algoritmo de acoplamiento de tipo de paso fraccionado. Se demostrará que esto se puede conseguir con una modificación simple del operador laplaciano de la ecuación de presión.

Con el uso del esquema propuesto se pretende obtener las siguientes mejoras:

• Permitir incorporar no linealidades en la solución de la estructura en las iteraciones del proceso de acoplamiento. Como consecuencia, las no linealidades que pueden aparecer en la estructura (como por ejemplo daño) no afectarán directamente el número de iteraciones del problema acoplado.
• Dar la posibilidad de resolver la estructura mediante esquemas de solución directos. Esto es fundamental por ejemplo en la solución de láminas delgadas, que frecuentemente incluyen rotaciones u otros grados de libertad no estándar.
• Poder utilizar de forma natural mallas independientes para cada dominio, sin necesidad de que concuerden en la interfase, siendo necesario únicamente proporcionar un mecanismo de transmisión de las variables adecuado.
• Conseguir resolver problemas en los que el fluido queda totalmente confinado (sin condiciones de contorno de presión en ningún punto) sin tener que introducir modificaciones en la solución de la estructura.
• Se permite utilizar el esquema propuesto usando directamente un código estructural existente, aunque será necesario introducir una pequeña modificación en el código del fluido.
• Dejar la posibilidad de usar diferentes esquemas de integración temporal en la estructura y el fluido. Sin embargo, en este contexto se debe tener cuidado para asegurar la continuidad de desplazamientos y aplicar condiciones de contorno de velocidad apropiadas en la interfase.

Por otro lado, la aplicabilidad del método propuesto queda limitada al uso de interpolaciones de igual orden para velocidad y presión, puesto que la formulación se plantea a partir de la presencia de estabilización para la presión. El desarrollo se hace aquí en el contexto de métodos de tipo paso fraccionado, aunque algunos de los resultados obtenidos se pueden extender a un contexto más amplio.

Teniendo esto presente, el artículo se organiza en la siguiente forma:

• Descripción del problema físico.
• Justificación del esquema FSI segregado aplicado al problema completo.
• Derivación de un nuevo algoritmo cuasi-implícito.
• Resolución de ejemplos de aplicación.

Todos los ejemplos presentados en el presente artículo se han resuelto usando Kratos , un código de solución multipropósito que se está desarrollando en CIMNE [23]  and [24] . Más información sobre el mismo se encuentra disponible en la dirección http://kratos.cimne.upc.es/kratoswiki/index.php/Main_Page .

2. Planteamiento del problema físico

El punto de partida de la formulación del problema acoplado es el equilibrio dinámico (conservación del momento lineal), expresado en forma diferencial como

$\rho {\frac {D{\boldsymbol {\mbox{v}}}}{Dt}}-\nabla \cdot {\boldsymbol {\sigma }}=$${\boldsymbol {\mbox{f}}}$

( 1)

En el contexto de los elementos finitos esta ecuación se escribe sobre la base de una discretización espacial. Para su uso en la práctica es conveniente expresar la derivada total explícitamente, tomando en cuenta la posibilidad de una deformación dinámica de los elementos finitos que constituyen la discretización. Esto se consigue mediante la denominada descripción Arbitraria Lagrangiana-Euleriana (ALE, véase [7] ), que permite reescribir la conservación de momento como

$\rho {\frac {\partial {\boldsymbol {\mbox{v}}}}{\partial t}}+$$\rho \left({\boldsymbol {\mbox{v}}}-{\boldsymbol {\mbox{v}}}_{M}\right)\cdot \nabla {\boldsymbol {\mbox{v}}}-$$\nabla \cdot {\boldsymbol {\sigma }}={\boldsymbol {\mbox{f}}}$

( 2)

donde v_M representa la velocidad de la malla, definida como la velocidad con que se mueven los nodos de la malla por el espacio, y ${\textstyle {\frac {\partial {\boldsymbol {\mbox{v}}}}{\partial t}}}$ es la derivada de las velocidades nodales calculadas en los puntos de la malla ALE en movimiento. Esta forma permite tener en cuenta descripciones tanto eulerianas como lagrangianas, además de cualquier posibilidad intermedia.

La Ecuación (2) deriva de principios básicos y no depende de ninguna hipótesis sobre cómo una parte del material se va a comportar ante cargas externas. Esto justifica su validez tanto en el dominio sólido como en el fluido. La elección del material se realiza a posteriori prescribiendo la relación o ley constitutiva entre las tensiones y una medida de las deformaciones o la velocidad de deformación.

Desde este punto de vista el acoplamiento entre un fluido y una estructura no presenta ningún problema conceptual, puesto que se puede definir una formulación común que permita resolver ambos dominios a la vez. Este sistema es el que se utiliza por ejemplo en el Método de Elementos Finitos y Partículas (PFEM) al tomar una discretización monolítica del sistema acoplado [25] o una fraccionada como en [26] .

Esto no obstante, cabe destacar que el procedimiento más habitual es utilizar descripciones distintas para cada dominio, siendo la euleriana preferida para el fluido y la lagrangiana la más habitual para describir el dominio estructural.

En general, el caso incompresible se trata mediante la introducción de la restricción ∇ · v  = 0 en el campo de velocidades (véanse por ejemplo [7]  and [13] ). El sistema de ecuaciones resultante está bien planteado en el continuo pero presenta problemas en su forma discreta a no ser que se elijan cuidadosamente los espacios de velocidad y presión (véase por ejemplo [7] ). En la práctica se prefiere utilizar espacios de interpolación del mismo orden para velocidad y presión, para lo cual es necesario introducir un término de estabilización de presión. En este artículo se asume en todo momento que existe un término de estabilización de este tipo y que se usa una interpolación de igual orden para velocidades y presiones. Como consecuencia, la formulación aquí presentada no será directamente aplicable a elementos basados en interpolaciones estables.

Las ecuaciones

${\boldsymbol {\mbox{r}}}^{f}\left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)=$${\boldsymbol {\mbox{f}}}^{f}-{\boldsymbol {\mbox{K}}}^{_{\ast }}\left({\boldsymbol {\mbox{v}}}-\right.$$\left.{\boldsymbol {\mbox{v}}}_{M}\right){\boldsymbol {\mbox{v}}}-$${\boldsymbol {\mbox{G}}}{\boldsymbol {\mbox{p}}}-\rho {\boldsymbol {\mbox{M}}}^{f}{\frac {\partial {\boldsymbol {\mbox{v}}}}{\partial t}}$

( 3)

y

$\Psi \left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)=$${\boldsymbol {\mbox{D}}}{\boldsymbol {\mbox{v}}}+{\frac {\tau }{\rho }}{\boldsymbol {\mbox{S}}}{\boldsymbol {\mbox{p}}}$

( 4)

representan la forma residual de las ecuaciones discretas de conservación de momento y de masa. El término ${\textstyle {\frac {\tau }{\rho }}{\boldsymbol {\mbox{S}}}}$ representa una estabilización adecuada de la ecuación de conservación de masa, siendo τ un factor escalar de estabilización. v y p representan los valores de velocidad y presión en el instante de tiempo n  + 1. Las matrices y vectores de las ecuaciones (3) y (4) se expresan siguiendo una notación estándar que se presenta por ejemplo en [7] o [27] . Se considera que cada una de las matrices se construye a partir de la discretización de Galerkin de los correspondientes operadores diferenciales, además de términos de estabilización, que van a depender del método se estabilización empleado. Si se consideran las funciones de test w , q las matrices se van a construir a partir de la discretización de los siguientes términos

${\boldsymbol {\mbox{K}}}^{\ast }\leftarrow ({\boldsymbol {\mbox{w}}},\rho ({\boldsymbol {\mbox{v}}}-$${\boldsymbol {\mbox{v}}}_{M})\cdot \nabla ))+({\boldsymbol {\mbox{w}}},\mu {\mbox{ }}\Delta {\mbox{ }}{\boldsymbol {\mbox{v}}})+$${\boldsymbol {\mbox{S}}}_{K}$

( 5)

${\boldsymbol {\mbox{G}}}\leftarrow ({\boldsymbol {\mbox{w}}},\nabla p)+$${\boldsymbol {\mbox{S}}}_{D}$

( 6)

${\boldsymbol {\mbox{D}}}\leftarrow \left(q,,,\nabla \cdot {\boldsymbol {\mbox{v}}}\right)+$${\boldsymbol {\mbox{S}}}_{G}$

( 7)

Las matrices S , S_K , S_G y S_D contienen los términos de estabilización y dependen del método elegido. Por ejemplo, para ASGS [28] , asumiendo funciones de interpolación lineales:

${\boldsymbol {\mbox{S}}}\leftarrow \left(\nabla q{\mbox{,}}\nabla p\right)$

( 8)

${\boldsymbol {\mbox{S}}}_{K}\leftarrow \tau \left(\left({\boldsymbol {\mbox{v}}}-\right.\right.$$\left.\left.{\boldsymbol {\mbox{v}}}_{M}\right)\cdot \nabla {\boldsymbol {\mbox{w}}}{\mbox{,}}\rho {\boldsymbol {\mbox{M}}}^{f}{\frac {\partial {\boldsymbol {\mbox{v}}}}{\partial t}}+\right.$$\left.\rho \left({\boldsymbol {\mbox{v}}}-{\boldsymbol {\mbox{v}}}_{M}\right)\cdot \nabla {\boldsymbol {\mbox{v}}}\right)+$${\overset {\mbox{ˆ}}{\tau }}\left(\nabla \cdot {\boldsymbol {\mbox{w}}}{\mbox{,}}\nabla \cdot {\boldsymbol {\mbox{v}}}\right)$

( 9)

${\boldsymbol {\mbox{S}}}_{G}\leftarrow \tau \left(\left({\boldsymbol {\mbox{v}}}-\right.\right.$$\left.\left.{\boldsymbol {\mbox{v}}}_{M}\right)\cdot \nabla {\boldsymbol {\mbox{w}}}{\mbox{,}}\nabla p\right)$

( 10)

${\boldsymbol {\mbox{S}}}_{D}\leftarrow \tau \left(\nabla q{\mbox{,}}{\boldsymbol {\mbox{M}}}^{f}{\frac {\partial {\boldsymbol {\mbox{v}}}}{\partial t}}+\right.$$\left.\left({\boldsymbol {\mbox{v}}}-{\boldsymbol {\mbox{v}}}_{M}\right)\cdot \nabla {\boldsymbol {\mbox{v}}}\right)$

( 11)

Donde ${\textstyle {\overset {\mbox{ˆ}}{\tau }}}$ representa un segundo parámetro de estabilización y (a , b ) el producto interno entre los vectores a y b . Sin embargo, se quiere destacar que no es necesario el uso de estabilización ASGS para formulación que se propone en este artículo. Lo fundamental en este caso es la presencia de la matriz S , que no aparece en el método de Galerkin pero es habitual en la mayoría de formulaciones estabilizadas, y que va a jugar un papel importante en nuestro método.

Para completar la presentación se va a introducir una formulación discreta análoga para el dominio estructural. Siempre que la estructura sea compresible se puede usar una ley constitutiva para relacionar tensiones y deformaciones. Bajo esta condición es posible describir completamente el comportamiento estructural mediante el uso de la Ecuación (2) , lo cual conduce a la forma residual

${\boldsymbol {\mbox{r}}}^{s}\left({\boldsymbol {\mbox{v}}}\right)=$${\boldsymbol {\mbox{f}}}^{s}-{\boldsymbol {\mbox{K}}}^{s}\left({\boldsymbol {\mbox{u}}}\right){\boldsymbol {\mbox{u}}}-$$\rho {\boldsymbol {\mbox{M}}}^{s}{\frac {\partial {\boldsymbol {\mbox{v}}}}{\partial t}}$

( 12)

donde un proceso de discretización equivalente al del fluido resultará en

$K^{s}\left({\boldsymbol {\mbox{u}}}\right){\boldsymbol {\mbox{u}}}\leftarrow \left({\boldsymbol {\mbox{w}}},,,\nabla \cdot \sigma \right)$

( 13)

siendo ${\textstyle {\boldsymbol {\mbox{u}}}={\boldsymbol {\mbox{u}}}\left({\boldsymbol {\mbox{v}}}\right)}$ el vector de desplazamientos nodales, relacionado con las velocidades v por el esquema de integración temporal elegido.

2.1. Formulación FSI

En este punto centramos nuestra atención en el problema representado en la figura 1 . El equilibrio de fuerzas y la continuidad de desplazamientos se debe verificar en el contorno entre el fluido y la estructura.

Figura 1. Discretización monolítica de los dominios del fluido y la estructura.

Para tal discretización los nodos que quedan en la interfase de acoplamiento pertenecen simultáneamente a ambos dominios. Puesto que su posición en el espacio es única, una discretización de este tipo garantiza que los desplazamientos y velocidades serán continuas en todo el dominio.

En el dominio estructural los desplazamientos y velocidades están fuertemente relacionadas entre sí y la velocidad de la malla coincide con la del cuerpo. Este no es el caso en el dominio fluido, donde la deformación de la malla se fija independientemente del campo de velocidades euleriano. Por lo tanto, se va a asumir explícitamente que el desplazamiento en la interfase se refleja instantáneamente en la deformación del dominio fluido o, en otras palabras, se asumirá una relación dada ${\textstyle \Delta {\boldsymbol {\mbox{x}}}=\Delta {\boldsymbol {\mbox{x}}}\left({\boldsymbol {\mbox{v}}}_{fs}\right)}$ y ${\textstyle {\boldsymbol {\mbox{v}}}_{M}={\boldsymbol {\mbox{v}}}_{M}\left(\Delta {\boldsymbol {\mbox{x}}}\right)}$ que representa la acción del algoritmo de deformación de la malla.

La condición de equilibrio se puede expresar como la imposición de que el residuo de la ecuación de momento sea cero.

${\boldsymbol {\mbox{r}}}\left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)=$${\boldsymbol {0}}$

( 14)

Como se ha visto, el residuo total en el conjunto del dominio se puede expresar mediante la combinación de las ecuaciones (3) y (12) . La malla de la figura 1 se ha definido de forma que cada elemento de la malla sea exclusivamente un elemento de fluido o un elemento estructural. Por lo tanto, en cada elemento estructural el residuo es gobernado por la ecuación (12) , mientras que cada elemento fluido se rige por la ecuación (3) . Como algunos nodos son compartidos entre ambos dominios, en sus correspondientes posiciones en la matriz del sistema se van a sumar contribuciones de ambos dominios.

Antes de continuar con la exposición de la formulación propuesta, es conveniente aclarar la notación que vamos a utilizar. Los residuos introducidos hasta ahora se van a caracterizar en la práctica por sus valores nodales. Puesto cada término del residuo está asociado a un grado de libertad de un nodo de la malla, las componentes del residuo se pueden organizar según la posición del nodo al que pertenecen. Tal y como se observa en la figura 1 , se van a agrupar los términos pertanecientes al interior del dominio fluido, la interfase y el interior del dominio sólido, que se identifican con los subíndices f , fs y s respectivamente. Esto nos permite escribir simbólicamente

${\boldsymbol {\mbox{r}}}\quad :=\left({\begin{array}{c}{\boldsymbol {\mbox{r}}}_{f}\\{\boldsymbol {\mbox{r}}}_{fs}\\{\boldsymbol {\mbox{r}}}_{s}\end{array}}\right)=\left({\begin{array}{c}{\boldsymbol {\mbox{r}}}_{f}\\{\boldsymbol {\mbox{r}}}_{fs}^{f}\\{\boldsymbol {0}}\end{array}}\right)+\left({\begin{array}{c}{\boldsymbol {0}}\\{\boldsymbol {\mbox{r}}}_{fs}^{s}\\{\boldsymbol {\mbox{r}}}_{s}\end{array}}\right)={\boldsymbol {\mbox{r}}}^{s}\oplus {\boldsymbol {\mbox{r}}}^{f}$

( 15)

donde el símbolo ⊕ denota el ensamblaje de elementos finitos y

${\boldsymbol {\mbox{r}}}^{f}\quad :=\left({\begin{array}{c}{\boldsymbol {\mbox{r}}}_{f}\\{\boldsymbol {\mbox{r}}}_{fs}^{f}\\{\boldsymbol {0}}\end{array}}\right){\mbox{;}}\quad {\boldsymbol {\mbox{r}}}^{s}\quad :=$$\left({\begin{array}{c}{\boldsymbol {0}}\\{\boldsymbol {\mbox{r}}}_{fs}^{s}\\{\boldsymbol {\mbox{r}}}_{s}\end{array}}\right)$

( 16)

Para continuar con el razonamiento se escriben los residuos explícitamente como

$\left({\begin{array}{c}{\boldsymbol {\mbox{r}}}_{f}\left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)\\{\boldsymbol {\mbox{r}}}_{fs}\left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)\\{\boldsymbol {\mbox{r}}}_{s}\left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)\end{array}}\right):=\left({\begin{array}{c}{\boldsymbol {\mbox{f}}}_{f}\\{\boldsymbol {\mbox{f}}}_{fs}^{f}+{\boldsymbol {\mbox{f}}}_{fs}^{s}\\{\boldsymbol {\mbox{f}}}_{s}\end{array}}\right)-\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\\{\boldsymbol {0}}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{p}}}\end{array}}\right)$ $-\left({\begin{array}{ccc}{\rho }_{f}{\boldsymbol {\mbox{M}}}_{f}&{\boldsymbol {0}}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\rho }_{f}{\boldsymbol {\mbox{M}}}_{fs}^{f}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {0}}&{\boldsymbol {0}}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{a}}}_{f}\\{\boldsymbol {\mbox{a}}}_{fs}\\{\boldsymbol {\mbox{a}}}_{s}\end{array}}\right)-\left({\begin{array}{ccc}{\boldsymbol {\mbox{K}}}_{11}^{_{\ast }}&{\boldsymbol {\mbox{K}}}_{12}^{_{\ast }}&{\boldsymbol {0}}\\{\boldsymbol {\mbox{K}}}_{21}^{_{\ast }}&{\boldsymbol {\mbox{K}}}_{22}^{_{\ast }}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {0}}&{\boldsymbol {0}}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{v}}}_{f}\\{\boldsymbol {\mbox{v}}}_{fs}\\{\boldsymbol {\mbox{v}}}_{s}\end{array}}\right)$ $-\left({\begin{array}{ccc}{\boldsymbol {0}}&{\boldsymbol {0}}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\rho }_{s}{\boldsymbol {\mbox{M}}}_{fs}^{s}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {0}}&{\rho }_{s}{\boldsymbol {\mbox{M}}}_{s}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{a}}}_{f}\\{\boldsymbol {\mbox{a}}}_{fs}\\{\boldsymbol {\mbox{a}}}_{s}\end{array}}\right)-\left({\begin{array}{ccc}{\boldsymbol {0}}&{\boldsymbol {0}}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {\mbox{K}}}_{11}^{s}&{\boldsymbol {\mbox{K}}}_{12}^{s}\\{\boldsymbol {0}}&{\boldsymbol {\mbox{K}}}_{21}^{s}&{\boldsymbol {\mbox{K}}}_{22}^{s}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{u}}}_{f}\\{\boldsymbol {\mbox{u}}}_{fs}\\{\boldsymbol {\mbox{u}}}_{s}\end{array}}\right)$

( 17)

donde la nueva variable a indica la aceleración. Los términos ${\textstyle {\boldsymbol {\mbox{M}}}_{fs}^{s}}$ y ${\textstyle {\boldsymbol {\mbox{M}}}_{fs}^{f}}$ siguen la notación introducida para los residuos y representan, respectivamente, la contribución de la estructura y del fluido a los términos de la matriz de masas que corresponden a la interfase FSI. Se asume que la matriz de masa es aglutinada (lumped ) en ambos dominios.

La ecuación de los residuos se debe completar con la condición de incompresibilidad descrita en la ecuación (4) , que proporciona una relación entre las velocidades de todos los nodos en el dominio fluido.

${\boldsymbol {\Psi }}\left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)=$$\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{v}}}_{f}\\{\boldsymbol {\mbox{v}}}_{fs}\end{array}}\right)+{\frac {\tau }{\rho }}{\boldsymbol {\mbox{S}}}\quad \left({\begin{array}{c}{\boldsymbol {\mbox{p}}}\end{array}}\right)={\boldsymbol {0}}$

( 18)

El problema acoplado no lineal se puede resolver ahora mediante la adecuada linealización de las ecuaciones (17) y (18) , que se puede expresar simbólicamente como

$\left({\begin{array}{c}-{\frac {\partial {\boldsymbol {\mbox{r}}}}{\partial {\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}}}\\-{\frac {\partial {\boldsymbol {\Psi }}}{\partial {\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}}}\end{array}}\right)\left({\begin{array}{c}d{\boldsymbol {\mbox{v}}}\\d{\boldsymbol {\mbox{p}}}\end{array}}\right)=\left({\begin{array}{c}{\boldsymbol {\mbox{r}}}\left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)\\{\boldsymbol {\Psi }}\left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)\end{array}}\right)$

( 19)

La construcción de tal sistema no entraña ninguna dificultad especial si se dispone de una formulación en velocidades para el fluido y la estructura en un mismo código. Es asimismo destacable que la solución obtenida de este modo es óptima en el sentido de que proporciona las mismas garantías teóricas que una formulación de elementos finitos estándar de un solo dominio.

Sin embargo es necesario pagar un precio para obtener tal resultado: la dificultad de resolver el problema se transforma en la de resolver el sistema lineal resultante. Esto es fácil de ver si se considera que el sistema que se debe resolver en cada iteración

• Es no simétrico debido a la presencia del operador convectivo.
• Está dotado de rigideces «viscosas» muy distintas en las partes estructurales y fluidas.
• Trata con variables que en general están escaladas de forma muy distinta.

Además, la inclusión de mallas no concordantes en la interfase de acoplamiento, de vital importancia en problemas reales, no se puede hacer de forma trivial con esta formulación.

A continuación se partirá del problema monolítico para derivar nuestra formulación acoplada mediante una técnica clásica de segregación, siguiendo el enfoque algebraico presentado en [29] .

Antes de proceder se hace necesario destacar dos puntos importantes. En primer lugar, se debe reconocer que las ecuaciones descritas hasta ahora tienen una dependencia no lineal con la posición de los nodo y/o con la velocidad de la malla. Tal dependencia se debe tener en cuenta cuando se pretende escribir los residuos y el operador tangente de la estructura. Los autores afirmamos que, sin embargo, el efecto del movimiento de la interfase de acoplamiento desde la posición predicha hasta la posición corregida tiene un impacto menor en la evaluación de los operadores tangentes en el dominio fluido. El razonamiento heurístico en el que se basa esta afirmación es que un cambio dv de la velocidad de la interfase (que se comporta como lagrangiana) implica una corrección dx  ≈ dvdt de su posición, que tenderá a cero a medida que dt  → 0. Además, aun en el caso de que esta aproximación no sea válida, se obtendrá igualmente la solución exacta, puesto que los residuos se calculan en todo momento exactamente, viéndose afectada únicamente la velocidad de convergencia.

En segundo lugar, se observa que la Ecuación (19) se ha expresado simbólicamente en términos de velocidades, una formulación que no es la habitual para el dominio estructural. Esto simplifica la justificación del algoritmo, pero no es un requerimiento para su implementación práctica. Todos los ejemplos que se presentan en este artículo han sido calculados usando una formulación estándar en desplazamientos para el sólido. El único requerimiento ineludible es ser capaz de evaluar el operador tangente a partir de la matriz de rigidez estática tangente, la matrices de masas y amortiguamiento y el correspondiente operador dinámico tangente, una capacidad que se requiere comúnmente en cualquier código de cálculo dinámico de estructuras.

3. Esquemas segregados para FSI

A continuación se trata la aplicación de enfoques segregados a la solución del problema FSI acoplado.

Los esquemas de segregación se basan en la idea clave de que la conservación del momento linear se puede expandir en una serie de Taylor alrededor de un valor auxiliar de velocidad ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}$ y un valor conocido de presión p_kn para obtener

${\boldsymbol {\mbox{r}}}\left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)\approx {\boldsymbol {\mbox{r}}}\left({\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}_{kn}\right)+$${\frac {\partial {\boldsymbol {\mbox{r}}}}{\partial {\boldsymbol {\mbox{p}}}}}\Delta {\boldsymbol {\mbox{p}}}+$${\frac {\partial {\boldsymbol {\mbox{r}}}}{\partial {\boldsymbol {\mbox{v}}}}}\left({\boldsymbol {\mbox{v}}}-\right.$$\left.{\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}\right)$

( 20)

con Δp  ≔ p  − p_kn . Claramente tal expansión será aceptable únicamente si Δp  ≈ 0 y ${\textstyle {\boldsymbol {\mbox{v}}}-{\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}\approx 0}$ , es decir, que el enfoque segregado proporcionará una solución adecuada solo para pasos de tiempo suficientemente pequeños (o en régimen estacionario). En la práctica el paso de tiempo se debe elegir de manera que los términos de inercia dominen sobre el resto en el dominio fluido.

En consecuencia, el esquema que se propone se diseñará bajo la suposición de que la inercia sea dominante en el fluido pero, como se verá más adelante, se evitará hacer la misma hipótesis para el sólido.

Continuando con la presentación del método, nos vamos a centrar en el problema expresado por la ecuación (20) , que representa una aproximación al problema original. La ecuación aproximada (20) se puede formular como un sistema (no lineal) de ecuaciones

${\boldsymbol {\mbox{r}}}\left({\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}_{kn}\right)=$$0$

( 21)

y

${\frac {\partial {\boldsymbol {\mbox{r}}}}{\partial {\boldsymbol {\mbox{p}}}}}\Delta {\boldsymbol {\mbox{p}}}+$${\frac {\partial {\boldsymbol {\mbox{r}}}}{\partial {\boldsymbol {\mbox{v}}}}}\left({\boldsymbol {\mbox{v}}}-\right.$$\left.{\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}\right)=$${\boldsymbol {0}}$

( 22)

La ventaja es que ahora la ecuación (22) se puede resolver simbólicamente para explicitar la relación entre la velocidad de fin de paso v y la variable auxiliar ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}$ , normalmente denominada velocidad de paso fraccionado, como

${\boldsymbol {\mbox{v}}}={\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}+$${\left({\frac {\partial {\boldsymbol {\mbox{r}}}}{\partial {\boldsymbol {\mbox{v}}}}}\right)}^{-1}{\frac {\partial {\boldsymbol {\mbox{r}}}}{\partial {\boldsymbol {\mbox{p}}}}}\Delta {\boldsymbol {\mbox{p}}}$

( 23)

Esta relación explícita se puede reemplazar simbólicamente en la ecuación de conservación de la masa para proporcionar una relación entre la corrección de presión y la variable auxiliar ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}$

${\boldsymbol {\Psi }}\left({\boldsymbol {\mbox{v}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)=$$0\rightarrow {\boldsymbol {\Psi }}\left({\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)=$$0$

( 24)

Esta observación nos permite diseñar el siguiente esquema segregado

• Tomar un valor para p_kn . Típicamente p_kn  = p_n para un paso fraccionado puro o ${\textstyle {\boldsymbol {\mbox{p}}}_{kn}={\boldsymbol {\mbox{p}}}_{n+1}^{it-1}}$ para un esquema predictor-corrector.
• Resolver la ecuación (21) para ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}$ . Hasta el momento no se ha introducido ningún tipo de separación, por lo que este paso requeriría una solución monolítica de la velocidad en ambos dominios.
• Obtener la corrección de presión Δp mediante la ecuación (24) .
• Calcular la velocidad de final de paso usando la ecuación (23) .

El concepto básico de este esquema es encontrar primero un valor de la velocidad ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}$ que esté en equilibrio con el valor de presión dado. A continuación se determina una corrección para la presión que garantice que la velocidad final tenga divergencia nula.

Es interesante destacar que hasta este punto se usa una única discretización para el conjunto del dominio. Esto implica que el único error en la solución se debe a la división de las operaciones realizada en la ecuación (20) .

Antes de continuar se proporcionan algunas relaciones explícitas entre las distintas derivadas.

Los residuos totales se obtienen mediante el ensamblaje de las contribuciones de las partes estructurales y fluidas tal como se expone en la ecuación (15) . La definición de la ecuación (12) afirma que la estructura no tiene una dependencia explícita de la presión. De ello se deduce que

${\frac {\partial {\boldsymbol {\mbox{r}}}^{s}}{\partial {\boldsymbol {\mbox{p}}}}}=$${\boldsymbol {0}}$

( 25)

Por otro lado, utilizando la ecuación (3) se obtiene

${\frac {\partial {\boldsymbol {\mbox{r}}}^{f}}{\partial {\boldsymbol {\mbox{p}}}}}=$$\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)$

( 26)

de donde se deduce que

${\frac {\partial {\boldsymbol {\mbox{r}}}}{\partial {\boldsymbol {\mbox{p}}}}}=$${\frac {\partial {\boldsymbol {\mbox{r}}}^{s}}{\partial {\boldsymbol {\mbox{p}}}}}\oplus {\frac {\partial {\boldsymbol {\mbox{r}}}^{f}}{\partial {\boldsymbol {\mbox{p}}}}}=$$\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\\{\boldsymbol {0}}\end{array}}\right)$

( 27)

Un procedimiento similar se puede aplicar para calcular ${\textstyle {\frac {\partial {\boldsymbol {\mbox{r}}}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}}$ como

${\frac {\partial {\boldsymbol {\mbox{r}}}^{f}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}=$${\boldsymbol {\mbox{K}}}^{_{\ast }}\left({\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}-\right.$$\left.{\boldsymbol {\mbox{v}}}_{M}\right)+{\rho }_{f}{\frac {\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{a}}}}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}{\boldsymbol {\mbox{M}}}^{f}$

( 28)

y

${\frac {\partial {\boldsymbol {\mbox{r}}}^{s}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}=$${\frac {\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{u}}}}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}{\boldsymbol {\mbox{K}}}\left({\boldsymbol {\mbox{u}}}\right)+$${\rho }_{s}{\frac {\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{a}}}}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}{\boldsymbol {\mbox{M}}}^{s}$

( 29)

Las derivadas de la aceleración ${\textstyle {\frac {\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{a}}}}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}}$ y la posición ${\textstyle {\frac {\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{u}}}}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}}$ con respeto a la velocidad dependen del esquema de integración temporal utilizado. Para los esquemas de integración más frecuentes se pueden establecer relaciones del tipo

${\frac {\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{a}}}}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}=$${\frac {1}{dt}}c_{a}{\boldsymbol {\mbox{I}}}$

( 30)

y

${\frac {\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{x}}}}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}=$$c_{v}dt{\boldsymbol {\mbox{I}}}$

( 31)

donde c_a y ${\textstyle c_{v}}$ son constantes escalares que dependen del esquema de integración. Para pasos de tiempo suficientemente pequeños los términos dependientes de la masa dominarán en las ecuaciones (28) y (29)

${\frac {\partial {\boldsymbol {\mbox{r}}}^{f}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}{\overset {dt\rightarrow 0}{\rightarrow }}{\frac {c_{a}{\rho }_{f}}{dt}}{\boldsymbol {\mbox{M}}}^{f}$

( 32)

${\frac {\partial {\boldsymbol {\mbox{r}}}^{s}}{\partial {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}}{\overset {dt\rightarrow 0}{\rightarrow }}{\frac {c_{a}{\rho }_{s}}{dt}}{\boldsymbol {\mbox{M}}}^{s}$

( 33)

Esta observación es importante, puesto que los términos descritos aparecen en la ecuación (23) con sus inversas. Desafortunadamente, si se considera que el dominio estructural es típicamente mucho más rígido que el fluido, se debe concluir que el paso de tiempo para el cual la inercia dominará sobre la rigidez será frecuentemente mucho más pequeño para la estructura que para el fluido. Esto implica que cuando se elige el paso de tiempo (como se hace normalmente) para obtener una solución buena (y computacionalmente eficiente) para el dominio fluido la ecuación de la estructura (29) no debería ser aproximada, mientras que la aproximación de la ecuación (24) es aceptable para el fluido. Con esta observación, la conservación de masa, expresada por la ecuación (24) se puede escribir como

$\left(\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {\mbox{A}}}_{11}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {\mbox{A}}}_{22}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)\right)\left({\begin{array}{c}\Delta {\boldsymbol {\mbox{p}}}\end{array}}\right)+{\frac {\tau }{\rho }}{\boldsymbol {\mbox{S}}}\left({\begin{array}{c}{\boldsymbol {\mbox{p}}}\end{array}}\right)=\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{c}{\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}_{f}\\{\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}_{fs}\end{array}}\right)$

( 34)

con

${\boldsymbol {\mbox{A}}}_{11}\quad :={\frac {dt}{c_{a}{\rho }_{f}}}{\boldsymbol {\mbox{M}}}_{f}^{-1}$

( 35)

y

${\boldsymbol {\mbox{A}}}_{22}\quad :={\left({\frac {c_{a}}{dt}}\left({\rho }_{f}{\boldsymbol {\mbox{M}}}_{fs}^{f}+{\rho }_{s}{\boldsymbol {\mbox{M}}}_{fs}^{s}\right)+c_{v}dt{\boldsymbol {\mbox{K}}}^{s}\right)}^{-1}$

( 36)

Al escribir la ecuación (36) se está cometiendo un pequeño abuso de notación: para ser consistente se debería invertir el conjunto de la matriz tangente estructural (incrementada por ${\textstyle {\rho }_{f}{\boldsymbol {\mbox{M}}}_{fs}^{f}}$ en los términos de la interfase) y entonces tomar A₂₂ como los términos de la matriz inversa que corresponden a la interfase de acoplamiento. Una alternativa presentada en [18] es condensar simbólicamente todos los términos en el interior del dominio estructural, lo que permite definir una matriz de la interfase equivalente, cuya inversa se podría usar directamente en la definición de (36) . En cualquier caso esta es una discusión puramente formal, puesto que lo que se pretende es justificar una aproximación de A₂₂ que evite tener que calcularla explícitamente.

Para su posterior uso se especifican también los residuos de la ecuación de conservación de masa:

${\boldsymbol {\mbox{r}}}_{masa}\quad :=\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{c}{\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}_{f}\\{\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}_{fs}\end{array}}\right)-\left(\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {\mbox{A}}}_{11}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {\mbox{A}}}_{22}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)+{\frac {\tau }{\rho }}{\boldsymbol {\mbox{S}}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{p}}}\end{array}}\right)$

( 37)

4. Algoritmo propuesto

El algoritmo para separar la presión como se describió en las secciones precedentes es aplicable a la solución de problemas en los que se utiliza un único dominio que incluye tanto el sólido como el fluido. A continuación se pretende escribir un algoritmo equivalente para un entorno dividido, donde se usan mallas diferentes para la estructura y el fluido. Teniendo esto presente se propone el siguiente algoritmo:

• Resolver la estructura partiendo de p  = p_n .
• Transferir el tamaño aproximado de los elementos, la densidad y propiedades mecánicas de la estructura al fluido.
• Transferir los desplazamientos de la estructura al fluido.
• Calcular la velocidad de la malla e imponer que la velocidad del fluido en el contorno de acoplamiento coincide con esta.
• Deformar la malla de forma acorde a los desplazamientos impuestos.
• Resolver la ecuación de conservación de momento en el fluido exclusivamente.
• Calcular las fuerzas que el fluido ejerce sobre la estructura (usando fuerzas nodales consistentes) y transferirlos a la estructura.
• Resolver la estructura.
• Transferir los desplazamientos de la estructura al fluido.
• Desplazar en el dominio fluido los nodos compartidos con la estructura (la malla no se debe deformar más que en la interfase).
• Resolver la ecuación de presión.
• Corregir las velocidades de acuerdo con la ecuación (23) .
• Verificar la convergencia y volver al paso 7 si es necesario.

En este proceso se distinguen dos bloques de resolución:

• Bucle 1: Los pasos 1 a 8 representan una iteración de punto fijo que resuelve «monolíticamente» la ecuación de momentos partiendo de un valor conocido de la presión.
• Bucle 2: Los pasos 7 a 14 representan un procedimiento iterativo que tiende a la satisfacción exacta de la restricción en la divergencia para el dominio fluido.

Un esquema en la forma descrita se puede implementar fácilmente partiendo de un código existente basado en el esquema de paso fraccionado. Sin embargo, la experiencia demuestra que el algoritmo resultante falla por la falta de convergencia del segundo bucle. Este comportamiento se ha observado anteriormente en el trabajo teórico de [17] , que lo resolvía empleando un esquema iterativo de resolución sin matrices.

El presente artículo pretende exponer un planteamiento alternativo para obtener el mismo resultado. El objetivo de esta formulación es satisfacer las siguientes condiciones:

• Resolver las no linealidades en la solución de la estructura dentro del bucle exterior (bucle 2).
• Permitir el uso de esquemas de solución directos en la estructura, lo cual es fundamental en casos como las láminas delgadas.
• Permitir el uso de mallas no concordantes en la interfase.
• Modificar el problema para el fluido de forma que sea posible resolver problemas completamente confinados (en los que no se imponen condiciones de presión en ningún punto).

Una vez identificado nuestro objetivo, podemos formular el problema en la forma siguiente: una vez encontrada una velocidad de paso fraccionado ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}}$ que resuelva «exactamente» el problema de la ecuación (21) (es decir, suponiendo que el bucle 1 converja satisfactoriamente), encontrar una presión p y una velocidad corregida v_fs en la interfase tal que

${\boldsymbol {\Psi }}\left({\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}_{f}{\mbox{,}}{\boldsymbol {\mbox{v}}}_{fs}\left({\boldsymbol {\mbox{p}}}\right){\mbox{,}}{\boldsymbol {\mbox{p}}}\right)=$${\boldsymbol {0}}\rightarrow {\boldsymbol {\Psi }}\left({\boldsymbol {\mbox{p}}}\right)=$$0$

( 38)

Nuestra propuesta para la solución de este problema es considerar la conservación de masa como un problema no lineal en p que se pretende resolver mediante un método de Newton aproximado

$-{\frac {\partial {\boldsymbol {\Psi }}}{\partial {\boldsymbol {\mbox{p}}}}}d{\boldsymbol {\mbox{p}}}=$${\boldsymbol {\Psi }}\left({\overset {\mbox{ˆ}}{\boldsymbol {\mbox{v}}}}{\mbox{,}}{\boldsymbol {\mbox{p}}}\right)$

( 39)

A pesar de que existen muchos métodos cuasi-Newton que permiten acelerar la convergencia, se intentará encontrar un método que se mantenga tan cerca como sea posible de las iteraciones de Newton-Raphson estándar. El punto de partida será por lo tanto la ecuación (37) que, se recuerda, se obtuvo asumiendo un comportamiento linealizado de la estructura alrededor de los valores de velocidades y presiones obtenidos al final del bucle 1.

$-{\frac {\partial {\boldsymbol {\Psi }}}{\partial {\boldsymbol {\mbox{p}}}_{n+1}}}=$$\left(\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {\mbox{A}}}_{11}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {\mbox{A}}}_{22}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)+{\frac {\tau }{\rho }}{\boldsymbol {\mbox{L}}}\right)$

( 40)

El uso de dicha matriz garantiza en teoría las mejores velocidades de convergencia (hasta garantizar convergencia en una iteración si la estructura es lineal). El punto clave, sin embargo, es que se pretende seguir asumiendo un comportamiento no lineal de la estructura, uno de los objetivos fijados. Desafortunadamente, la matriz tangente exacta no se conoce explícitamente y sería difícil de construir en un contexto desacoplado, aun sin tener en cuenta su alto coste de construcción. Será necesario por lo tanto emplear una aproximación del tipo

${\boldsymbol {\mbox{A}}}_{aprox}\approx -{\frac {\partial {\boldsymbol {\Psi }}}{\partial {\boldsymbol {\mbox{p}}}_{n+1}}}$

( 41)

Vamos a empezar considerando el caso no estabilizado que se obtiene fijando τ  = 0 en la ecuación (34) .

Cuando una formulación exclusivamente de fluido se utiliza en el bucle 2, el esquema de resolución no tiene ninguna información asociada a la presencia de la estructura. Esto se refleja en la matriz de iteración

${\boldsymbol {\mbox{A}}}_{aprox}=\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\frac {dt}{c_{a}{\rho }_{f}}}{\boldsymbol {\mbox{M}}}_{f}^{-1}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\frac {dt}{c_{a}{\rho }_{f}}}{\left({\boldsymbol {\mbox{M}}}_{fs}^{f}\right)}^{-1}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)$

( 42)

obtenida sencillamente omitiendo el término ${\textstyle {\frac {c_{a}}{dt}}{\rho }_{s}{\boldsymbol {\mbox{M}}}_{fs}^{s}+}$$c_{v}dt{\boldsymbol {\mbox{K}}}$ en el cálculo de A₂₂ en la interfase. Hay que destacar que en una implementación en un único dominio el término ${\textstyle {\boldsymbol {\mbox{M}}}_{fs}^{f}}$debe valer cero para asegurar la aplicación exacta de las condiciones de contorno de velocidad. Remarcablemente, este comportamiento es reproducido satisfactoriamente por la ecuación (40) , puesto que A₂₂  → 0 para el caso que estructuras infinitamente rígidas. También es importante destacar que en el contorno de acoplamiento el término no desaparece puesto que las velocidades no son prescritas sino que se obtienen como resultado del análisis. Introduciendo la notación abreviada ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{L}}}}}$ para indicar la matriz laplaciana discreta,

${\overset {\mbox{ˆ}}{\boldsymbol {\mbox{L}}}}\quad :=$$\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {\mbox{M}}}_{f}^{-1}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\left({\boldsymbol {\mbox{M}}}_{fs}^{f}\right)}^{-1}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)$

( 43)

Por supuesto, el uso de cualquier operador tangente aproximado va a permitir encontrar la solución exacta en convergencia . Sin embargo, la convergencia no está garantizada a no ser que el operador sea una aproximación suficientemente buena al operador lineal original. La aproximación de la ecuación (43) parece ser suficiente para algunas familias de problemas (con estructuras rígidas o pesadas) pero no en otros casos.

Se debe observar también que:

• Para estructuras flexibles y pasos de tiempo grandes el término ${\textstyle c_{v}dt{\boldsymbol {\mbox{K}}}}$ de la ecuación (36) no es negligible.
• Para estructuras ligeras ${\textstyle {\boldsymbol {\mbox{M}}}_{s}^{-1}/dt}$ no es negligible.

Llegados a este punto, encontrar una solución a nuestro problema es relativamente fácil. Se requiere simplemente introducir otra vez en el problema los términos que se han omitido en el acoplamiento Dirichlet-Neumann. Afortunadamente, bajo la hipótesis de pasos de tiempo pequeños esto se puede hacer de forma relativamente sencilla. La idea es que el operador tangente de la ecuación (40) se puede aproximar como

$\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {\mbox{A}}}_{11}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {\mbox{A}}}_{22}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)\approx {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{L}}}}+$$\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {0}}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\left({\rho }_{s}{\boldsymbol {\mbox{M}}}_{fs}^{s}+c_{v}dt{\boldsymbol {\mbox{K}}}\right)}^{-1}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)$

( 44)

donde se ha hecho uso de la definición en la ecuación (43) . El punto interesante es que se puede construir una aproximación mejor al operador tangente con

$\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {0}}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\left({\rho }_{s}{\boldsymbol {\mbox{M}}}_{fs}^{s}+c_{v}dt{\boldsymbol {\mbox{K}}}\right)}^{-1}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)$

( 45)

Se observa que el término adicional:

• Actúa únicamente en la interfase de acoplamiento
• Es mayor cuando la estructura es ligera (en aquellos casos donde la convergencia del esquema Dirichlet-Neumann tiende a fallar)
• Tiende a cero cuando la rigidez de la estructura es dominante (cuando Dirichlet-Neumann converge satisfactoriamente)

Desafortunadamente, el término ${\textstyle {\left({\rho }_{s}{\boldsymbol {\mbox{M}}}_{fs}^{s}+c_{v}dt{\boldsymbol {\mbox{K}}}\right)}^{-1}}$ no se puede usar en la práctica, puesto que requiere la inversión de la matriz de rigidez de la estructura. Además, incluso si se deseara invertir la matriz, asumiendo que el dominio está dividido, la información perteneciente al dominio estructural no está disponible al resolver el fluido a no ser que se proporcione explícitamente.

Nuestros esfuerzos se van a centrar en encontrar una forma de modelar la contribución de la estructura sobre el dominio fluido sin tener que evaluar explícitamente la inversa de la matriz de rigidez estructural. Esto va a permitir definir exactamente la mínima cantidad de información que se requiere transferir de la estructura al fluido.

Una observación atractiva, de acuerdo con el espíritu de las técnicas de paso fraccionado, es que para pasos de tiempo pequeños ${\textstyle {\left({\rho }_{s}{\boldsymbol {\mbox{M}}}_{fs}^{s}+c_{v}dt{\boldsymbol {\mbox{K}}}\right)}^{-1}\rightarrow {\left({\rho }_{s}{\boldsymbol {\mbox{M}}}_{fs}^{s}\right)}^{-1}}$ . La contribución de la interfase a la matriz tangente puede por lo tanto ser modelada como

${\boldsymbol {\mbox{L}}}_{fs}\quad :=\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {0}}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\left({\rho }_{s}{\boldsymbol {\mbox{M}}}_{fs}^{s}\right)}^{-1}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)$

( 46)

Centrándonos en las ecuaciones (46) y (42) podemos ver que el término L_fs difiere únicamente de la matriz laplaciana estándar en los términos de la matriz diagonal en el centro. Esto justifica denominar L_fs como laplaciano de interfase . Teniendo en cuenta que los bloques de la matriz central se aproximan como diagonales, se puede concluir que para cada término de la diagonal con índice II existe un parámetro escalar β_II tal que

${\beta }_{II}{\left({\frac {{\rho }_{f}c_{a}}{dt}}{\boldsymbol {\mbox{M}}}_{fs}^{f}\right)}_{II}^{-1}=$${\left({\frac {{\rho }_{s}c_{a}}{dt}}{\boldsymbol {\mbox{M}}}_{fs}^{s}\right)}_{II}^{-1}$

( 47)

o, una vez invertido,

${\beta }_{II}={\frac {{\frac {{\rho }_{f}c_{a}^{f}}{dt}}{\left({\boldsymbol {\mbox{M}}}_{fs}^{f}\right)}_{II}}{\left({\frac {{\rho }_{s}c_{a}^{s}}{dt}}{\left({\boldsymbol {\mbox{M}}}_{fs}^{s}\right)}_{II}\right)}}$

( 48)

Teniendo en cuenta que la masa es proporcional al tamaño de elemento h_I asociado al nodo I   (${\textstyle M_{II}\propto h_{I}^{3}}$ ) y bajo la hipótesis de estar utilizando el mismo esquema de integración temporal en ambos lados se obtiene la estimación

${\beta }_{II}\propto {\frac {{\rho }_{f}{\left(h_{f}^{3}\right)}_{I}}{{\rho }_{s}{\left(h_{s}^{3}\right)}_{I}}}$

( 49)

Lo cual permite definir una matriz diagonal β tal que los términos sean nulos fuera de la interfase y dados por la ecuación (49) en el contorno de acoplamiento. Una vez obtenida dicha aproximación, la matriz laplaciana de la interfase se puede construir sin necesidad de más comunicación usando la matiz ${\textstyle {\boldsymbol {\mbox{M}}}_{fs}^{f}}$ , que es conocida por el código de resolución del fluido. Esto proporciona la ecuación

${\boldsymbol {\mbox{L}}}_{fs}\quad :=\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {0}}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {\beta }}{\left({\rho }_{f}{\boldsymbol {\mbox{M}}}_{fs}^{f}\right)}^{-1}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)$

( 50)

Esta observación nos proporciona una ventaja crucial, puesto que es conceptualmente más simple dar un valor a β que estimar directamente la matriz laplaciana de interfase o construir la matriz laplaciana acoplada . Se debe remarcar también que únicamente se requiere una estimación de β , puesto que el resultado correcto se obtiene en convergencia.

Una ventaja importante de el enfoque propuesto es la facilidad de estimar el tamaño ${\textstyle h_{i}^{s}}$ en cada nodo de interfase. Dicha cantidad geométrica se puede transferir al dominio fluido, proporcionando la posibilidad de tener en cuenta fácilmente mallas no concordantes.

4.1. Un enfoque mejorado

Para completar la discusión queremos discutir la hipótesis de paso de tiempo pequeño para la estructura. En este sentido es importante observar que la inercia va a dominar sobre los otros términos cuando el paso de tiempo se aproxima al que se usaría para un esquema explícito. Puesto que sería interesante poder usar pasos de tiempo mayores, se debería tener en cuenta de algún modo el efecto de la rigidez (local) de la estructura cuando se evalúa β . Se propone partir de la observación heurística de que dicha rigidez está relacionada con la velocidad c a la cual una onda generada en la interfase se propagará hacia el interior del dominio estructural. El módulo de elasticidad asociado a una onda que se propaga con una velocidad c se puede escribir como k  = ρc² . Por lo tanto, se puede estimar que la rigidez se va a comportar como

${\boldsymbol {\mbox{K}}}_{diag}\propto {\frac {{\rho }_{s}c^{2}}{h_{s}}}$

( 51)

Esta estimación es válida para cualquier sólido, puesto que hasta el momento no se ha asociado la velocidad de la onda a las características mecánicas (o geométricas) de la estructura. Esto significa que la expresión anterior puede dar un criterio para modelar la rigidez en estructuras sólidas y en otros tipos de estructura, como láminas o membranas.

Para poder utilizar esta observación debemos retroceder a la ecuación (45) y desarrollar una aproximación distinta. La idea es aproximar el término que evalúa dicha ecuación como la suma de dos términos en la interfase

${\boldsymbol {\mbox{A}}}_{aprox}=\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {0}}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\left({\rho }_{s}{\boldsymbol {\mbox{M}}}_{fs}^{s}+c_{v}dt{\boldsymbol {\mbox{K}}}\right)}^{-1}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)\approx {\boldsymbol {\mbox{L}}}_{fs}+$${\boldsymbol {\mbox{L}}}_{k}$

( 52)

de los cuales uno es el laplaciano de interfase y el otro es un nuevo término que describe la estructura, es decir

${\boldsymbol {\mbox{L}}}_{k}\quad :=\left({\begin{array}{cc}{\boldsymbol {\mbox{D}}}_{f}&{\boldsymbol {\mbox{D}}}_{fs}\end{array}}\right)\left({\begin{array}{cc}{\boldsymbol {0}}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\left(c_{v}dt{\boldsymbol {\mbox{K}}}\right)}^{-1}\end{array}}\right)\left({\begin{array}{c}{\boldsymbol {\mbox{G}}}_{f}\\{\boldsymbol {\mbox{G}}}_{fs}\end{array}}\right)$

( 53)

Evidentemente el nuevo término no se puede usar todavía a no ser que se haga una simplificación adicional para poder evaluarlo fácilmente. Se destaca que este nuevo término, de forma similar a la matriz laplaciana de interfase L_fs , existe exclusivamente en la interfase de acoplamiento.

En lo que sigue se va a tener en cuenta que una variación de presión Δp aplicada al contorno se transforma en una fuerza normal actuando en el mismo contorno por la acción del operador gradiente G . Ahora es fácil de demostrar que una variación de la fuerza en un nodo I en la interfase de acoplamiento introducida por una variación de presión en el mismo nodo se puede escribir directamente como una función del área de influencia del nodo (definida como la porción de la interfase alrededor del nodo) mediante la relación

$\Delta {\boldsymbol {\mbox{f}}}_{I}={\boldsymbol {\mbox{G}}}_{IJ}\Delta p_{J}=$$A_{I}\Delta p_{I}{\boldsymbol {\mbox{n}}}_{I}\approx h_{s}^{2}\Delta p_{I}{\boldsymbol {\mbox{n}}}_{I}$

( 54)

Usando la aproximación de la ecuación (51) se obtiene

$\Delta {\boldsymbol {\mbox{v}}}={\left(c_{v}dt{\boldsymbol {\mbox{K}}}\right)}^{-1}{\boldsymbol {\mbox{G}}}_{fs}\Delta {\boldsymbol {\mbox{p}}}\approx {\boldsymbol {\mbox{n}}}{\frac {h_{s}}{\rho c^{2}dt}}\Delta p$

( 55)

Aplicando el teorema de la divergencia (estamos asumiendo que la única variación de la velocidad sucede en el contorno, y en la dirección de la normal del mismo contorno) podemos escribir

${\boldsymbol {\mbox{D}}}_{fs}\Delta {\boldsymbol {\mbox{v}}}\approx A\left({\boldsymbol {\mbox{n}}}\cdot {\boldsymbol {\Delta }}{\boldsymbol {\mbox{v}}}\right){\boldsymbol {\mbox{n}}}$

( 56)

Combinando finalmente las ecuaciones (55) y (56) k, y considerando que cada término de la matriz de masa diagonal cercano a la interfase se puede escribir como ${\textstyle M_{II}=A_{I}{h_{f}}_{I}}$ , podemos proponer la siguiente aproximación mejorada para L_k :

${\boldsymbol {\mbox{L}}}_{k}\approx {\frac {1}{{\frac {h_{f}}{h_{s}}}{\rho }_{s}c^{2}dt}}{\boldsymbol {\mbox{M}}}_{f}$

( 57)

que, una vez más, se puede evaluar fácilmente una vez se dispone de una estimación del tamaño de los elementos y de la rigidez de la estructura. Se observa que el segundo término propuesto para la interfase de la matriz es diagonal y remarcablemente similar al término que aparecería (en todo el dominio) si se modelase el fluido como compresible.

El último paso es proporcionar una fórmula que relacione la velocidad del fluido con los parámetros del material. Para un medio elástico semi-infinito, caracterizado por un módulo de Young E y un coeficiente de Poisson ν la velocidad de propagación de una onda de presión se puede tomar como

$c=\left({\frac {E}{{\rho }_{s}}}\right){\frac {\left(1-\nu \right)}{\left(1+\nu \right)\left(1-2\nu \right)}}$

( 58)

Se pueden desarrollar fácilmente formulaciones alternativas para estructuras delgadas (vigas o láminas), para las cuales existen fórmulas analíticas que relacionan la variación de presión con deformaciones de flexión en la dirección normal.

Por ejemplo, para estructuras 2D de tipo viga con un canto t , la rigidez k se puede estimar considerando una viga simplemente apoyada de longitud h_s . Esto resulta en la aproximación

${\rho }_{s}c_{viga}^{2}=h_{s}k=32E{\frac {t^{3}}{h_{s}^{3}}}$

( 59)

Obviamente, esta fórmula es válida mientras el tamaño del elemento h_s sea superior al canto de la viga t .

La aproximación correspondiente para una lámina 3D se puede obtener de una forma parecida mediante la solución analítica de una lámina simplemente apoyada de espesor t y dimensión h_s .

Con este último punto se han aportado todos los ingredientes necesarios para estimar la matriz tangente modificada para la ecuación de presiones. Uno de los ejemplos en la última sección se destinará a discutir la efectividad de las distintas aproximaciones.

5. Comentarios y observaciones

Antes de proceder a los ejemplos se considera necesario destacar algunos detalles interesantes.

La formulación que se propone se basa en la definición de una matriz tangente aproximada para la ecuación de presión. El algoritmo relaja la condición de incompresibilidad de forma consistente en la interfase FSI, lo cual permite tener en cuenta la deformabilidad de la estructura. La formulación también es consistente en el sentido de que, en convergencia la condición de incompresibilidad se impone correctamente.

Se puede establecer un paralelo entre este enfoque y otras técnicas de acoplamiento basadas en las condiciones de contorno de Robin [21] . Simplificando al extremo, tales enfoques se basan en introducir una cierta masa y amortiguamiento en la interfase FSI. Lo que se quiere remarcar es que el efecto de tal acoplamiento es modificar la ecuación de conservación del momento para tener en cuenta la deformación de la estructura. Además, esta modificación se suele introducir en un esquema de resolución monolítico (con un único dominio).

La técnica propuesta realiza una operación conceptualmente similar pero introduce el efecto de la deformabilidad de la estructura en la ecuación de conservación de masa . El desarrollo que se hace en este documento resulta en un operador que se añade a los términos de conservación de masa. Incluso teniendo en cuenta que el término en cuestión se ha obtenido en el contexto de un esquema de segregación de presión, existe la posibilidad de introducir la misma modificación en un esquema monolítico. Creemos, aunque no se ha verificado todavía, que las ventajas derivadas de la modificación de la contribución de la interfase tal como se presenta aquí se pueden obtener también en el enfoque monolítico.

El segundo punto que se quiere remarcar es que en este trabajo no se han utilizado técnicas de aceleración para mejorar la convergencia del acoplamiento fluido-estructura. Existe una extensa literatura que trata el uso de métodos de Newton-Krylov y técnicas de Aitken para mejorar la convergencia de los algoritmos FSI, siendo los primeros preferidos por la comunidad matemática y las segundas comúnmente utilizadas en problemas de ingeniería. Nada impide utilizar este tipo de técnicas en nuestro bucle de acoplamiento y creemos que la aceleración sería efectiva (aunque, dada la buena convergencia del método sin acelerar, no tenemos claro si los resultados de Newton-Krylov serían mejores que los de Aitken). En cualquier caso la discusión de las ventajas de las distintas técnicas de aceleración queda lejos del objetivo de este trabajo, que es demostrar que simplemente con la modificación que se propone para la ecuación de presión se obtiene un método suficientemente bueno como para enfrentarse a problemas fuertemente acoplados.

6. Ejemplos de aplicación

A continuación se incluye una serie de ejemplos de validación con los que se pretende demostrar que el método no solamente permite encontrar una solución, sino que además esta representa una buena aproximación a la solución exacta. Recientemente se han publicado múltiples trabajos con el objetivo de proporcionar soluciones de referencia para estudios futuros. Algunos ejemplos sencillos se presentan en [30] para verificar la capacidad de un esquema de resolución para reproducir correctamente las características de flujos simples. Otro benchmark se está desarrollando actualmente en una colaboración entre distintos grupos en Alemania. Sin embargo, los resultados finales no se han dado a conocer al público en el momento de escribir este documento, lo que los convierte en inadecuados para nuestro propósito. Los ejemplos en [30] , por otro lado, se centran básicamente en en poner a prueba la capacidad del método de solución para reproducir correctamente la transferencia de datos (equilibrio de fuerzas y continuidad de desplazamientos) en la interfaz FSI. El objetivo del presente trabajo, por otro lado, es presentar una nueva estrategia de acoplamiento, bajo la hipótesis de que la transferencia de datos se hace correctamente. Para minimizar otras fuentes de error, para los ejemplos se utilizan mallas concordantes en el contorno de acoplamiento, sin que ello sea un requisito para la aplicación del método en general. En lo que sigue vamos a empezar con un ejemplo muy simple que nos permita evaluar el impacto de la elección de una aproximación a la matriz tangente de presión. A continuación se aplicará el esquema propuesto a la simulación del efecto del flujo sobre un obstáculo flexible y el hinchado de un globo. Finalmente, se validará el código con un ejemplo más aplicado, comparándolo con resultados experimentales para un problema en el campo de flujo sanguíneo.

6.1. Probando la matriz de interfase

La estrategia de solución propuesta se basa en la aproximación del operador de presión acoplado, que relaciona incrementos de presión con incrementos de velocidad en un sistema lineal, por uno más simple obtenido a partir de modelar el efecto de la estructura sobre el dominio fluido. El concepto es reemplazar el problema lineal por otro no lineal para el cual se puede hacer una buena estimación del operador tangente.

Para discutir la efectividad de las distintas opciones que se plantean para aproximar el operador tangente se ha preparado una problema geométricamente simple, pero que requiere una técnica fuertemente acoplada para encontrar la solución. Las principales características del problema se proporcionan en la figura 2 . Una representación de la configuración deformada del problema en la fases descendiente y ascendiente de la solución (puesto que la estructura presenta un movimiento oscilatorio) se puede ver en la figura 3 .

Figura 2. Modelo de referencia: estructura delgada a la izquierda – estructura gruesa a la derecha.

Figura 3. Distintas fases del movimiento: descendiente a la izquierda – ascendente a la derecha.

El ejemplo fue resuelto utilizando una formulación con laplaciano discreto en el dominio fluido. El parámetro de estabilización τ ha sido tomado como τ  = dt en todos los ejemplos para asegurar que no juega un papel en mejorar la convergencia de las pruebas.

Las densidades de la estructura y el fluido se han tomado como ρ_s  = ρ_f  = 1.000 kg/m³ , el módulo de Young se ha fijado en E  = 200e 6 N/m² y la viscosidad en ν  = 0, 01. Se han considerado dos cantos distintos para la estructura: primero un espesor t  = 0, 01m , que modela una viga flexible, y posteriormente un canto t  = 0, 1, que asegura un comportamiento mucho más rígido.

La teoría desarrollada predice que la estructura delgada va a suponer mayores dificultades para obtener convergencia debido a la mayor flexibilidad de la estructura y su menor masa. Asimismo, se espera que la inclusión del término laplaciano en la interfase va a ser suficiente para conseguir convergencia, mientras que la versión mejorada con L_fs  + L_k va a facilitar aún más la convergencia. Además, se espera que la incorporación de la fórmula de la ecuación (59) para estimar la rigidez estructural va a suponer una mejora adicional. Finalmente, se espera que la convergencia va a ser mejor para la estructura más rígida y que, en este caso, la diferencia entre los métodos se va a reducir por la menor influencia de la flexibilidad de la estructura.

Una última observación es que se espera que el número total de iteraciones se va a reducir cuando disminuya el paso de tiempo.

Para verificar esta afirmación se han probado sistemáticamente los dos casos (estructura gruesa y delgada) utilizando tres opciones distintas para la evaluación de la tangente. El número total de iteraciones requeridas para obtener la convergencia de 1.000 pasos de tiempo se dan en los Cuadros 1 y 2 .

Tabla 1. Caso de la estructura delgada: iteraciones totales para 1.000 pasos de tiempo.

dt	Lfs	${\textstyle \left({\boldsymbol {\mbox{L}}}_{fs}+{\boldsymbol {\mbox{L}}}_{k}\right)}$	${\textstyle \left({\boldsymbol {\mbox{L}}}_{fs}+{\boldsymbol {\mbox{L}}}_{k}^{viga}\right)}$
0,0005	6.528	5.967	3.861
0,001	6.727	6.109	4.503
0,002	7.456	7.435	6.318

Tabla 2. Caso de la estructura gruesa: iteraciones totales para 1.000 pasos de tiempo.

dt	Lfs	${\textstyle \left({\boldsymbol {\mbox{L}}}_{fs}+{\boldsymbol {\mbox{L}}}_{k}\right)}$	${\textstyle \left({\boldsymbol {\mbox{L}}}_{fs}+{\boldsymbol {\mbox{L}}}_{k}^{viga}\right)}$
0,0005	3.812	3.557	3.497
0,001	3.413	3.378	3.376
0,002	6.903	6.876	6.854

Los resultados obtenidos confirman satisfactoriamente las predicciones teóricas. El uso del término laplaciano L_fs es suficiente para la convergencia. La velocidad de convergencia aumenta asimismo al incluir el término L_k . También es interesante observar como, de acuerdo con lo esperado, una mejor estimación de la rigidez de la estructura reduce el número de iteraciones, una característica que puede ser importante en el análisis de estructuras de lámina.

6.2. Validación con un problema estándar

A continuación se utiliza un problema benchmark bastante extendido, propuesto originalmente en [8] y reproducido en múltiples trabajos de investigación como [14]  and [31] . En la experiencia de los autores, este problema es un desafío para los esquemas que resuleven la velocidad y la presión de forma alternada. El dominio del modelo se muestra en la figura 4 .

Figura 4. Imagen del dominio de cálculo.

La estructura es una ménsula con las siguientes propiedades:

• Módulo de Young: 2, 3 × 10⁶  N/m²
• Coeficiente de Poisson: 0, 45
• Densidad: 1500,0 kg/m³

El fluido es newtoniano con:

• Viscosidad: 0,145 kg/m s
• Densidad: 956,0 kg/m³

La velocidad de entrada, con sentido de izquierda a derecha, tiene una distribución parabólica que empieza en cero en el fondo del canal y alcanza un máximo en la parte superior. La velocidad horizontal se incrementa según la ley

$v={\frac {0{\mbox{,}}06067}{2}}\left(1-cos\left(1-{\frac {\pi t}{10}}\right)\right)\left[m/s\right]$

( 60)

para t  < 10s   , tomando el valor ${\textstyle v=0{\mbox{,}}06067}$  [m/s] en instantes posteriores.

La discretización de la estructura se ha obtenido usando un único elemento estructural a lo largo del espesor de la ménsula. Esto es posible gracias al uso del elemento triangular descrito en [32] extendido al rango geométricamente no lineal mediante una formulación corrotacional [33] . El elemento es exacto para la flexión de viga, lo cual garantiza una correcta representación de la deformación de la viga. El número total de elementos en el fluido está alrededor de 15, 000, con un tamaño de malla de 0,005 m cerca de la estructura (fig. 5 ).

Figura 5. Solución para el tiempo t  = 30 s.

Los resultados para el desplazamiento del extremo de la viga se muestran en la figura 6 .

Figura 6. Desplazamiento del extremo de la viga.

El análisis fue realizado utilizando un pasos de tiempo de 0, 1 s y 0, 05 s, hasta un tiempo total de 30 s. El tiempo de cálculo total ha sido en todo caso inferior a 30 minutos en un Intel core 2 duo de 2,0 GHz. El número de iteraciones requeridas para cada paso de tiempo se resume en la figura 7 . Para este ejemplo el término laplaciano es suficiente para garantizar una buena convergencia y la importancia del término L_k parece ser negligible. Por otro lado, el total de iteraciones mejora drásticamente al reducir el paso de tiempo, como se esperaba en base a las predicciones teóricas.

Figura 7. Recuento de iteraciones.

Finalmente se quiere destacar que el uso del elemento triangular elegido para la discretización de la estructura conduce a un sistema lineal rígido. Este hecho no representa un problema para la formulación aquí presentada, puesto que se utiliza un algoritmo directo para resolver la estructura, mientras que para resolver la ecuación de Poisson de presión se prefiere un método iterativo.

6.3. Hinchado de un globo

Una prueba muy interesante para verificar los métodos acoplados para la resolución de problemas FSI es el hinchado de un globo estanco. El ejemplo elegido aquí se ha extraído de [16] , que a su vez se inspira en los trabajos de [34] y [35] . La idea es muy simple: dado un flujo de entrada, la masa de fluido almacenada en el volumen delimitado por el globo está determinada para todo instante de tiempo. Así, se puede dar una medida de la validez de los resultados obtenidos comparando la masa obtenida por el modelo computacional con el valor analítico en cada instante. Antes de describir en detalle los resultados se quiere remarcar que la velocidad es fijada en todos los puntos del contorno del fluido. Esto representa un verdadero desafío, puesto que el problema no va a quedar bien definido a no ser que la flexibilidad de la estructura se tenga en cuenta. Esto implica que los métodos desacoplados más comunes, para los cuales durante la solución del fluido no se dispone de información referente al comportamiento de la estructura, no van a poder conseguir un resultado.

En el enfoque desacoplado que se propone, se proporciona una estimación de la flexibilidad de la estructura al código de resolución del fluido (que en particular se usa para resolver las presiones) a través de los términos de interfase. Esto es suficiente para convertir la matriz de presión en invertible, lo cual permite obtener un resultado.

En el ejemplo que se va a resolver se han utilizado los siguientes parámetros, extraídos de [16] :

• Diámetro del globo ϕ  = 2 m
• Espesor de la estructura circundante t  = 2 mm
• Módulo de Young: E  = 1.000 N/m²
• Coeficiente de Poisson: 0, 4
• Densidad de la estructura ρ_s  = 100, 0 kg/m³
• Densidad del fluido ρ_f  = 1, 0 kg/m³
• Viscosidad del fluido μ_f  = 0, 00001 kg/m s

La velocidad de entrada es variable con una dependencia temporal del tipo ${\textstyle (v\left(t\right)=1{\mbox{,}}0+sin\left(\pi t\right)}$ .

Los resultados de la simulación se muestran en la figura 8 .

Figura 8. Hinchado del globo en distintos instantes de tiempo.

En lugar de representar la variación total de volumen, es conveniente representar la variación de volumen en el tiempo (comparada con su valor analítico) puesto que esta cantidad es más sensible a errores de acoplamiento. Los resultados obtenidos para este valor se muestran en la figura 9 .

Figura 9. Derivada de la variación de volumen.

El total de iteraciones requeridas se muestra en la figura 10 .

Figura 10. Iteraciones requeridas hasta convergencia.

Se observa que el algoritmo descrito proporciona una buena solución cuando el bucle interno converge hasta una tolerancia suficiente. Sin embargo, es interesante ver que si la iteración no converge totalmente, o si se relaja la tolerancia, también se obtiene una solución, que es parecida a la que se obtendría si la interfase fuera menos rígida. Es de esperar que, a pesar de imponer un límite en el número de iteraciones, de forma que no se alcance la tolerancia deseada, el algoritmo descrito convergerá a una solución estable aproximada. Para verificar este fenómeno el número máximo de iteraciones se ha limitado a 5. Como se esperaba, se obtiene una buena solución, como se puede ver en la figura 11 .

Figura 11. Derivada de la variación de volumen cuando se limita el número máximo de iteraciones a 5.

Este hecho representa una propiedad muy importante del algoritmo, puesto que el problema se está modificando en cada paso en otro más fácil de resolver, mientras que la solución del problema original se obtiene en convergencia. Esto implica que la estabilidad del problema es prioritaria, siendo la exactitud de la solución únicamente requerida una vez ha convergido, lo cual está en acuerdo con el espíritu de los métodos basados en la introducción de una compresibilidad artificial en la interfase de acoplamiento.

6.4. Validación experimental con una válvula aórtica

En esta sección se presenta un caso de validación experimental. La geometría se muestra en la figura 12 . Consiste en un canal con una membrana delgada que puede pivotar alrededor de su extremo, permitiendo su rotación. En su parte central, el canal se expande en una cavidad mayor, modelada aquí como un semicírculo, que forma el seno aórtico. Se pueden consultar los detalles de la configuración del experimento y los datos asociados en [36] . La membrana tiene un espesor de 1, 5⁻⁴  m y una densidad de 1.000 kg/m³ . Para el fluido, se ha considerado una densidad de 1.000 kg/m³ y una viscosidad constante de 0,0043 kg/m s. El fluido entra en el dominio por el extremo izquierdo de la figura 12 con una velocidad dada por la ley:

$V\left(t\right)=-27{\mbox{,}}61+65{\mbox{,}}25t-48{\mbox{,}}68t^{2}+$$11{\mbox{,}}04t^{3}+0{\mbox{,}}25t^{4}$

( 61)

donde t indica tiempo. Dichos valores se han obtenido por interpolación de resultados experimentales. En el panel inferior de la figura 13 se se representa gráficamente la velocidad de entrada. En el extremo derecho del modelo se impone una presión p  = 0. En el resto del dominio se impone una condición de no deslizamiento.

Figura 12. Geometría del modelo de validación de la válvula aórtica. Dimensiones en metros.

Figura 13. Perfiles interpolados de velocidad comparados con perfiles experimentales en cinco instantes de tiempo, nombrados de acuerdo con el panel inferior derecho. La línea discontinua representa el experimento y la continua los resultados numéricos.

Los resultados para el perfil del flujo se representan en la figura 13 . Se muestra el perfil instantáneo de velocidad para cinco instantes de tiempo, comparado con los resultados experimentales. La imagen del extremo inferior derecho de la figura representa el flujo de entrada en función del tiempo. Se puede ver que, en general, los resultados experimentales concuerdan con los numéricos, excepto en algunas zonas del seno aórtico. Las causas más razonables de este fenómeno estarían relacionadas con la aproximación del flujo de entrada y el modelo usado para pivotar la membrana. Finalmente, la figura 14 muestra el ángulo de apertura de la membrana en función del tiempo. Una vez más, se observa una buena concordancia entre los resultados numéricos y los experimentales. La parte plana de la curva que se observa para los resultados experimentales entre los instantes 0, 8 s y 0, 9 s corresponde a la inversión del flujo e indica el contacto entre la membrana y las paredes del canal. Este contacto es sobrestimado por el modelo numérico.

Figura 14. Ángulo de apertura obtenido (línea continua) comparado con los resultados experimentales (linea discontinua).

7. Conclusiones

En este documento se ha descrito un procedimiento simple que permite resolver problemas FSI en el marco de los esquemas de segregación de presión. El método presentado es computacionalmente asequible y por lo tanto apto para la solución eficiente de problemas de gran tamaño. La idea básica es tener en cuenta la relación entre la presión y la deformación de la estructura como un problema no lineal. Dicho problema no lineal se resuelve mediante el uso de un método de Newton aproximado. A continuación se deriva una forma explícita para el operador tangente aproximado, que se demuestra que toma la forma de un término adicional en la interfase que se debe añadir a la ecuación laplaciana de presión. Se demuestra que la adición del término propuesto no requiere el uso de métodos iterativos para resolver el sistema planteado en cada iteración de Newton. Esto supone una ventaja importante cuando se quiere trabajar con láminas o sistemas estructurales no lineales, para los cuales son más adecuados los métodos directos de solución.

Agradecimientos

Este trabajo ha sido parcialmente respaldado por los proyectos XPRES y Consolider Sedurec, del Programa Nacional de I+D del Ministerio de Educación y Ciencia español.

References

1. [1] C. Felippa, K.C. Park, C. Farhat; Partitioned analysis of coupled mechanical systems; Comput. Methods Appl. Mech. Eng., 190 (2001), pp. 3247–3270
2. [2] C. Farhat, P. Geuzaine, G. Brown; Application of a three-field non-linear fluid-structure formulation to the prediction of the aeroelastic parameters of an F-16 fighter; Comput. Fluids, 32 (2003), pp. 3–29
3. [3] S. Piperno, C. Farhat, B. Larroutorou; Partitioned procedures for the transient solution of coupled aeroelastic problems - part1 - model problem, theory and two dimensional applications; Comput. Methods Appl. Mech. Eng., 124 (1995), pp. 79–112
4. [4] S. Piperno, C. Farhat; Partitioned procedures for the transient solution of coupled aeroelastic problems - part2 - energy transfer analysis and three dimensional applications; Comput. Methods Appl. Mech. Eng., 190 (1995), pp. 3147–3170
5. [5] M. Lesoinne, C. Farhat; Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations; Comput. Methods Appl. Mech. Eng., 134 (1996), pp. 71–90
6. [6] Forster Ch. Robust methods for fluid-structure interaction with stabilised finite elements. PhD thesis, Institut fur Baustatik und Baudynamik, Universitat Stuttgart, 2007.
7. [7] J. Donea, A. Huerta; Finite Element Methods for Flow Problems; J. Wiley and Sons (2002)
8. [8] D.P. Mok, W.A. Wall; Partitioned analysis schemes for the transient interaction of incompressible flows and nonlinear flexible structures; Trends in computational structural mechanics, Barcelona (2001)
9. [9] U. Kuttler, W.A. Wall, Fixed-point fluid-structure interaction solvers with dynamic relaxation, Journal Computational Mechanics.
10. [10] J. Vierendeels, L. Lanoye, J. Degroote, P. Verdonck; Implicit coupling of partitioned fluid-structure interaction problems with reduced order models; Comput. Struct., 85 (2007), pp. 970–976
11. [11] R. Lohner; Applied CFD Techniques: An introduction Based on FEM; John Wiley and Sons (2001)
12. [12] S. Turek; Efficient solvers for incompressible flow problems: An algorithmic approach in view of computational aspects; Springer-Verlag (1998)
13. [13] O.C. Zienkiewicz, R.L. Taylor, P. Nithiarasu; The Finite Element Method for Fluid Dynamics; ELSEVIER Butterworth-Heinemann (2005)
14. [14] R. Rossi, Light-weight structures - Numerical Analysis and Coupling Issues. PhD thesis, University of Padova, Italy, 2005.
15. [15] E. Oñate, R. Rossi. Analysis of some partitioned algorithms for fluid-structure interaction. Engineering Computations, accepted for publication, 2009.
16. [16] Y. Bazilevs, V.M. Calo, T.J.R. Hughes, Y. Zhang. Isogeometric fluid-structure interaction: Theory, algorithms, and computations. Technical report, ICES, 2008.
17. [17] M.A. Fernandez, J.-F. Gerbeau, C. Grandmont; A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid; Int. J. Numerical Methods Eng., 69 (2007), pp. 794–821
18. [18] S. Badia, A. Quaini, A. Quarteroni; Splitting methods based on algebraic factorization for fluid-structure interaction; SIAM J. Sci. Comput., 30 (2008), pp. 1778–1805
19. [19] P. Causin, J-F. Gerbeau, F. Nobile; Added-mass effect in the design of partitioned algorithms for fluid-structure problems; Comput. Methods Appl. Mech. Eng., 194 (2005), pp. 4506–4527
20. [20] Ch. Forster, W.A. Wall, E. Ramm; Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows; Comput. Methods Appl. Mech. Eng., 196 (2007), pp. 1278–1293
21. [21] S. Badia, F. Nobile, C. Vergara; Fluid-structure partitioned procedures based on robin transmission conditions; J. Comput. Phys., 227 (2008), pp. 7027–7051
22. [22] W. Dettmer, D. Peric; A computational framework for fluid-structure interaction: Finite element formulation and applications; Comput. Methods Appl. Mech. Eng., 195 (2006), pp. 5754–5779
23. [23] P. Dadvand, R. Rossi, E. Oñate; An object-oriented environment for developing finite element codes for multi-disciplinary applications; Arch. Comput. Methods Eng., 17 (2010), pp. 253–297
24. [24] P. Davdand. framework for developing finite element codes for multi-disciplinary applications. 2007.
25. [25] S.R. Idelsohn, J. Marti, A. Limache, E. Oñate; Unified lagrangian formulation for elastic solids and incompressible fluids. application to fluid-structure interaction problems via the pfem; Comput. Methods Appl. Mech. Eng., 197 (2008), pp. 1762–1776
26. [26] E. Oñate, S.R. Idelsohn, M.A. Celigueta, R. Rossi; Advances in the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows; Comput. Methods Appl. Mech. Eng., 197 (2008), pp. 1777–1800
27. [27] R. Codina; Stabilized finite element approximation of transient incompressible flows using orthogonal subscales; Comput. Methods Appl. Mech. Eng., 191 (2002), pp. 4295–4321
28. [28] R. Codina; A stabilized finite element method for generalized stationary incompressible flows; Comput. Methods Appl. Mech. Eng., 190 (20–21) (2001), pp. 2681–2706
29. [29] R. Codina; Pressure stability in fractional step finite elements methods for incompressible flows; J. Comput. Phys. (2001)
30. [30] K.-J. Bathe, G.A. Ledezma; Benchmark problems for incompressible fluid flows with structural interactions; Comput. Struct., 85 (2007), pp. 628–644
31. [31] G. Valdes, Nonlinear Analysis of Orthotropic Membrane and Shell Structures Including Fluid-Structure Interaction. PhD thesis, CIMNE, 2007.
32. [32] C.A. Felippa; A study of optimal membrane triangles with drilling freedoms; Comput. Methods Appl. Mech. Eng., 192 (2003), pp. 2125–2168
33. [33] P. Khosravi, R. Ganesan, R. Sedaghati; Corotational non-linear analysis of thin plates and shells using a new shell element; IJNME, 69 (2007), pp. 859–885
34. [34] U. Kuttler, C. Forster, W.A. Wall; Efficient approaches for the fluid structure interaction with fully enclosed incompressible flow domains.; ECCOMAS CFD, Egmond aan Zee (2006)
35. [35] T.E. Tezduyar, S. Sathe; Modelling of fluid-structure interactions with the space-time finite elements: Solution techniques; Int. J. Numerical Methods Fluids, 54 (2007), pp. 855–900
36. [36] J.M.A. Stijnen, J. de Hart, P.H.M. Bovendeerd, F.N. Van de Vosse; Evaluation of a ficticious domain mehtod for predicting the dynamic response of mechanical heart valves; J. Fluids Struct., 19 (2004), pp. 835–850