A scalable preconditioning framework for stabilized contact mechanics with hydraulically active fractures

Resumen Técnico: Un Marco de Preacondicionamiento Escalable para la Mecánica de Contacto Estabilizada con Fracturas Hidráulicamente Activas

Planteamiento del Problema
El artículo aborda el desafío computacional de simular la física acoplada de la mecánica de contacto con fricción y el flujo de fluidos dentro de redes de fracturas embebidas en un medio poroso. Este problema acoplado es crítico para aplicaciones de subsuelo como la producción de energía geotérmica, la fracturación hidráulica y el secuestro de CO2. Los procesos físicos están estrechamente acoplados: la apertura y el deslizamiento de la fractura impulsan el flujo de fluidos, mientras que las variaciones en la presión del fluido perturban el estado de esfuerzos e influyen en la mecánica de contacto.

Los autores emplean un Modelo de Fractura Discreta (DFM, por sus siglas en inglés) donde las fracturas se representan explícitamente. La discretización utiliza elementos finitos continuos de bajo orden para el campo de desplazamiento mecánico y un esquema de volúmenes finitos centrado en la celda para la presión del fluido. Las restricciones de contacto se imponen mediante multiplicadores de Lagrange centrados en la celda (que representan fuerzas de tracción) que comparten la misma malla que el campo de presión, eliminando la necesidad de interpolación entre campos. Sin embargo, esta formulación específica no es uniformemente estable inf-sup (inestable LBB) y requiere una estabilización algebraica global.

El sistema lineal resultante, derivado de la linealización Newton-Krylov del problema no lineal, es una matriz jacobiana de bloque $3 \times 3$, grande, dispersa y no simétrica. La estructura de la matriz es dinámica, cambiando su naturaleza (de un punto de silla simétrico a uno no simétrico indefinido) a medida que la red de fracturas evoluciona a través de regiones de adherencia (stick), deslizamiento (slip) y apertura (open) durante la simulación. Los resolvedores directos estándar son computacionalmente prohibitivos para problemas a gran escala, lo que hace necesario el uso de resolvedores iterativos robustos y escalables.

Metodología
La contribución central es el diseño de un marco de preacondicionamiento escalable basado en la partición de bloques informada por la física y técnicas de estado del arte de Multigrid Algebraico (AMG). El marco opera restringiendo el sistema multifísico totalmente acoplado a subproblemas de física única, resolviéndolos de forma aproximada y luego prolongando la solución de vuelta al sistema completo.

Se proponen y analizan dos secuencias de reducción (enfoques) distintas:

1. El enfoque $t-p-u$ (Tracción-Presión-Desplazamiento):

   • Estrategia: Elimina simultáneamente las variables de tracción ($t$) y presión ($p$) para reducir el sistema a una única ecuación de equilibrio para el desplazamiento ($u$).
   • Complemento de Schur: Esto resulta en un complemento de Schur $S = A + C_1 \tilde{H}^{-1}C_2 - Q_1 T^{-1}Q_2$.
   • Aproximación: El bloque de tracción $H$ se aproxima mediante una matriz de bloque diagonal $\tilde{H}$, y el bloque de presión $T$ se aproxima mediante su diagonal $\tilde{T}$. La inversa del complemento de Schur $S$ resultante se aplica de forma inexacta utilizando un operador AMG (específicamente GAMG de PETSc).
   • Características: Este enfoque es computacionalmente más barato por iteración (una llamada a AMG) pero depende de la definición positiva de $S$. El artículo señala que $S$ puede volverse indefinido y no simétrico dependiendo del estado de la fractura (específicamente cuando existen regiones abiertas), lo que podría comprometer la efectividad del resolvedor AMG.

2. El enfoque $t-u-p$ (Tracción-Desplazamiento-Presión):

   • Estrategia: Primero elimina la tracción ($t$) para formar un complemento de Schur de primer nivel $S_1 = A + C_1 \tilde{H}^{-1}C_2$ (una matriz de elasticidad con rigidez ficticia). Luego, se elimina el desplazamiento ($u$) para formar un complemento de Schur de segundo nivel $S_2 = T - Q_2 S_1^{-1} Q_1$.
   • Aproximación: La inversa de $S_1$ se aplica utilizando AMG. El complemento de Schur de segundo nivel $S_2$ se aproxima mediante una matriz de esfuerzo fijo de tipo diagonal $S_K$, que puede calcularse algebraicamente o mediante estimaciones de volumen local.
   • Características: Este enfoque requiere dos llamadas a AMG por cada aplicación del precondicionador (mayor costo por iteración). Sin embargo, ofrece mayor robustez porque el complemento de Schur de primer nivel $S_1$ está garantizado teóricamente como definido positivo, independientemente del estado de la fractura. La posible indefinitud se confina al complemento de Schur de segundo nivel, el cual es simétrico y puede manejarse eficazmente mediante resolvedores directos anidados.

Contribuciones Clave

• Diseño del Marco: Los autores presentan un marco de preacondicionamiento unificado capaz de manejar la naturaleza no simétrica, indefinida y variable en el tiempo de la jacobiana de bloque $3 \times 3$ proveniente de la mecánica de contacto estabilizada con fracturas hidráulicamente activas.
• Análisis Comparativo: El artículo proporciona una comparación teórica y numérica rigurosa de las secuencias de reducción $t-p-u$ y $t-u-p$. Establece que, si bien el enfoque $t-p-u$ es teóricamente eficiente, el enfoque $t-u-p$ es superior en robustez para simulaciones complejas del mundo real donde los estados de las fracturas son impredecibles.
• Demostración de Escalabilidad: El trabajo demuestra que los precondicionadores propuestos logran una escalabilidad algorítmica, manteniendo un número casi constante de iteraciones de Krylov a medida que aumenta el tamaño del problema y a medida que cambia la relación entre las incógnitas de la fractura y las de la matriz (bulk).

Resultados Numéricos
El marco fue probado en tres conjuntos de experimentos:

1. Robustez (Prueba 1): Un problema de fractura única reveló que, mientras el enfoque $t-p-u$ funciona bien con resolvedores directos, puede fallar en la convergencia cuando se usa AMG si el complemento de Schur se vuelve indefinido. El enfoque $t-u-p$ mantuvo la convergencia en todas las configuraciones probadas.
2. Escalabilidad Débil (Prueba 2): Las pruebas en mallas refinadas uniformemente (hasta ~760,000 incógnitas) y variando las proporciones de variables 2D-a-3D mostraron que el conteo de iteraciones del resolvedor lineal se mantuvo estable (promediando ~17–24 iteraciones de GMRES) independientemente del refinamiento de la malla o del tamaño relativo de la red de fracturas.
3. Aplicación del Mundo Real (Prueba 3): Una simulación de un pozo inclinado que intersecta nueve fracturas en un escenario de fracturación hidráulica (más de 1 millón de incógnitas) confirmó la eficiencia computacional. A pesar de que el resolvedor no lineal requirió muchos pasos de conjunto activo (active-set) e iteraciones de Newton debido a las complejas transiciones de adherencia/deslizamiento/apertura, el resolvedor lineal se mantuvo estable y eficiente, con un promedio de 24.2 iteraciones de GMRES por paso.

Significancia y Reivindicaciones
El artículo afirma que el marco de preacondicionamiento propuesto aborda con éxito el cuello de botella computacional de simular el acoplamiento de contacto con fricción y flujo de fluidos en medios fracturados. Al aprovechar la partición de bloques basada en la física y las técnicas de multigrid, el método logra la escalabilidad necesaria para simulaciones paralelas a gran escala.

Los autores concluyen modestamente que, aunque el enfoque $t-p-u$ ofrece una ligera ventaja en eficiencia teórica y menor costo de aplicación, el enfoque $t-u-p$ es la elección preferida para aplicaciones prácticas debido a su superior robustez contra la indefinitud de las matrices de complemento de Schur que surgen durante la compleja evolución de las fracturas. El trabajo establece una base para resolver problemas geomecánicos cada vez más grandes y complejos, con trabajos futuros dirigidos hacia implementaciones de alto rendimiento y extensiones al flujo multifásico.