Discretización por elementos finitos del problema termo-difusivo-mecánico con grandes deformaciones_prueba_latex

Discretización por elementos finitos del problema termo-difusivo-mecánico con grandes deformaciones

Ángel Ortiz-Toranzo and Ignacio Romero ${}^{1,2,}$

${}^{1}$ Departamento de Ingeniería Mecánica, Universidad Politécnica de Madrid, Madrid, 28006, Spain

${}^{2}$ IMDEA Materials, Madrid, 28906, Spain

Resumen

En algunos problemas relevantes de la ingeniería, el campo mecánico está fuertemente acoplado con el de temperatura y el de transporte de masa. La solución de estos problemas es compleja sobre todo cuando las deformaciones son grandes. Por ello, a menudo se acude a aproximaciones numéricas que pueden encontrar todos los campos involucrados en el problema acoplado, sus efectos y sus interacciones. En este artículo se describe una discretización por el método de los elementos finitos del problema acoplado de difusión, temperatura y deformación incluyendo el rango no lineal de deformaciones. La formulación tiene en cuenta todos los acoplamientos posibles entre los tres campos de estudio y se proporcionan todos los detalles necesarios para su completa implementación llenando un hueco de la literatura de estos métodos. Asimismo, el artículo describe la fundamentación termodinámica de los problemas acoplados termo-difusivo-mecánicos, insistiendo en la derivación de las ecuaciones de balance y restricciones que siguen de la segunda ley de la termodinámica. Los resultados del artículo serán de interés para investigadores que necesiten implementar las ecuaciones de problemas acoplados en códigos de elementos finitos y, en particular, para aplicaciones en el modelado de baterías, del comportamiento de metales bajo los efectos del hidrógeno, de geles, etc.

In certain relevant engineering problems, the mechanical field is strongly coupled with the temperature and mass transport fields. The solution to these problems is complex, especially when deformations are large. Therefore, numerical approximations are often used to find all the fields involved in the coupled problem, their effects, and their interactions. This article describes a finite element method discretization of the coupled

problem of diffusion, temperature, and deformation, including the nonlinear range of deformations. The formulation takes into account all possible couplings between the three fields of study, and all the necessary details for its complete implementation are provided, filling a gap in the literature of these methods. Additionally, the article describes the thermodynamic foundation of coupled thermo-diffusive-mechanical problems, emphasizing the derivation of balance equations and constraints that follow from the second law of thermodynamics. The results of the article will be of interest to researchers who need to implement the equations of coupled problems in finite element codes, particularly for applications in battery modeling, the behavior of metals under the effects of hydrogen, gels, etc.

(¹)

(²) #1

OPEN ACCESS

Received: 21/07/2024

Accepted: 29/10/2024

DOI

10.23967/j.rimni.2024.10.56364

Keywords:

Problemas acoplados elementos finitos problema termo-difusivo-mecánico grandes deformaciones implementación termo-mecánica tensión-difusión

Keywords:

Coupled problems finite elements thermo-diffusive-mechanical problem large deformations implementation thermo-mechanics stress-diffusion

1 Introducción

Los problemas de acoplamiento de múltiples campos que involucran la transmisión de calor, el transporte de masa y la deformación mecánica son trascendentes en diversas aplicaciones de ingeniería y tecnología avanzada, tales como las relacionadas con materiales compuestos reforzados con fibras [1,2], dispositivos de conversión y almacenamiento de energía [3,4], tejidos biológicos [5], mezclas porosas de geomateriales como el hormigón [6-8], procesos de fragilización del acero por el hidrógeno, relevante en aplicaciones de ingeniería como las relacionadas con el almacenamiento y transporte de petróleo, gas e hidrógeno [9], estudios relacionados con los geles poliméricos [10-12], en el desarrollo de portadores para la administración de fármacos y en el estudio de tejidos [13], y en otras aplicaciones en el ámbito de la biología y medicina [14]. El acoplamiento de los campos térmico, difusivo y mecánico es la característica común de todos estos problemas y, por lo tanto, constituyen el objetivo de multitud de estudios, tanto para poder predecir el comportamiento y posibles problemas en el desempeño de estructuras y componentes de máquinas como para el desarrollo de nuevos materiales.

Durante los últimos años se han realizado multitud de investigaciones centradas en los problemas de acoplamiento entre tensiones mecánicas, transporte de masa y campo térmico. En líneas generales estos problemas se pueden conceptualizar como un sólido deformable poroso que actúa de huésped y que se encuentra embebido en líquidos, mezclas gaseosas o átomos o iones de soluto que pueden migrar dentro o fuera del sólido huésped induciendo dilatación o contracción y modificación del campo de temperatura. Existe una extensa literatura focalizada en los problemas de acoplamiento de dos o más campos, desde los estudios de Gibbs sobre la difusión de fluidos en sólidos, las investigaciones de Biot sobre el asentamiento de suelos bajo consolidación (similar en muchos casos al proceso de exprimir agua en un medio poroso elástico, seguido de sus estudios relativos a problemas termoviscoelásticos acoplados) hasta los estudios recientes sobre geles poliméricos y materiales elastoméricos [11,12,15,16], materiales porosos [17,18], fenómenos difusivos en los procesos de fragilización de aceros por hidrógeno [19,20] donde se presenta una implementación totalmente acoplada del campo mecánico y fenómenos de difusión y tejidos biológicos [5,21].

Mientras tanto, promovido por el desarrollo de baterías y celdas de combustible, el problema de la difusión en sólidos acoplado con el comportamiento térmico y mecánico ha sido ampliamente estudiado [22-24]. Con la llegada de los vehículos de tecnología eléctrica, las baterías de ion litio se han convertido en el estándar tecnológico en la motorización de este tipo de vehículos [25,26],

además de ser los principales sistemas de almacenamiento para equipos eléctricos y electrónicos de tamaño reducido [27], implantes biomédicos y componentes microelectrónicos autoalimentados [28], aplicaciones aeroespaciales [29]. Además, las baterías de ion litio, tienen aplicaciones en redes de distribución, nivelación de carga y energía renovables [30], pudiendo considerarse una tecnología con un papel fundamental en el camino hacia un futuro energético más limpio. Por todo ello, desarrollar baterías con alta capacidad energética y larga vida útil es un desafío tecnológico trascendente, siendo el daño producido en los electrodos por el proceso de litiación el cuello de botella en el desarrollo de baterías de alta densidad y uno de los retos tecnológicos presentes.

El objetivo de este estudio es presentar de manera completa el modelo y la discretización por elementos finitos del problema acoplado de tres campos, mecánico, térmico y transporte de masa, en el rango no lineal de deformación, teniendo en cuenta todos los acoplamientos posibles entre los campos pero limitando la respuesta mecánica al comportamiento elástico. Un paso previo para alcanzar el objetivo propuesto es la definición completa de los fundamentos termodinámicos de dichos problemas no lineales. Para ello, el artículo presenta de manera resumida, pero completa, estos fundamentos y los emplea para formular métodos numéricos que resuelven los tres campos, manteniendo en todo momento una conexión de la discretización y la formulación continua de estos problemas fuertemente acoplados.

El presente artículo está organizado de la siguiente forma: en la sección 2 se describirán las ecuaciones que rigen el comportamiento acoplado de los cuerpos deformables elásticos bajo los efectos de la temperatura y el transporte de masa, considerando un modelo de difusión de una sola especie. El modelo tendrá en cuenta los acoplamientos bidireccionales existentes entre los tres campos: mecánico, térmico y difusión de masa. En dicha sección, se presenta un conjunto de relaciones constitutivas para la energía libre de Helmholtz, el flujo másico de acuerdo con los postulados del principio de disipación y el flujo de calor. En la sección 3, se presenta la discretización espacio-temporal de las ecuaciones diferenciales parciales resultantes y una vía para su resolución numérica basada en la formulación de elementos finitos tipo Galerkin. En la sección 4, se presentarán algunos ejemplos de aplicación del modelo planteado. Para ilustrar la generalidad del modelo se resolverán dos problemas completamente diferentes: por un lado un ejemplo de geles y por otro se mostrará una simulación relacionada con el comportamiento de las baterías de ion litio; ambos ejemplos utilizarán modelos donde los campos termo-difusivo-mecánicos están fuertemente acoplados. Por último, se finalizará con unas conclusiones en la sección 5.

2 Planteamiento del problema

En esta sección se presenta el problema mecánico acoplado con el fenómeno de transporte de masa y con el de conducción de calor. El propósito, por tanto, es desarrollar un modelo matemático, termodinámicamente consistente, para el medio continuo con grandes deformaciones que permita describir la interacción termo-difusivo-mecánica para sólidos capaces de absorber especies químicas. Posteriormente, la solución del modelo se aproximará mediante un modelo numérico cuya implementación se describirá en detalle.

En lo que sigue, se trabajará con un cuerpo que se considerará medio continuo, definido matemáticamente como un conjunto de partículas ${\textstyle {\mathcal {B}}=\left\{P_{1},P_{2},P_{3}\ldots \right\}}$, de tal manera que existe un conjunto de aplicaciones inyectivas ${\textstyle {\mathcal {K}}=\{\chi \}}$ que transforman el conjunto ${\textstyle {\mathcal {B}}}$ en conjuntos abiertos del espacio ${\textstyle \mathbb {R} ^{3}}$. Esto es, para toda partícula ${\textstyle P\in {\mathcal {B}}}$

${\textstyle {\boldsymbol {x}}=\chi (P)\in \mathbb {R} ^{3};\quad \chi ({\mathcal {B}})={\mathcal {O}}\subseteq \mathbb {R} ^{3}}$, abierto, ${\textstyle y\chi ^{-1}({\mathcal {O}})={\mathcal {B}}}$,

cumpliéndose que para cualquier pareja de aplicaciones ${\textstyle \chi _{1},\chi _{2}\in {\mathcal {K}}}$ su composición

${\textstyle \varphi _{12}:\chi _{1}({\mathcal {B}})\subseteq \mathbb {R} ^{3}\rightarrow \chi _{2}({\mathcal {B}})\subseteq \mathbb {R} ^{3},\quad \varphi _{12}=\chi _{2}\circ \chi _{1}^{-1}}$,

es diferenciable. Se denomina configuración a cada una de las infinitas aplicaciones ${\textstyle \chi }$. Como es habitual, se tomará como configuración de referencia, ${\textstyle \chi _{ref}}$, la configuración sin deformar. Con frecuencia se utiliza la expresión "configuración" para referirse al conjunto imagen de una de ellas, es decir, ${\textstyle \chi ({\mathcal {B}})}$; su uso es tan extendido que aquí también se utilizará.

Las configuraciones de referencia y una cualquiera deformada definen una aplicación entre ${\textstyle {\mathcal {B}}_{ref}}$ y ${\textstyle {\mathcal {B}}_{def}}$ que se denomina deformación (véase figura 1):

$\varphi :{\mathcal {B}}_{ref}\rightarrow {\mathcal {B}}_{def},\quad \varphi :=\chi _{def}\circ \chi _{ref}^{-1}.$

Se define un movimiento como una familia de configuraciones ${\textstyle \varphi _{t}}$ que depende de un parámetro ${\textstyle t}$, el tiempo, de forma diferenciable. Esto es, si ${\textstyle {\mathcal {K}}({\mathcal {B}})}$ es el espacio de configuraciones posibles de un cuerpo continuo ${\textstyle {\mathcal {B}}}$, entonces ${\textstyle \varphi _{t}}$ es una curva en ${\textstyle {\mathcal {K}}({\mathcal {B}})}$.

Figura 1: Campo de deformación de un medio continuo

El gradiente de deformación se define como ${\textstyle {\boldsymbol {F}}=D{\boldsymbol {\varphi }}}$ y se adoptará la hipótesis de que admite una descomposición multiplicativa en tres partes, ${\textstyle {\boldsymbol {F}}^{m},{\boldsymbol {F}}^{s}}$ y ${\textstyle {\boldsymbol {F}}^{t}}$, tal que:

${\textstyle {\boldsymbol {F}}={\boldsymbol {F}}^{m}{\boldsymbol {F}}^{s}{\boldsymbol {F}}^{t}\quad }$ con ${\textstyle \quad {\boldsymbol {F}}^{s}=\lambda ^{s}\mathbf {1} ;{\boldsymbol {F}}^{t}=\lambda ^{t}\mathbf {1} ;\quad \lambda ^{s}>0\quad \lambda ^{t}>0}$.

En esta expresión ${\textstyle {\boldsymbol {F}}^{s}}$ está asociado a la deformación del material debida a las moléculas absorbidas por la estructura base, ${\textstyle \lambda ^{s}}$ representa el estiramiento correspondiente, ${\textstyle {\boldsymbol {F}}^{t}}$ está relacionado con la deformación del material debido al efecto térmico, ${\textstyle \lambda ^{t}}$ representa su correspondiente estiramiento y ${\textstyle {\boldsymbol {F}}^{m}}$ está asociado a la deformación mecánica.

Esta descomposición multiplicativa del tensor gradiente de deformación se le atribuye a Flory [31] y su uso está muy extendido [32,33]. Existen diversas situaciones en las que también se utiliza

una descomposición multiplicativa del tensor gradiente de deformación, tales como termoelasticidad, elastoplasticidad y biomecánica.

El jacobiano de la transformación se define como ${\textstyle J=\operatorname {det} ({\boldsymbol {F}})}$, continuando con la notación indicada, se tiene que ${\textstyle J^{m}=\operatorname {det} \left({\boldsymbol {F}}^{m}\right),J^{s}=\operatorname {det} \left({\boldsymbol {F}}^{s}\right)}$ y ${\textstyle J^{t}=\operatorname {det} \left({\boldsymbol {F}}^{t}\right)}$. Consecuentemente, a la vista de la descomposición multiplicativa del gradiente de deformación planteada en (4), el jacobiano de la transformación puede escribirse:

${\textstyle J=J^{m}J^{s}J^{t}\quad }$ con ${\textstyle \quad J^{s}=\left(\lambda ^{s}\right)^{3},J^{t}=\left(\lambda ^{t}\right)^{3}}$.

Si se supone que el volumen de un mol de moléculas de soluto, que se denotará como ${\textstyle \Omega }$, es constante y que además todo el cambio por unidad de volumen de referencia se debe al cambio de concentración de las moléculas de soluto en la estructura base, entonces:

${\textstyle J^{s}=1+\Omega g\quad \Rightarrow \quad {\dot {J}}^{s}=\Omega {\dot {g}}}$,

siendo ${\textstyle g}$ la concentración del soluto por unidad de volumen de referencia de la estructural huésped en seco. A la vista de la relación entre ${\textstyle J^{s}}$ y ${\textstyle \lambda ^{s}}$ planteada en (5), la primera ecuación de (6) puede expresarse también como

${\textstyle \lambda ^{s}=(1+\Omega g)^{\frac {1}{3}}}$.

A continuación se acometerá la descripción matemática de las leyes fundamentales de la física que rigen el comportamiento de un medio continuo: las ecuaciones de balance.

2.1 Forma local de las ecuaciones de balance

La ecuación de balance de masa, que determina la evolución de la concentración de soluto en la estructura base, junto con la ecuación de la cantidad de movimiento y con la ecuación de balance de energía establecen el sistema local de ecuaciones diferenciales en derivadas parciales que, conjuntamente con la restricción impuesta por el segundo principio de la termodinámica, gobiernan el problema acoplado termo-difusivo-mecánico. Estos principios son válidos siempre y para cualquier material.

Para un cuerpo arbitrario con configuración de referencia ${\textstyle {\mathcal {B}}_{ref}}$ y contorno ${\textstyle \Gamma }$, la forma local (material) de las tres ecuaciones de balance puede escribirse como [34,35]:

donde ${\textstyle {\boldsymbol {P}}}$ representa el primer tensor de Piola Kirchhoff, ${\textstyle {\boldsymbol {B}}}$ la descripción material del campo vectorial de fuerzas másicas por unidad de masa, ${\textstyle \rho _{0}}$ la densidad, ${\textstyle {\boldsymbol {V}}}$ la velocidad material ${\textstyle {\boldsymbol {V}}={\frac {d}{dt}}{\boldsymbol {\varphi }},{\dot {\boldsymbol {V}}}={\frac {d}{dt}}{\boldsymbol {V}}}$ es la aceleración material, ${\textstyle E}$ es la energía interna por unidad de volumen, ${\textstyle {\boldsymbol {Q}}}$ el vector flujo de calor y ${\textstyle Q}$ es el calor generado por las fuentes internas por unidad de volumen en la configuración de referencia y por unidad de tiempo, ${\textstyle {\boldsymbol {h}}}$ es el flujo másico, ${\textstyle \mu }$ el potencial químico y ${\textstyle g}$ la concentración. Los símbolos GRAD ( ${\textstyle \cdot }$ ) y ${\textstyle \operatorname {DIV} (\cdot )}$ denotan, respectivamente, el gradiente y la divergencia materiales.

Las ecuaciones correspondientes al balance de la cantidad de movimiento, balance de masa y balance de la energía están complementadas por el segundo principio de la termodinámica cuya

expresión material es:

${\textstyle {\dot {s}}\geq {\frac {Q}{\Theta }}-\operatorname {DIV} \left({\frac {\boldsymbol {Q}}{\Theta }}\right)}$.

Multiplicando ambos términos de la ecuación (9) por ${\textstyle \Theta }$ y definiendo ${\textstyle {\boldsymbol {G}}=-{\frac {1}{\Theta }}\operatorname {GRAD} (\Theta )}$, se llega a la expresión

${\textstyle {\boldsymbol {Q}}\cdot {\boldsymbol {G}}-{\boldsymbol {h}}\cdot \operatorname {GRAD} (\mu )\geq 0}$.

Las ecuaciones (8) definen un problema inicial de valores de contorno con las siguientes condiciones de frontera:

$\left.{\begin{array}{rrrrrr}{\boldsymbol {\varphi }}={\overline {\boldsymbol {\varphi }}}&{en}\Gamma _{\varphi },&\Theta ={\bar {\Theta }}&{en}\Gamma _{\Theta },&{\begin{array}{r}\mu \end{array}}={\bar {\mu }}&{en}\Gamma _{\mu },\\{\boldsymbol {PN}}={\overline {\boldsymbol {t}}}={\overline {\boldsymbol {t}}}&{en}\Gamma _{t},&{\boldsymbol {QN}}={\overline {\boldsymbol {q}}}_{n}&{en}\Gamma _{q},&{\boldsymbol {h}}{\boldsymbol {N}}={\overline {\boldsymbol {h}}}_{n}&{en}\Gamma _{h},\end{array}}\right\}\quad \mathrm {t} \in [0,\mathrm {~T} ],$

y con las siguientes condiciones iniciales:

${\textstyle {\boldsymbol {\varphi }}={\boldsymbol {\varphi }}_{0},\quad {\boldsymbol {V}}={\boldsymbol {V}}_{0},\quad \Theta =\Theta _{0},\quad \mu =\mu _{0},\quad }$ para ${\textstyle t=0}$,

siendo ${\textstyle {\boldsymbol {N}}}$ la normal unitaria hacia el exterior de ${\textstyle \Gamma }$. Se ha utilizado la notación ${\textstyle ({\bar {\cdot }})}$ y ${\textstyle (\cdot )_{0}}$ para valores conocidos de las funciones en sus respectivos dominios y el contorno ${\textstyle \Gamma }$ se ha descompuesto como:

${\textstyle \Gamma ={\overline {\Gamma _{\varphi }\cup \Gamma _{t}}}={\overline {\Gamma _{\Theta }\cup \Gamma _{q}}}={\overline {\Gamma _{\mu }\cup \Gamma _{h}}}}$

con ${\textstyle \quad \Gamma _{\varphi }\cap \Gamma _{t}=\Gamma _{\Theta }\cap \Gamma _{q}=\Gamma _{\mu }\cap \Gamma _{h}=\emptyset }$.

Nótese que la variable fundamental del problema de transporte de masa es el potencial electroquímico ${\textstyle \mu }$ y no la concentración ${\textstyle g}$, que es su variable conjugada. Esto está motivado porque físicamente no se puede imponer la concentración en la frontera; tan sólo el potencial electroquímico. Sin embargo, sí que se puede imponer el flujo de masa en el contorno y no existen dispositivos para imponer el flujo de potencial electroquímico. Ver, por ejemplo [36].

En virtud del postulado II de la termodinámica se considerará la existencia de una función, denominada entropía ${\textstyle s}$, de los parámetros extensivos, definida para todos los estados de equilibrio y que tiene la propiedad siguiente: los valores que toman los parámetros extensivos, en ausencia de ligaduras internas, son aquellos que maximizan la entropía respecto al conjunto de los estado de equilibrio ligados. Además se supondrá que el estado termodinámico estará unívocamente determinado por el tensor gradiente de deformación ${\textstyle {\boldsymbol {F}}}$, la entropía ${\textstyle s}$ y la concentración ${\textstyle g}$; la energía interna, que es una función de estado termodinámico, podrá escribirse como ${\textstyle E=E({\boldsymbol {F}},s,g)}$. Mediante la transformada de Legendre pueden obtenerse otros potenciales, definiéndose la energía libre ${\textstyle \Psi }$ como:

${\textstyle \Psi ({\boldsymbol {F}},\Theta ,g)=E({\boldsymbol {F}},s,g)-s\Theta }$,

y el potencial gran canónico ${\textstyle {\mathfrak {G}}}$ :

${\textstyle {\mathfrak {G}}({\boldsymbol {F}},\Theta ,\mu )=\Psi ({\boldsymbol {F}},\Theta ,g)-g\mu }$.

Teniendo en cuenta que la entropía ${\textstyle s}$ puede definirse como:

${\textstyle s=-{\frac {\partial {\mathfrak {G}}}{\partial \Theta }}}$,

aplicando la regla de la cadena en el primer término de la segunda ecuación de (8) se sigue que:

$\underbrace {\frac {\partial E}{\partial F}} _{\boldsymbol {P}}{\dot {\boldsymbol {F}}}+\underbrace {\frac {\partial E}{\partial s}} _{\Theta }{\dot {s}}+\underbrace {\frac {\partial E}{\partial g}} _{\mu }{\dot {g}}={\boldsymbol {P}}:{\dot {\boldsymbol {F}}}+Q-\operatorname {DIV} ({\boldsymbol {Q}})+\mu {\dot {g}}-{\boldsymbol {h}}\operatorname {GRAD} (\mu ),$

o equivalentemente:

${\textstyle \Theta {\dot {s}}=Q-\operatorname {DIV} ({\boldsymbol {Q}})-{\boldsymbol {h}}\operatorname {GRAD} (\mu )}$.

Previamente a establecer la formulación débil del problema, resulta necesario caracterizar dos tipos de funciones. En primer lugar se tienen las funciones de prueba, que son aquellas que pertenecen al espacio de funciones de cuadrado integrable cuya derivada también es de cuadrado integrable (espacio de Hilbert ${\textstyle H^{1}}$ ) y cumplen la primera de las condiciones de contorno establecidas en (11). Así se tiene que los espacios de funciones de prueba para los tres campos del problema son:

${\textstyle {\mathcal {S}}_{1}=\left\{\varphi \mid {\boldsymbol {\varphi }}\in \left[H^{1}({\mathcal {B}})\right]^{3},{\boldsymbol {\varphi }}={\bar {\varphi }}\right.}$ en ${\textstyle \left.\Gamma _{\varphi }\right\}}$,

${\textstyle {\mathcal {S}}_{2}=\left\{\Theta \mid \Theta \in H^{1}({\mathcal {B}}),\Theta ={\bar {\Theta }}\right.}$ en ${\textstyle \left.\Gamma _{\Theta }\right\}}$,

${\textstyle {\mathcal {S}}_{3}=\left\{\mu \mid \mu \in H^{1}({\mathcal {B}}),\mu ={\bar {\mu }}\right.}$ en ${\textstyle \left.\Gamma _{\mu }\right\}}$,

y las funciones de ponderación:

${\begin{aligned}&{\mathcal {W}}_{1}=\left\{\delta \varphi \mid \delta \varphi \in \left[H^{1}({\mathcal {B}})\right]^{3},\delta \varphi =0{en}\Gamma _{\varphi }\right\},\\&{\mathcal {W}}_{2}=\left\{\delta \Theta \mid \delta \Theta \in H^{1}({\mathcal {B}}),\delta \Theta =0{en}\Gamma _{\Theta }\right\},\\&{\mathcal {W}}_{3}=\left\{\delta \mu \mid \delta \mu \in H^{1}({\mathcal {B}}),\delta \mu =0{en}\Gamma _{\mu }\right\}.\end{aligned}}$

Denominando ${\textstyle {\boldsymbol {S}}}$ el segundo tensor de Piola-Kirchhoff, la formulación débil correspondiente a la ecuación del balance de la cantidad de movimiento, primera ecuación de (8), puede escribirse como:

$\left\langle {\boldsymbol {S}},{\boldsymbol {F}}^{T}\operatorname {GRAD} (\delta {\boldsymbol {\varphi }})\right\rangle +\left\langle \rho _{0}{\dot {\boldsymbol {V}}},\delta {\boldsymbol {\varphi }}\right\rangle =\left\langle \rho _{0}{\boldsymbol {B}},\delta {\boldsymbol {\varphi }}\right\rangle +\langle {\boldsymbol {T}},\delta {\boldsymbol {\varphi }}\rangle _{\Gamma _{t}}.$

En esta ecuación, y a partir de ahora, se ha utilizado la notación ${\textstyle \langle \cdot ,\cdot \rangle _{\Xi }}$ para indicar el producto ${\textstyle L_{2}}$ de dos campos vectoriales o tensoriales sobre un conjunto ${\textstyle \Xi }$. Cuando no se indica ningún conjunto se sobreentenderá que el producto escalar se extiende a toda la configuración de referencia. Más adelante se considerará el caso en el que las fuerzas inerciales son despreciables frente al resto. En esta situación la forma débil de la ecuación del balance de cantidad de movimiento se simplifica a:

$\left\langle {\boldsymbol {S}},{\boldsymbol {F}}^{T}\operatorname {GRAD} (\delta {\boldsymbol {\varphi }})\right\rangle =\left\langle \rho _{0}{\boldsymbol {B}},\delta {\boldsymbol {\varphi }}\right\rangle +\langle {\boldsymbol {T}},\delta {\boldsymbol {\varphi }}\rangle _{\Gamma _{t}}.$

Respecto a la ecuación correspondiente al balance de la energía, segunda ecuación de (8), multiplicando dicha ecuación por una función de ponderación ${\textstyle \delta \Theta \in {\mathcal {W}}_{2}}$, integrando y aplicando el teorema de Gauss, se tiene que:

${\textstyle \langle \Theta {\dot {s}},\delta \Theta \rangle -\langle {\boldsymbol {Q}},\operatorname {GRAD} (\delta \Theta )\rangle +\langle {\boldsymbol {h}}\cdot \operatorname {GRAD} (\mu ),\delta \Theta \rangle =\langle Q,\delta \Theta \rangle -\langle {\boldsymbol {Q}}\cdot {\boldsymbol {N}},\delta \Theta \rangle _{\Gamma _{q}}}$,

que es la formulación débil del balance de la energía en coordenadas materiales. En cuanto al balance de masa, tercera ecuación de (8), multiplicando dicha ecuación por una función de ponderación ${\textstyle \delta \mu \in }$ ${\textstyle {\mathcal {W}}_{3}}$ e integrando por partes se llega a:

$\langle {\dot {g}},\delta \mu \rangle =\langle -\operatorname {DIV} ({\boldsymbol {h}}),\delta \mu \rangle$

Á. Ortiz-Toranzo and I. Romero,

Utilizando la relación ${\textstyle \operatorname {DIV} ({\boldsymbol {h}}\delta \mu )=\operatorname {DIV} ({\boldsymbol {h}})\delta \mu +{\boldsymbol {h}}\operatorname {GRAD} (\delta \mu )}$ y aplicando el teorema de la divergencia, se obtiene:

$\langle {\dot {g}},\delta \mu \rangle -\langle {\boldsymbol {h}},\operatorname {GRAD} (\delta \mu )\rangle =-\langle {\boldsymbol {h}}\cdot {\boldsymbol {N}},\delta \mu \rangle _{\Gamma _{h}}.$

Obteniéndose así la forma débil del balance de masa.

2.2 Ecuaciones constitutivas del material

En el apartado precedente, las leyes de balance junto con la descripción cinemática habitual de los medios continuos [35] permiten establecer las ecuaciones del equilibrio del problema termo-difusivomecánico planteado. Sin embargo, es necesario definir una regla de comportamiento del material que permita vincular la cinética del problema con el equilibrio del cuerpo, a fin de poder resolver las ecuaciones de gobierno establecidas. Esta regla de comportamiento del material es lo que se denomina habitualmente modelo constitutivo.

Los modelos constitutivos elásticos más ampliamente estudiados en Mecánica de Medios continuos son los modelos hiperelásticos también llamados modelos elásticos de Green. Estos modelos son aquellos en los que se supone la existencia de una función de energía libre de Helmholtz ${\textstyle \Psi }$ definida por unidad de volumen, de forma que el primer tensor de tensiones de Piola-Kirchhoff se puede calcular simplemente derivando la función de energía libre de Helmholtz. Para el modelado del comportamiento del material termo-difusivo-mecánico, se plantea extender el modelo hiperelástico de forma que también exista un funcional, la energía libre de Helmholtz, cuyas derivadas den lugar a los campos conjugados.

Para facilitar el modelado de los procesos termo-difusivo-mecánicos y sus acoplamientos, supondremos que la energía libre estará construida por una expresión aditiva de una parte que establece la relación de acoplamiento térmico-difusivo ( ${\textstyle \Psi ^{t-d}}$ ) y otra parte correspondiente al acoplamiento termomecánico ( ${\textstyle \Psi ^{t-m}}$ ), otra difusivo-mecánico ( ${\textstyle \Psi ^{d-m}}$ ) además de las contribuciones puramente mecánica, térmica y difusiva ( ${\textstyle \Psi ^{m},\Psi ^{t},\Psi ^{d}}$ ). Según esta hipótesis simplificadora, la energía libre tendría el siguiente aspecto:

${\textstyle \Psi ({\boldsymbol {F}},\Theta ,g)=\Psi ^{m}({\boldsymbol {F}})+\Psi ^{t}(\Theta )+\Psi ^{d}(\mu )+\Psi ^{t-d}(\Theta ,g)+\Psi ^{t-m}({\boldsymbol {F}},\Theta )+\Psi ^{d-m}({\boldsymbol {F}},g)}$.

Considerando en este caso que el material es hiperelástico, la energía libre de Helmholtz se expresaría ${\textstyle \Psi =\Psi ({\boldsymbol {F}}({\boldsymbol {X}},t),\Theta ,g,{\boldsymbol {X}})}$, y el primer tensor de tensiones de Piola-Kirchhoff, se escribiría:

${\textstyle {\boldsymbol {P}}={\frac {\partial \Psi ({\boldsymbol {F}},\Theta ,g,{\boldsymbol {X}})}{\partial {\boldsymbol {F}}}}}$.

La definición de un material hiperelástico dada por la ecuación (27) no tiene en cuenta el principio de objetividad del modelo constitutivo, que establece que las relaciones constitutivas deben de ser válidas para cualquier observador. Denominando ${\textstyle {\boldsymbol {C}}}$ el tensor derecho de Cauchy-Green y ${\textstyle {\boldsymbol {C}}}$ el tensor derecho de Cauchy-Green, una forma alternativa de definición de material hiperelástico, que incluye intrínsecamente el principio de objetividad, es:

${\textstyle {\boldsymbol {S}}=2{\frac {\partial \Psi ({\boldsymbol {C}},\Theta ,g,{\boldsymbol {X}})}{\partial {\boldsymbol {C}}}}}$,

La ecuación constitutiva correspondiente a la conducción de calor establece que el flujo de calor depende linealmente del gradiente térmico:

${\textstyle {\boldsymbol {Q}}=-{\boldsymbol {K}}\operatorname {=} \Theta }$,

donde ${\textstyle {\boldsymbol {K}}}$ es el tensor material de conductividad térmica, que en el caso, como se va a suponer, de un medio isótropo puede escribirse como ${\textstyle K=K\mathbf {1} }$, siendo ${\textstyle K}$ la constante de conductividad térmica.

En cuanto a la ecuación constitutiva correspondiente al flujo difusivo, vendrá dada por la Ley de Fick, que en este contexto establece que el flujo del soluto depende linealmente del gradiente del potencial químico ([34], ec. 66.23):

${\textstyle {\boldsymbol {h}}=-{\boldsymbol {M}}\operatorname {=} \mu }$,

siendo ${\textstyle {\boldsymbol {M}}}$ el tensor de movilidad. Considerando medio isotrópico, en tensor de movilidad puede escribirse ${\textstyle {\boldsymbol {M}}=m{\boldsymbol {I}}}$ siendo ${\textstyle m}$ el coeficiente de movilidad (escalar) y la ecuación anterior puede expresarse de manera más sencilla:

${\textstyle {\boldsymbol {h}}=-m\operatorname {=} \mu }$,

donde se supondrá que el coeficiente de movilidad dependerá de la temperatura según la expresión:

${\textstyle m={\frac {Dg}{R\Theta J}}}$,

siendo ${\textstyle D}$ el coeficiente de difusión, ${\textstyle R}$ es la constante universal de los gases ideales y ${\textstyle \Theta }$ la temperatura absoluta.

El potencial químico podrá calcularse como derivada de primer orden del potencial (26), mediante la expresión:

${\textstyle \mu ={\frac {\partial \Psi }{\partial g}}}$.

En cuanto a la capacidad térmica:

${\textstyle c=-\Theta {\frac {\partial ^{2}\Psi }{\partial \Theta ^{2}}}}$.

Teniendo en cuenta que ${\textstyle \Theta }$ es positiva, la positividad estricta de la capacidad térmica para todo valor de ${\textstyle \Theta }$ es equivalente a expresar que la función entropía es estrictamente creciente con la temperatura ${\textstyle \Theta }$ y , asimismo, equivale a indicar que ${\textstyle {\frac {\partial \Psi ^{2}}{\partial \Theta ^{2}}}<0}$ para todo valor de ${\textstyle \Theta }$, esto es, que ${\textstyle \Psi }$ es estrictamente cóncava considerando la variable ${\textstyle \Theta }$.

En lo que respecta a las derivadas segundas del potencial termodinámico correspondiente a la energía libre de Helmholtz, se define el tensor de acoplamiento termo-mecánico, ${\textstyle {\boldsymbol {M}}}$, y su correspondiente parte simétrica, ${\textstyle {\boldsymbol {M}}^{s}}$ :

${\textstyle {\boldsymbol {M}}=\Theta {\frac {\partial ^{2}\Psi }{\partial {\boldsymbol {F}}\partial \Theta }},\quad {\boldsymbol {M}}^{s}=\Theta {\frac {\partial ^{2}\Psi }{\partial \Theta \partial {\boldsymbol {C}}}}}$.

Para el cálculo de ${\textstyle {\boldsymbol {M}}}$ se tendrá en cuenta que:

${\textstyle {\boldsymbol {M}}^{s}=\Theta {\frac {\partial ^{2}\Psi }{\partial \Theta \partial C}}={\frac {\Theta }{2}}{\frac {\partial {\boldsymbol {S}}}{\partial \Theta }}={\frac {\Theta }{2}}{\frac {\partial \left({\boldsymbol {F}}^{-1}{\boldsymbol {P}}\right)}{\partial \Theta }}={\frac {\Theta }{2}}{\boldsymbol {F}}^{-1}{\frac {\partial ^{2}\Psi }{\partial \Theta \partial {\boldsymbol {F}}}}={\frac {1}{2}}{\boldsymbol {F}}^{-1}{\boldsymbol {M}}}$,

consecuentemente

${\textstyle {\boldsymbol {M}}=2{\boldsymbol {F}}{\boldsymbol {M}}^{s}}$.

Se definirá un segundo tensor de acoplamiento, difusivo-mecánico, que vendrá dado por:

${\textstyle {\boldsymbol {L}}={\frac {\partial ^{2}\Psi }{\partial {\boldsymbol {F}}\partial \mu }}={\frac {\partial {\boldsymbol {P}}}{\partial \mu }}}$.

En este caso, al existir tres campos acoplados, surge un nuevo acoplamiento, en este caso un escalar: el coeficiente de acoplamiento difusivo-térmico:

$T=\Theta {\frac {\partial ^{2}\Psi }{\partial \Theta \partial g}}.$

3 Discretización de las ecuaciones y linealización

A continuación, se presentará la formulación del problema algebraico discreto asociado a las ecuaciones que gobiernan el problema de sólidos con grandes deformaciones y acoplamiento termo-difusivo-mecánico y que resulta de la discretización espacial y temporal de dichas ecuaciones.

Para ello, se partirá de la formulación débil del problema expresado mediante las ecuaciones (22), (23) y (25). En primer lugar se dividirá el dominio en el que está definido el problema, ${\textstyle \Upsilon }$, en ${\textstyle n}$ elementos o subdominios ${\textstyle \Upsilon ^{(h)}}$ y se reemplazan los espacios de dimensión infinita ${\textstyle {\mathcal {S}}_{1},{\mathcal {S}}_{2}}$ y ${\textstyle {\mathcal {S}}_{3}}$ a los cuales pertenecen las soluciones ${\textstyle \varphi ,\Theta }$ y ${\textstyle \mu }$, por los espacios de dimensión finita ${\textstyle {\mathcal {S}}_{1}^{h},{\mathcal {S}}_{2}^{h}}$ y ${\textstyle {\mathcal {S}}_{3}^{h}}$, tales que ${\textstyle {\mathcal {S}}_{1}^{h}\subset {\mathcal {S}}_{1},{\mathcal {S}}_{2}^{h}\subset {\mathcal {S}}_{2}}$ y ${\textstyle {\mathcal {S}}_{3}^{h}\subset {\mathcal {S}}_{3}}$ a los cuales pertenecen las funciones aproximadas ${\textstyle \varphi ^{h},\Theta ^{h}}$ y ${\textstyle \mu ^{h}}$, respectivamente. Asimismo, los espacios de dimensión infinita ${\textstyle {\mathcal {W}}_{1},{\mathcal {W}}_{2}}$ y ${\textstyle {\mathcal {W}}_{3}}$, a los cuales pertenecen las funciones de ponderación ${\textstyle \delta \varphi }$, ${\textstyle \delta \Theta }$ y ${\textstyle \delta \mu }$, se reemplazan por los espacios de dimensión finita ${\textstyle {\mathcal {W}}_{1}^{h},{\mathcal {W}}_{2}^{h}}$ y ${\textstyle {\mathcal {W}}_{3}^{h}}$, tales que ${\textstyle {\mathcal {W}}_{1}^{h}\subset {\mathcal {W}}_{1},{\mathcal {W}}_{2}^{h}\subset }$ ${\textstyle {\mathcal {W}}_{2}}$ y ${\textstyle {\mathcal {W}}_{3}^{h}\subset {\mathcal {W}}_{3}}$, a los cuales pertenecen las funciones aproximadas ${\textstyle \delta \varphi ^{h},\delta \Theta ^{h}}$ y ${\textstyle \delta \mu ^{h}}$, respectivamente.

Suponiendo que cada miembro de ${\textstyle {\mathcal {S}}_{1}^{h},{\mathcal {S}}_{2}^{h}}$ y ${\textstyle {\mathcal {S}}_{3}^{h}}$ puede expresarse de la forma:

${\textstyle {\boldsymbol {\varphi }}^{h}={\boldsymbol {\varphi }}_{v}^{h}+{\boldsymbol {\varphi }}_{g}^{h},\quad \Theta ^{h}=\Theta _{v}^{h}+\Theta _{g}^{h},\mu ^{h}=\mu _{v}^{h}+\mu _{g}^{h}}$,

donde ${\textstyle \varphi _{v}^{h}\in {\mathcal {W}}_{1}^{h},\Theta _{v}^{h}\in {\mathcal {W}}_{2}^{h}}$ y ${\textstyle \mu _{v}^{h}\in {\mathcal {W}}_{3}^{h}}$ satisfacen las condiciones de contorno ${\textstyle \varphi _{g}^{h}={\boldsymbol {\varphi }}^{h},\Theta _{g}^{h}=\Theta ^{h}\mathrm {y} }$ ${\textstyle \mu _{g}^{h}=\mu ^{h}}$ en ${\textstyle \Gamma _{\varphi },\Gamma _{\Theta }\mathrm {y} \Gamma _{\mu }}$, respectivamente.

Cualquiera de las funciones de ponderación ${\textstyle \delta \varphi ^{h},\delta \Theta ^{h}}$ y ${\textstyle \delta \mu ^{h}}$, al pertenecer a espacios lineales de dimensión finita, se pueden escribir como combinación de ${\textstyle n}$ funciones de forma, ${\textstyle N_{i}\operatorname {con} i=1,2,\ldots ,n}$, linealmente independientes. La expresión discretizada correspondiente al balance de la cantidad de movimiento se obtiene sustituyendo en su formulación débil ${\textstyle \delta {\boldsymbol {\varphi }}}$ por su aproximación ${\textstyle \delta {\boldsymbol {\varphi }}^{h}}$, expresada como combinación lineal de las funciones de forma ${\textstyle N_{i}}$, quedaría:

$\left\langle {\boldsymbol {S}},{\boldsymbol {F}}^{T}\left(\sum _{a}\nabla N^{a}\delta {\boldsymbol {\varphi }}^{a}\right)\right\rangle =\left\langle \rho _{0}{\boldsymbol {B}},\sum _{a}N^{a}\delta {\boldsymbol {\varphi }}^{a}\right\rangle +\left\langle {\boldsymbol {T}},\sum _{a}N^{a}\delta {\boldsymbol {\varphi }}^{a}\right\rangle _{\Gamma _{t}}.$

Del mismo modo, sustituyendo ${\textstyle \delta \Theta }$ por ${\textstyle \delta \Theta ^{h}=\sum _{a=1}^{N_{nd}}N^{a}({\boldsymbol {X}})\delta \Theta ^{a}}$ en la ecuación correspondiente a la formulación débil del balance de la energía, se obtiene la ecuación discretizada del balance de la energía:

${\textstyle \left\langle \Theta {\dot {s}},\sum _{a}N^{a}\delta \Theta ^{a}\right\rangle -\left\langle {\boldsymbol {Q}},\sum _{a}\nabla N^{a}\delta \Theta ^{a}\right\rangle +\left\langle {\boldsymbol {h}}\cdot \operatorname {GRAD} (\mu ),\sum _{a}N^{a}\delta \Theta ^{a}\right\rangle }$

${\textstyle =\left\langle Q,\sum _{a}N^{a}\delta \Theta ^{a}\right\rangle -\left\langle {\boldsymbol {Q}}\cdot {\boldsymbol {N}},\sum _{a}N^{a}\delta \Theta ^{a}\right\rangle _{\Gamma _{q}}}$.

De forma similar, a partir de la formulación débil del balance de masa, sustituyendo ${\textstyle \delta \mu }$ por su aproximación correspondiente ${\textstyle \delta \mu ^{h}}$ expresada como combinación lineal de las funciones de forma, esto

es, ${\textstyle \delta \mu ^{h}=\sum _{a=1}^{N_{nd}}N^{a}({\boldsymbol {X}})\delta \mu ^{a}}$, se obtiene la formulación discretizada de la expresión correspondiente al balance de masa:

$\left\langle {\dot {g}},\sum _{a}N^{a}\delta \mu ^{a}\right\rangle -\left\langle {\boldsymbol {h}},\sum _{a}\nabla N^{a}\delta \mu ^{a}\right\rangle =-\left\langle {\boldsymbol {h}}\cdot {\boldsymbol {N}},\sum _{a}N^{a}\delta \mu ^{a}\right\rangle _{\Gamma _{h}}.$

Se ha realizado la discretización espacial de las ecuaciones de evolución por aplicación del método de elementos finitos, obteniendo así las ecuaciones (41), (42) y (43) procedentes de la aplicación del balance de la cantidad de movimiento, balance de la energía y principio de conservación de la masa respectivamente. Sin embargo, en las ecuaciones (42) y (43) aparecen las derivadas temporales ${\textstyle {\dot {\Theta }},{\dot {\boldsymbol {F}}}}$ y ${\textstyle {\dot {g}}}$, por lo que dichas ecuaciones tambien deben ser discretizadas con respecto al tiempo. Para ello se realizará una partición del intervalo de integración ${\textstyle [0,T]}$ en ${\textstyle N}$ subintervalos del tipo ${\textstyle \left[t_{n},t_{n+1}\right]}$, con ${\textstyle 0=t_{0}<t_{1}<\ldots t_{N}=T}$ tal que ${\textstyle \Delta t=t_{n+1}-t_{n}}$ y se utilizará el método implícito de Euler, según el cual:

${\textstyle {\boldsymbol {F}}_{n+1}\approx {\boldsymbol {F}}_{n}+\Delta t\cdot {\dot {\boldsymbol {F}}}_{n+1}}$,

${\textstyle s_{n+1}\approx s_{n}+\Delta t\cdot {\dot {s}}_{n+1}}$,

${\textstyle g_{n+1}\approx g_{n}+\Delta t\cdot {\dot {g}}_{n+1}}$,

donde los subíndices ${\textstyle n}$ y ${\textstyle n+1}$ denotan el instante de la evaluación, ya sea ${\textstyle t_{n}\mathrm {o} t_{n+1}}$. Según lo anterior, sustituyendo los valores de ${\textstyle {\dot {\Theta }},{\dot {\boldsymbol {F}}}}$ y ${\textstyle {\dot {g}}}$, de las expresiones (44), en la ecuación (42) quedaría:

${\begin{aligned}&\left\langle \Theta {\frac {s_{n+1}-s_{n}}{\Delta t}},\sum _{a}N^{a}\delta \Theta ^{a}\right\rangle -\left\langle {\boldsymbol {Q}},\sum _{a}\nabla N^{a}\delta \Theta ^{a}\right\rangle +\\&\quad +\left\langle {\boldsymbol {h}}\cdot \operatorname {GRAD} (\mu ),\sum _{a}N^{a}\delta \Theta ^{a}\right\rangle =\left\langle Q,\sum _{a}N^{a}\delta \Theta ^{a}\right\rangle -\left\langle {\boldsymbol {Q}}\cdot {\boldsymbol {N}},\left(\sum _{a}N^{a}\delta \Theta ^{a}\right)\right\rangle _{\Gamma _{q}}.\end{aligned}}$

Sustituyendo el valor de ${\textstyle {\dot {g}}}$, de la última expresión de (44), en la ecuación (43) puede escribirse:

${\textstyle \left\langle {\frac {g_{n+1}-g_{n}}{\Delta t}},\sum _{a}N^{a}\delta \mu ^{a}\right\rangle -\left\langle {\boldsymbol {h}},\sum _{a}\nabla N^{a}\delta \mu ^{a}\right\rangle =-\left\langle {\boldsymbol {h}}\cdot {\boldsymbol {N}},\sum _{a}N^{a}\delta \mu ^{a}\right\rangle _{\Gamma _{h}}}$.

Las ecuaciones de conservación-balance expuestas, constituyen un sistema de ecuaciones no lineal que deben ser previamente linealizadas en aras de acometer su resolución numérica por las técnicas habituales de resolución de sistemas de ecuaciones lineales basadas en el método de Newton Raphson. Los detalles del proceso de linealización pueden encontrarse en el Apéndice A.

4 Aplicaciones

El objetivo de esta sección es ilustrar las posibilidades que ofrece el modelo planteado de cara a resolver diferentes problemas acoplados, destacando además la generalidad del método propuesto que es capaz de resolver problemas no lineales complejos que involucran tres campos, o únicamente uno o dos campos, bloqueando la variación de las variables correspondientes mediante la imposición de las correspondientes condiciones de contorno. En todos los ejemplos se ha empleado una discretización espacial con elementos hexaédricos trilineales con una formulación de deformación supuesta (la formulación "B-barra") para los términos de la energía mecánica.

4.1 Ejemplo: geles elastoméricos con dos campos acoplados

Como primer ejemplo, para ilustrar la confluencia del modelo de tres campos a un modelo de dos campos, se resolverá un problema de tensión-difusión acoplado, imponiendo en cada punto el valor de

la temperatura. Se trata de los problemas resueltos mediante modelos de dos campos tensión-difusión tradicionales [32], variacionales [11], y sobre los que hay estudios experimentales [37].

En particular, se considerará un tubo cilíndrico de gel, con diámetro exterior ${\textstyle D}$, espesor ${\textstyle t}$ y longitud ${\textstyle h}$ al que se le sujeta una de sus bases mientras se inyecta fluido por la cara opuesta mediante la imposición de un potencial químico, utilizando el modelo de tres campos presentado bloqueando las variaciones térmicas. La ecuación del potencia químico utilizada para establecer la condición de contorno propuesta es:

${\textstyle \mu (t)=\mu ^{0}+\mu _{0}\exp \left(-t/t_{d}\right)\quad \operatorname {con} t_{d}=60}$.

Para el modelado del comportamiento del material termo-difusivo-mecánico, se va a plantear un modelo constitutivo hiperelástico, de forma que el problema termo-difusivo-mecánico puede describirse a partir un funcional, energía libre de Helmholtz, dada por la expresión:

${\begin{aligned}\Psi _{R}=&\mu _{0}g+R\Theta g\left(\ln \left({\frac {\Omega g}{1+\Omega g}}\right)+\chi (\Theta )\left({\frac {1}{1+\Omega g}}\right)\right)\\&+{\frac {1}{2}}G(\Theta )\left(3\left({\bar {\lambda }}^{2}-1\right)-2\ln J\right)+J^{s}\left({\frac {1}{2}}K\left(\ln J^{e}\right)^{2}\right)+c_{0}\left(\left(\Theta {-\Theta }_{0}\right)-\Theta \ln \left({\frac {\Theta }{\Theta }}\right)\right)\\&-3\alpha K\left(\Theta {-\Theta }_{0}\right)\ln J-3\beta K\left(g-g_{0}\right)\ln J.\end{aligned}}$

En la expresión anterior ${\textstyle \mu _{0}}$ es el potencial químico de referencia, ${\textstyle g}$ es la concentración de fluido en moles de fluido por unidad de volumen de estructura base seca correspondiente al fenómeno difusivo, ${\textstyle R}$ es la constante universal de los gases ideales, ${\textstyle \theta }$ es la temperatura absoluta, ${\textstyle \Omega }$ denota el volumen de un mol de moléculas de fluido, ${\textstyle \chi }$ es un parámetro adimensional que establece una medida de la entalpía de la mezcla conocido como parámetro de Flory-Huggins, ${\textstyle G(\theta )=N_{R}k_{B}\theta }$ es el módulo de cizalladura, que se considerará dependiente de la temperatura, de la constante de Bolzman ( ${\textstyle k_{B}}$ ) y del número de cadenas de polímero por volumen de referencia, ${\textstyle K}$ es el módulo de rigidez volumétrica de la estructura base, ${\textstyle \alpha }$ y ${\textstyle \beta }$ son los coeficientes de expansión térmica y química respectivamente, y ${\textstyle {\bar {\lambda }}}$ representa el alargamiento efectivo, definido como [34]:

${\textstyle {\bar {\lambda }}{\stackrel {def}{=}}{\sqrt {{\frac {1}{3}}\left(\lambda _{1}^{2}+\lambda _{2}^{2}+\lambda _{3}^{2}\right)}}={\sqrt {{\frac {1}{3}}\operatorname {tr} {\boldsymbol {B}}}}={\sqrt {{\frac {1}{3}}\operatorname {tr} {\boldsymbol {C}}}}}$,

donde ${\textstyle \lambda _{i}i=1,2,3}$, son los estiramientos principales, y se ha tenido en cuenta que ${\textstyle {\boldsymbol {B}}={\boldsymbol {F}}{\boldsymbol {F}}^{T}}$ (tensor izquierdo de Cauchy-Green), y la propiedad conmutativa de la traza del producto.

El material utilizado se trata de un gel y sus propiedades son las mismas que [32], ampliadas con el coeficiente de expansión térmica y el coeficiente de expansión debida a la difusión. Las propiedades del material se resumen en la tabla ${\textstyle 1[10,32]}$.

Tabla 1: Propiedades material para los ejemplos de la sección 4 relativos a geles

Propiedad	Valor
${\textstyle \mu ^{0}}$	${\textstyle 0\mathrm {~J} /\mathrm {mol} }$
${\textstyle K}$	100 MPa
${\textstyle \Omega }$	${\textstyle 10^{-4}\mathrm {~m} ^{3}/\mathrm {mol} }$
D	${\textstyle 5\cdot 10^{-9}\mathrm {~m} ^{2}/\mathrm {s} }$
${\textstyle N_{r}}$	${\textstyle 2.43/\mathrm {m} ^{3}}$

(Continued)

Tabla 1 (continued)

Propiedad	Valor
${\textstyle G}$	1 MPa (a la temp. de ref.)
${\textstyle \chi _{L}}$	0.1 (adimensional)
${\textstyle \chi _{H}}$	0.7 (adimensional)
${\textstyle \Delta }$	5.0 K
${\textstyle c_{0}}$	${\textstyle 10^{7}\mathrm {~J} /\left(\mathrm {m} ^{3}\mathrm {~K} \right)}$
${\textstyle K}$	${\textstyle 0.5\mathrm {~W} /(\mathrm {m} \mathrm {K} )}$
${\textstyle \alpha }$	${\textstyle 10^{-5}\mathrm {~K} ^{-1}}$
${\textstyle \beta }$	${\textstyle 6\cdot 10^{-5}\mathrm {~mm} ^{3}/\mathrm {mol} }$

Se supondrá una dependencia térmica del parámetro de Flory-Huggins dada por la siguiente expresión:

${\textstyle \chi (\Theta )={\frac {1}{2}}\left(\chi _{L}+\chi _{H}\right)-{\frac {1}{2}}\left(\chi _{L}-\chi _{H}\right)\tanh \left({\frac {\Theta {-\Theta }_{T}}{\Delta }}\right)}$,

donde ${\textstyle \Theta _{T}}$ es la temperatura de transición del gel, ${\textstyle \chi _{L}}$ es el valor del parámetro de Flory-Huggins a temperaturas inferiores a la temperatura de transición del gel y ${\textstyle \chi _{H}}$ es el valor de dicho parámetro a temperaturas superiores y ${\textstyle \Delta }$ es la anchura de la zona de transición de temperatura entre ${\textstyle \chi _{L}}$ y ${\textstyle \chi _{H}}$.

A partir de la energía libre de Helmholtz propuesta, en primer lugar, puede obtenerse el segundo tensor de Piola-Kirchhoff:

${\textstyle {\boldsymbol {S}}=2{\frac {\partial \Psi }{\partial {\boldsymbol {C}}}}=G\left({\boldsymbol {I}}-{\boldsymbol {C}}^{-1}\right)+KJ^{s}\ln {\boldsymbol {J}}^{e}{\boldsymbol {C}}^{-1}-3\alpha K\left(\theta {-\theta }_{0}\right){\boldsymbol {C}}^{-1}-3\beta K\left(g-g_{0}\right){\boldsymbol {C}}^{-1}}$,

y a partir de éste, el primer tensor de Piola-Kirchhoff y el tensor de tensiones de Cauchy pueden obtenerse mediante las conocidas relaciones ${\textstyle {\boldsymbol {P}}={\boldsymbol {F}}{\boldsymbol {S}}}$ y ${\textstyle {\boldsymbol {\sigma }}=J^{-1}{\boldsymbol {F}}{\boldsymbol {S}}{\boldsymbol {F}}^{T}}$, quedando:

${\textstyle {\boldsymbol {P}}=G\left({\boldsymbol {F}}-{\boldsymbol {F}}^{T}\right)+{\boldsymbol {J}}^{s}K\mathbf {J} ^{\mathrm {e} }{\boldsymbol {F}}^{-T}-3\alpha K\left(\theta {-\theta }_{0}\right){\boldsymbol {F}}^{-T}-3\beta K\left(g-g_{0}\right){\boldsymbol {F}}^{-T}}$,

${\textstyle {\boldsymbol {\sigma }}=J^{-1}\left[G({\boldsymbol {B}}-{\boldsymbol {I}})+KJ^{s}\ln J^{e}{\boldsymbol {I}}-3\alpha K\left(\theta {-\theta }_{0}\right){\boldsymbol {I}}-3\beta K\left(g-g_{0}\right){\boldsymbol {I}}\right]}$.

El potencial químico se obtiene por derivación respecto a la concentración de la expresión de energía libre de Helmholtz propuesta:

${\textstyle \mu ={\frac {\partial \Psi }{\partial g}}=\mu _{0}+R\Theta \left[\ln(1-\phi )+\phi {+\chi }\phi ^{2}\right]+\Omega {\frac {1}{2}}K\left(\ln J^{e}\right)^{2}-K\Omega \ln J^{e}-3\beta K\ln J}$,

siendo:

${\textstyle \phi ={\frac {1}{1+\Omega g}}}$.

En lo que respecta a las derivadas segundas del potencial termodinámico correspondiente a la energía libre de Helmholtz, la parte simétrica del tensor de acoplamiento termo-mecánico, puede obtenerse como:

${\textstyle {\boldsymbol {M}}^{s}=\Theta {\frac {\partial ^{2}\Psi }{\partial \Theta \partial {\boldsymbol {C}}}}={\frac {\Theta }{2}}{\frac {1}{2}}{\frac {\partial {\boldsymbol {S}}}{\partial \Theta }}={\frac {\Theta }{2}}{\frac {\partial }{\partial \Theta }}\left[G\left({\boldsymbol {I}}-{\boldsymbol {C}}^{-1}\right)+KJ^{s}\ln J^{e}{\boldsymbol {C}}^{-1}\right]=}$

${\begin{aligned}&={\frac {\Theta }{2}}\left[{\frac {\partial G}{\partial \Theta }}\left({\boldsymbol {I}}-{\boldsymbol {C}}^{-1}\right)+K{\frac {\partial J^{s}}{\partial \Theta }}\ln J^{e}{\boldsymbol {C}}^{-1}+KJ^{s}{\frac {\partial \ln J^{e}}{\partial \Theta }}{\boldsymbol {C}}^{-1}\right]=\\&={\frac {\Theta }{2}}\left[{\frac {\partial G}{\partial \Theta }}\left({\boldsymbol {I}}-{\boldsymbol {C}}^{-1}\right)+K\Omega {\frac {\partial g}{\partial \Theta }}\ln J^{e}{\boldsymbol {C}}^{-1}+K{\boldsymbol {J}}^{s}\left({\frac {-\Omega {\frac {\partial g}{\partial \Theta }}}{(1+\Omega g)}}\right){\boldsymbol {C}}^{-1}\right].\end{aligned}}$

El tensor de acoplamiento termo-mecánico, ${\textstyle {\boldsymbol {M}}}$, se obtendrá mediante la ecuación (37) a partir de su parte simétrica (56).

En cuanto al tensor material de acoplamiento difusivo-mecánico, definido en la ecuación (38), y el coeficiente de acoplamiento difusivo-térmico, definido en la ecuación (39), se calcularán de forma numérica en la implementación computacional del modelo.

La figura 2 y la figura 3 muestran, respectivamente, los valores de la concentración y la tensión de von Mises. En todos los casos el fluido se inyecta en el fluido desde el exterior por la cara contraria a su base, a la que se han impedido los desplazamientos. Los tubos con geometrías 1 y 2 se pandean como resultado del hinchamiento presentado, 8 y 6 pétalos. El tubo con la geometría 3 es más grueso y no se pandea. Todo ello coincidiendo con los resultados experimentales [37] y una simulación de referencia [11].

Figura 2: De arriba a abajo, evolución de la concentración para ${\textstyle t=0,300,600,900\mathrm {~s} }$. De izquierda a derecha, geometrías 1,2 y 3 de la tabla 2 .

Figura 3: De arriba a abajo, evolución de la tensión de von Mises ${\textstyle [\mathrm {Pa} ]}$ para ${\textstyle t=0,300,600,900\mathrm {~s} }$. De izquierda a derecha, geometrías 1,2 y 3 de la tabla 2

4.2 Ejemplo: geles elastoméricos con tres campos acoplados

Para ilustrar el acoplamiento triple de los campos mecánico, térmico y difusivo se presenta un segundo ejemplo de los tres campos fuertemente acoplados describiéndose la respuesta del sólido como se indica en la sección 2.

En este ejemplo se ha utilizado una corona cilíndrica con las dimensiones correspondientes a la geometría 1 indicada en la tabla 2, las propiedades del material siguen siendo las indicadas en la tabla 1. En cuanto a las condiciones de contorno, se ha retenido el desplazamiento de los puntos de la superficie exterior impidiendo la deformación de la misma pero permitiendo los desplazamientos de los puntos interiores y se ha inyectado fluido por una de las bases mediante la imposición del potencial químico indicado en la expresión (47).

En la figura 4 y la figura 5 se muestran la temperatura y la tensión de von Mises en varios instantes de la simulación. En la figura 5 se ha practicado una sección radial en la corona cilíndrica para poder ver los puntos interiores y un factor de escala que permita apreciar su desplazamiento.

Tabla 2: Dimensiones tubo cilíndrico de diámetro exterior ${\textstyle D}$, altura ${\textstyle h}$ y espesor ${\textstyle t}$

Geometría	${\textstyle \mathrm {D} (\mathrm {mm} )}$	${\textstyle \mathrm {t} (\mathrm {mm} )}$	${\textstyle \mathrm {h} (\mathrm {mm} )}$
1	4.636	0.206	0.6
2	3.191	0.118	0.6
3	3.191	0.355	0.6

Figura 4: Evolución de la temperatura ${\textstyle [\mathrm {K} ]}$ para ${\textstyle t=0,100,200,300\mathrm {~s} }$

Figura 5: Evolución de la tensión de von Mises [Pa] para ${\textstyle t=0,100,200,300\mathrm {~s} }$. Se ha aplicado un factor de escala para apreciar la deformación de los puntos que no están situados en la superficie exterior

4.3 Ejemplo: con tres campos acoplados aplicado a baterías de ion litio

En aras de ilustrar la flexibilidad del modelo presentado para resolver problemas que involucran deformación mecánica, transmisión de calor y transporte de masa, con todos los campos acoplados, se planteará un tercer ejemplo, tomando como base el ejemplo del artículo [38] referente a baterías de ion litio, extendiéndolo para incorporar el campo térmico. Se analizará un cuarto de corona cilíndrica, cuya geometría y condiciones de contorno se muestran en la figura 6. Como condición de contorno se mantendrá un potencial prefijado en la superficie exterior utilizando una expresión similar a la empleada en [39]:

${\textstyle \mu (t)=\mu ^{0}+\left(\mu ^{0}-\mu ^{f}\right){\frac {t}{T}}}$,

de forma que el potencial químico aumenta linealmente desde un valor ${\textstyle \mu ^{0}}$ hasta un valor ${\textstyle \mu ^{f}}$, siendo ${\textstyle \mu ^{0}=-5.800\mathrm {~kJ} /\mathrm {mol} ,\mu ^{f}=5.800\mathrm {~kJ} /\mathrm {mol} }$ y ${\textstyle T=30}$ minutos.

Figura 6: Cuarto de corona cilíndrica: geometría y discretización

Para el modelado del comportamiento del material termo-difusivo-mecánico, se utilizará la siguiente energía libre de Helmholtz, construida como una expresión aditiva de una parte que establece la relación de acoplamiento térmico-difusivo ${\textstyle \left(\Psi ^{t-d}\right)}$, otra parte correspondiente al acoplamiento termo-mecánico ${\textstyle \left(\Psi ^{t-m}\right)}$ y otra relativa al acoplamiento difusivo-mecánico ( ${\textstyle \Psi ^{d-m}}$ ), además de las partes puramente térmica, mecánica y difusiva ( ${\textstyle \Psi ^{t},\Psi ^{m},\Psi ^{d}}$, respectivamente).

${\begin{aligned}\Psi _{R}=&\underbrace {\mu _{0}g} _{\Psi ^{d}}+\underbrace {{\frac {1}{2}}G\left(3\left({\bar {\lambda }}^{2}-1\right)-2\ln J\right)+J^{s}\left({\frac {1}{2}}K\left(\ln J^{e}\right)^{2}\right)} _{\Psi ^{m}}+\underbrace {c_{0}\left(\left(\Theta {-\Theta }_{0}\right)-\Theta \ln \left({\frac {\Theta }{\Theta }}\right)\right.} _{\psi ^{t}}\\&\times \underbrace {R\Theta g\left(\ln \left({\frac {\Omega g}{1+\Omega g}}\right)-1\right)} _{\Psi ^{t-d}}-\underbrace {3\alpha k\left(\Theta {-\Theta }_{0}\right)\ln J} _{\Psi ^{t-m}}-\underbrace {3\beta k\left(g-g_{0}\right)\ln J} _{\Psi ^{d-m}}.\end{aligned}}$

La expresión anterior es similar a la planteada en la expresión (48), con algunas salvedades. Por un lado, se ha considerado que el módulo elasticidad transversal no tiene dependencia de la temperatura para el rango térmico previsto en el problema en cuestión. Además, el término de acoplamiento térmico difusivo se ha visto modificado, desapareciendo el término que involucraba el parámetro de Flory-Huggins en el los ejemplos de geles, quedando para dicho acoplamiento una expresión similar a la utilizada en [11].

La expresión correspondientes al potencial químico quedaría:

${\textstyle \mu ={\frac {\partial \Psi }{\partial g}}=\mu _{0}+R\Theta [\ln(1-\phi )+\phi ]+\Omega {\frac {1}{2}}K\left(\ln J^{e}\right)^{2}-K\Omega \ln J^{e}-3\beta K\ln J}$.

Manteniéndose invariables respecto a los ejemplos anteriores las expresiones analíticas correspondientes a los tensores de Piola-Kirchhoff y tensor de tensiones de Cauchy ((51)-(53)) y a las

correspondientes a las derivadas segundas del la energía libre de Helmholtz. Las propiedades del material del material del electrodo pueden verse en la tabla 3 [38,40,41].

Tabla 3: Propiedades material

Propiedad	Valor
${\textstyle \mu ^{0}}$	${\textstyle 0\mathrm {~J} /\mathrm {mol} }$
${\textstyle K}$	${\textstyle 6.94\cdot 10^{4}\mathrm {MPa} }$
${\textstyle \Omega }$	${\textstyle 3.5\cdot 10^{-6}\mathrm {~m} ^{3}/\mathrm {mol} }$
D	${\textstyle 10^{-14}\mathrm {~m} ^{2}/\mathrm {s} }$
${\textstyle G}$	${\textstyle 3.97\cdot 10^{4}\mathrm {MPa} }$
${\textstyle c_{0}}$	${\textstyle 0.5\cdot 10^{7}\mathrm {~J} /\left(\mathrm {m} ^{3}\mathrm {~K} \right)}$
${\textstyle K}$	${\textstyle 0.5\mathrm {~W} /(\mathrm {m} \mathrm {K} )}$
${\textstyle \alpha }$	${\textstyle 10^{-5}\mathrm {~K} ^{-1}}$
${\textstyle \beta }$	${\textstyle 6\cdot 10^{-5}\mathrm {~mm} ^{3}/\mathrm {mol} }$

En la figura 7 puede verse la evolución de la deformación y de la temperatura. En la figura 8 puede verse la tensión de von Mises obteniéndose valores máximos de la tensión de von Mises similares a los presentados en diversos estudios [42-44].

Figura 7: Evolución de la temperatura ${\textstyle [\mathrm {K} ]}$ para ${\textstyle t=0,300,900,1000\mathrm {~s} }$

Figura 8: Evolución de la tensión de von Mises [Pa] para ${\textstyle t=0,300,900,1000\mathrm {~s} }$

5 Resumen y conclusiones

En este artículo se ha presentado una discretización por el método de los elementos finitos del problema de tres campos térmico-difusivo-mecánico. El problema incluye la transmisión de calor, el transporte de masa y el equilibrio mecánico, y todos los acoplamientos posibles han sido tenidos en cuenta. Además, la formulación es adecuada para grandes deformaciones.

A. Ortiz-Toranzo and I. Romero Discretización por elementos finitos del problema termo-difusivo-mecánico con grandes deformaciones, Rev. int. métodos numér. cálc. diseño ing. (2024). Vol.40, (4), 47

El modelo de tres campos se ha presentando dentro del marco de la termo-mecánica de medios continuos y se ha mostrado cómo las ecuaciones constitutivas derivan de una energía libre, confiriendo al planteamiento una gran generalidad. Después, el artículo describe la formulación por el método de los elementos finitos de una discretización espacial y temporal de los tres problemas acoplados. Como novedad adicional frente a otros trabajos de la literatura, se describe con detalle la implementación de las ecuaciones, incluyendo todos los térmicos de la linealización necesaria para su resolución por el método de Newton-Raphson.

Se ha empleado la formulación obtenida para simular problemas de geles elastoméricos y un ejemplo relacionado con baterías de ion litio. En este último ejemplo se consideran los procesos de difusión de iones en el interior del ánodo que lo deforman, provocando variaciones térmicas y del campo de tensiones.

El artículo proporciona los fundamentos teóricos para la descripción de los problemas termo-difusivo-mecánicos de medios continuos y su aproximación numérica. Este marco, junto con todos los detalles que se proporcionan para la implementación de sus ecuaciones, serán de gran utilidad para los investigadores que deseen desarrollar formulaciones de elementos finitos para este tipo problemas, unificando y complementando la información disponible en la literatura.

Un aspecto relevante que no se trata con profundidad en este trabajo, por limitar su extensión y sus objetivos, es la posibilidad de emplear métodos de resolución de ecuaciones no lineales distintos al de Newton-Raphson (cuasi-Newton, BFGS, gradiente conjugado no lineal, etc.). Los desarrollos de este trabajo y en particular los detalles del Apéndice Apéndice A pueden servir como guía para la formulación de dichos métodos más eficientes.

Agradecimientos/Acknowledgement: No aplica.

Fuentes de financiación/Funding Statement: I.R. agradece la financiación parcial de los proyectos TED2021-130255B-C32 y PID2021-128812OB-I00 del Ministerio de Ciencia e Innovación de España.

Aportaciones de los autores/Author Contributions: Ángel Ortiz-Toranzo: software, redacción del manuscrito, análisis, investigación. Ignacio Romero: conceptualización, metodología, revisión y edición del manuscrito, supervisión. Todos los autores revisaron los resultados y aprobaron la versión final del manuscrito.

Disponibilidad de datos y materiales/Availability of Data and Materials: Los autores proporcionarán información y datos sobre las simulaciones a aquellas personas que los soliciten por correo electrónico.

Consideraciones éticas/Ethics Approval: No aplican.

Conflictos de interés/Conflicts of Interest: Los autores declaran que no tienen intereses financieros que pudieran haber influido en el trabajo reportado en este artículo.

References

1. Pan Y, Zhong Z. A nonlinear constitutive model of unidirectional natural fiber reinforced composites considering moisture absorption. J Mech Phys Solids. 2014;69:132-42. doi:10.1016/j.jmps.2014.04.007.
2. Gangwar A, Kumar V, Yaylacı M, Panda S. Computational modelling and mechanical characteristics of polymeric hybrid composite materials: an extensive review. Arch Comput Methods Eng. 2024;1-21. doi:10.1007/s11831-024-10097-4.
3. Huang P, Guo Z. Li-ion distribution and diffusion-induced stress calculations of particles using an imagebased finite element method. Mech Mater. 2021;157(7):103843. doi:10.1016/j.mechmat.2021.103843.
4. Wang M, Ye H, Zhai C, Yang G. Stress-dependent chemo-mechanical performance of amorphous si anodes for li-ion batteries upon lithiation. ACS Appl Energy Mater. 2021;4(12):14718-26. doi:10.1021/ acsaem.1c03394.
5. Ateshian GA. On the theory of reactive mixtures for modeling biological growth. Biomech Model Mechanobiol. 2007;6:423-45.
6. Klemczak B, Smolana A, Jedrzejewska A. Modeling of heat and mass transfer in cement-based materials during cement hydration-a review. Energies. 2024;17(11):2513. doi:10.3390/en17112513.
7. Zhang M, Lin D, He Z, Yang R. Modeling the ionic diffusion coefficient of unsaturated hardened cement paste: a micromechanical approach. Cement Concr Res. 2024;177(5):107415. doi:10.1016/ j.cemconres.2023.107415.
8. Zhao H, Shan YL, Ding DJ, Shi LG, Wen DY, Cao WK, et al. Experiment study of stress and pore pressure in setting cement paste. In: ARMA US Rock Mechanics/Geomechanics Symposium, 2022; Santa Fe, New Mexico, USA, ARMA-2022. doi:10.56952/ARMA-2022-0386.
9. Elmukashfi E, Tarleton E, Cocks AC. A modelling framework for coupled hydrogen diffusion and mechanical behaviour of engineering components. Comput Mech. 2020;66(1):189-220. doi:10.1007/s00466-020-01847-9.
10. Chester SA, Anand L. A thermo-mechanically coupled theory for fluid permeation in elastomeric materials: application to thermally responsive gels. J Mech Phys Solids. 2011;59(10):1978-2006. doi:10.1016/ j.jmps.2011.07.005.
11. Romero I, Andrés EM, Ortiz-Toranzo Á. Variational updates for general thermo-chemomechanical processes of inelastic solids. Comput Methods Appl Mech Eng. 2021;385(1):114013. doi:10.1016/j.cma. 2021.114013.
12. Brunner F, Seidlhofer T, Ulz MH. A numerical model for chemo-thermo-mechanical coupling at large strains with an application to thermoresponsive hydrogels. Comput Mech. 2024;74(3):1-28. doi:10.1007/ s00466-024-02443-x.
13. Chan G, Mooney DJ. New materials for tissue engineering: towards greater control over the biological response. Trends Biotechnol. 2008;26(7):382-92. doi:10.1016/j.tibtech.2008.03.011.
14. Jayakumar A, Jose VK, Lee J-M. Hydrogels for medical and environmental applications. Small Methods. 2020;4(3):1900735. doi:10.1002/smtd. 201900735.
15. Chester SA. A constitutive model for coupled fluid permeation and large viscoelastic deformation in polymeric gels. Soft Matter. 2012;8(31):8223-33. doi:10.1039/c2sm25372k.
16. Chester SA, Anand L. A coupled theory of fluid permeation and large deformations for elastomeric materials. J Mech Phys Solids. 2010;58(11):1879-906. doi:10.1016/j.jmps.2010.07.020.
17. Jiang C, Zhong Z. A thermo-electro-chemo-mechanically coupling theory considering species diffusion and electrochemical reaction. Int J Solids Struct. 2023;267(17):112147. doi:10.1016/j.ijsolstr.2023.112147.
18. Keilegavlen E, Berge R, Fumagalli A, Starnoni M, Stefansson I, Varela J, et al. PorePy: an opensource software for simulation of multiphysics processes in fractured porous media. Comput Geosci. 2021;25(1):243-65. doi:10.1007/s10596-020-10002-5.
19. Ogita G, Matsumoto K, Mochizuki M, Mikami Y, Ito K. Evaluation of hydrogen-induced cracking behavior in duplex stainless steel by numerical simulation of stress and diffusible hydrogen distribution at the microstructural scale. ISIJ Int. 2021;61(4):1135-42.
20. Sheng Z, Prahl U. Numerical simulation of hydrogen distribution around a crack tip in a high-Mn steel. Mater Today Commun. 2021;28:102647. doi:10.1016/j.mtcomm.2021.102647.
21. Peradzyński Z. Diffusion of calcium in biological tissues and accompanying mechanochemical effects. Arch Mech. 2010;62(6):423-40.
22. Zhu X, Chen Y, Chen H, Luan W. The diffusion induced stress and cracking behaviour of primary particle for Li-ion battery electrode. Int J Mech Sci. 2020;178:105608. doi:10.1016/j.ijmecsci.2020.105608.
23. Malaeke H, Asghari M. A mathematical formulation for analysis of diffusion-induced stresses in micropolar elastic solids. Arch Appl Mech. 2023;93(8):3093-111. doi:10.1007/s00419-023-02427-y.
24. Wu Y, Guo Z-S. Concentration distribution and stresses in porous electrodes with particle-particle contact. J Electrochem Soc. 2021;168(9):090507.
25. Martins LS, Guimarães LF, Botelho Junior AB, Tenório JAS, Espinosa DCR. Electric car battery: an overview on global demand, recycling and future approaches towards sustainability. J Environ Manage. 2021;295:113091. doi:10.1016/j.jenvman.2021.113091.
26. Thakur AK, Prabakaran R, Elkadeem MR, Sharshir SW, Arıcı M, Wang C, et al. A state of art review and future viewpoint on advance cooling techniques for Lithium-ion battery system of electric vehicles. J Energy Storage. 2020;32:101771. doi:10.1016/j.est.2020.101771.
27. Morita Y, Saito Y, Yoshioka T, Shiratori T. Estimation of recoverable resources used in lithiumion batteries from portable electronic devices in Japan. Resour Conserv Recycling. 2021;175:105884. doi:10.1016/j.resconrec.2021.105884.
28. Wang S, Cheng T, Zhang YZ, Wu X, Xiao S, Lai W-Y. Deformable lithium-ion batteries for wearable and implantable electronics. Appl Phys Rev. 2022;9(4). doi:10.1063/5.0117252.
29. Clarke M, Alonso JJ. Lithium-ion battery modeling for aerospace applications. J Aircr. 2021;58(6):1323-35. doi:10.2514/1.C036209.
30. Stecca M, Ramirez Elizondo L, Batista Soeiro T, Bauer P, Palensky P. A comprehensive review of the integration of battery energy storage systems into distribution networks. IEEE Open J Indus Electron Soc. 2020;1:46-65. doi:10.1109/OJIES.2020.2981832.
31. Flory PJ. Statistical mechanics of swelling of network structures. J Chem Phy. 1950;18(1):108-11. doi:10.1063/1.1747424.
32. Chester SA, Di Leo CV, Anand L. A finite element implementation of a coupled diffusion-deformation theory for elastomeric gels. Int J Solids Struct. 2015;52:1-18. doi:10.1016/j.ijsolstr.2014.08.015.
33. Chester SA. Gel mechanics: a thermo-mechanically coupled theory for fluid permeation in elastomeric material. Procedia IUTAM. 2015;12(5):10-9. doi:10.1016/j.piutam.2014.12.003.
34. Gurtin ME, Fried E, Anand L. The mechanics and thermodynamics of continua. New York, USA: Cambridge University Press; 2010.
35. Holzapfel G. Nonlinear solid mechanics. In: A continuum approach for engineering. English. Second print, New York, USA: John Wiley & Sons; 2001.
36. Anand L. A thermo-mechanically-coupled theory accounting for hydrogen diffusion and large elasticviscoplastic deformations of metals. Int J Solids Struct. 2011;48(6):962-71. doi:10.1016/j.ijsolstr. 2010. 11.029.
37. Lee H, Zhang J, Jiang H, Fang NX. Prescribed pattern transformation in swelling gel tubes by elastic instability. Phys Rev Lett. 2012;108(21):214304. doi:10.1103/PhysRevLett.108.214304.
38. Dal H, Miehe C. Computational electro-chemo-mechanics of lithium-ion battery electrodes at finite strains. Computat Mech. 2015;55:303-25.
39. Di Leo CV, Rejovitzky E, Anand L. A Cahn-Hilliard-type phase-field theory for species diffusion coupled with large elastic deformations: application to phase-separating Li-ion electrode materials. J Mech Phys Solids. 2014;70(5):1-29. doi:10.1016/j.jmps.2014.05.001.
40. Zhang X, Chumakov S, Li X, Klinsmann M, Kim SU, Linder C, et al. An electro-chemo-thermomechanical coupled three-dimensional computational framework for lithium-ion batteries. J Electrochem Soc. 2020;167(16):160542.
41. Malave V. Chemo-elastic behavior of reconstructed Li-ion battery cathode particles with phase transformation: a numerical and analytical investigation (Tesis Doctoral). Colorado School of Mines: CO, USA; 2014.
42. Seo JH, Chung M, Park M, Han SW, Zhang X, Sastry AM. Generation of realistic particle structures and simulations of internal stress: a numerical/AFM Study of ${\textstyle \mathrm {LiMn} _{2}\mathrm {O} _{4}}$ Particles. ECS Meet Abstrs, 2010; MA2010-03:625. doi:10.1149/MA2010-03/1/625.
43. Zhang X, Sastry AM, Shyy W. Intercalation-induced stress and heat generation within single lithium-ion battery cathode particles. J Electrochem Soc. 2008;155(7):A542. doi:10.1149/1.2926617.
44. Xiao X, Wu W, Huang X. A multi-scale approach for the stress analysis of polymeric separators in a lithiumion battery. J Power Sourc. 2010;195(22):7649-60. doi:10.1016/j.jpowsour.2010.06.020.

Apéndice A. Linealización ecuaciones

La linealización de los términos de la ecuación (41), correspondiente al balance de la cantidad de movimiento, quedaría:

${\begin{aligned}\operatorname {Lin} \left[\left\langle {\boldsymbol {S}},{\boldsymbol {F}}^{T}\left(\sum _{a}\nabla N^{a}\delta {\boldsymbol {\varphi }}^{a}\right)\right\rangle \right]=&\left\langle {\frac {\partial {\boldsymbol {S}}}{\partial {\boldsymbol {C}}}}2{\boldsymbol {F}}^{T}\left(\sum _{b}\nabla N^{b}\Delta {\boldsymbol {\varphi }}^{b}\right),{\boldsymbol {F}}^{T}\left(\sum _{a}\nabla N^{a}\delta {\boldsymbol {\varphi }}^{a}\right)\right\rangle +\\&+\left\langle {\frac {\partial {\boldsymbol {S}}}{\partial {\boldsymbol {\Theta }}}}\left(\sum _{b}N^{b}\Delta \Theta ^{b}\right),{\boldsymbol {F}}^{T}\left(\sum _{a}\nabla N^{a}\delta {\boldsymbol {\varphi }}^{a}\right)\right\rangle +\\&+\left\langle {\frac {\partial {\boldsymbol {S}}}{\partial \mu }}\left(\sum N^{b}\Delta \mu ^{b}\right),{\boldsymbol {F}}^{T}\sum _{a}\nabla N^{a}\delta {\boldsymbol {\varphi }}^{a}\right\rangle +\left\langle {\boldsymbol {S}},\Delta {\boldsymbol {F}}^{T}\left(\sum _{a}\nabla N^{a}\delta {\boldsymbol {\varphi }}^{a}\right)\right\rangle .\end{aligned}}$

Teniendo en cuenta la definición del tensor de acoplamiento difusivo-mecánico ${\textstyle {\boldsymbol {L}}}$ (véase (38)), y del tensor de acoplamiento termo-mecánico ${\textstyle {\boldsymbol {M}}}$ (según la expresión (35)), la linealización de la expresión (60) puede expresarse:

${\begin{aligned}&\underbrace {\sum _{a,b}\delta {\boldsymbol {\varphi }}^{a}\left\langle {\boldsymbol {F}}\mathbb {C} {\boldsymbol {F}}^{T}\nabla N^{a},\nabla N^{b}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}} _{{\boldsymbol {\varphi }}-{\boldsymbol {\varphi }}}+\underbrace {\sum _{a,b}\delta {\boldsymbol {\varphi }}^{a}\left\langle N^{b},{\frac {\boldsymbol {M}}{\Theta }}\nabla N^{a}\right\rangle \Delta \Theta ^{b}} _{{\boldsymbol {\varphi }}-{\boldsymbol {b}}}+\\&\underbrace {\sum _{a,b}\delta {\boldsymbol {\varphi }}^{a}\left\langle \nabla N^{a},{\boldsymbol {S}}\nabla N^{b}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}} _{{\boldsymbol {\varphi }}-{\boldsymbol {\varphi }}}+\underbrace {\sum _{a,b}\delta {\boldsymbol {\varphi }}^{a}\left\langle N^{b},{\boldsymbol {L}}\nabla N^{a}\right\rangle \Delta \mu ^{b}} _{\varphi {-\mu }}.\end{aligned}}$

La linealización de la ecuación (42) correspondiente al balance de la energía resulta en los siguientes bloques:

${\begin{aligned}\operatorname {Lin} \left[\left\langle \Theta {\dot {s}},\sum _{a}N^{a}\delta \Theta ^{a}\right\rangle \right]=&\left\langle \sum _{a}N^{a}\delta \Theta ^{a},\Delta \Theta {\dot {s}}\right\rangle +\left\langle \sum _{a}N^{a}\delta \Theta ^{a},\Theta \Delta {\dot {s}}\right\rangle =\\=&\left\langle \sum _{a}N^{a}\delta \Theta ^{a},\Delta \Theta {\dot {s}}\right\rangle +\\&+\left\langle \sum _{a}N^{a}\delta \Theta ^{a},\Theta {\frac {1}{\Delta t}}\left[-{\frac {\partial ^{2}{\mathfrak {G}}}{\partial {\boldsymbol {F}}\partial \Theta }}\Delta {\boldsymbol {F}}-{\frac {\partial ^{2}{\mathfrak {G}}}{\partial ^{2}\Theta }}\Delta \Theta -{\frac {\partial ^{2}{\mathfrak {G}}}{\partial \Theta \partial \mu }}\Delta \mu \right]\right\rangle =\end{aligned}}$

${\begin{aligned}=&\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle N^{a},N^{b}{\dot {S}}\right\rangle \Delta \Theta ^{b}} _{\Theta {-\Theta }}+\underbrace {\sum _{\Theta {-\varphi }}} _{{}_{a,b}\delta \Theta ^{a}\left\langle N^{a},\Theta {\frac {1}{\Delta t}}{\frac {\partial s}{\partial {\boldsymbol {F}}}}\nabla N^{b}\right\rangle \delta {\boldsymbol {\varphi }}^{b}}+\\&+\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle N^{a},\Theta {\frac {1}{\Delta t}}{\frac {\partial s}{\partial \Theta }}N^{b}\right\rangle \Delta \Theta ^{b}} _{\Theta {-}b}+\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle N^{a},\Theta {\frac {1}{\Delta t}}{\frac {\partial s}{\partial \mu }}N^{b}\right\rangle \Delta \mu ^{b}} _{\Theta {-\mu }}\end{aligned}}$

${\textstyle \operatorname {Lin} \left[\left\langle -{\boldsymbol {Q}},\sum _{a}\nabla N^{a}\delta \Theta ^{a}\right\rangle \right]=\left\langle \Delta \left({\boldsymbol {F}}^{-1}K_{0}{\boldsymbol {F}}^{-T}\right)\mathrm {GRAD} \Theta ,\sum _{a}\nabla N^{a}\delta \Theta ^{a}\right\rangle =}$

${\begin{aligned}&\left\langle {\boldsymbol {F}}^{-1}K_{0}{\boldsymbol {F}}^{-T}\operatorname {GRAD} \Delta \Theta ,\sum _{a}\nabla N^{a}\delta \Theta ^{a}\right\rangle +\\&\left\langle \Delta \left({\boldsymbol {F}}^{-1}K_{0}{\boldsymbol {F}}^{-T}\right)\operatorname {GRAD\Theta ,\sum _{a}\nabla N^{a}\delta \Theta ^{a}\rangle =} \right.\\=&\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle \nabla N^{a},{\boldsymbol {F}}^{-1}K_{0}{\boldsymbol {F}}^{-T}\nabla N^{b}\right\rangle \Delta \Theta ^{b}+} _{\Theta {-\Theta }}\\&\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle \nabla N^{a},K_{0}{\frac {\partial {\boldsymbol {C}}^{-1}}{\partial {\boldsymbol {C}}}}2{\boldsymbol {F}}^{T}\nabla N^{b}\operatorname {GRAD\Theta } \right\rangle \Delta {\boldsymbol {\varphi }}^{b}} _{\Theta -{\boldsymbol {\varphi }}}.\end{aligned}}$

${\textstyle \operatorname {Lin} \left[\left\langle {\boldsymbol {h}}\cdot \operatorname {GRAD} (\mu ),\left(\sum _{a}N^{a}\delta \Theta ^{a}\right)\right\rangle \right]=}$

${\begin{aligned}&\left\langle \left({\frac {\partial {\boldsymbol {h}}}{\partial {\boldsymbol {C}}}}\Delta {\boldsymbol {C}}+{\frac {\partial {\boldsymbol {h}}}{\partial \mu }}\Delta \mu +{\frac {\partial {\boldsymbol {h}}}{\partial \Theta }}\Delta \Theta \right)\operatorname {GRAD} (\mu )+{\boldsymbol {h}}\operatorname {GRAD} (\Delta \mu ),\sum _{a}N^{a}\delta \Theta ^{a}\right\rangle =\\&=\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {gD}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\operatorname {GRAD} (\mu ){\boldsymbol {F}}^{-1}\nabla N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}} _{\Theta {-\varphi }}-\end{aligned}}$

$-\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {2D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {\partial g}{\partial {\boldsymbol {C}}}}\operatorname {GRAD} (\mu ){\boldsymbol {F}}^{T}\nabla N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta \varphi ^{b}} _{\Theta {-\varphi }}-$

$-\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {\partial g}{\partial \mu }}\operatorname {GRAD} (\mu )N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta \mu ^{b}} _{\Theta {-\mu }}-$

$-\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {Dg}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\nabla N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta \mu ^{b}} _{\Theta {-\mu }}-$

${\begin{aligned}&-\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\operatorname {GRAD} (\mu ){\frac {\partial g}{\partial \Theta }}N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta \Theta ^{b}} _{\Theta {-\Theta }}+\\&+\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {Dg}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {1}{\Theta }}\operatorname {GRAD} (\mu )N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta \Theta ^{b}} _{\Theta {-\Theta }}-\\&-\underbrace {\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {Dg}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\operatorname {GRAD} (\mu )\nabla N^{b},N^{a}\right\rangle \Delta \mu ^{b}} _{\Theta {-\mu }}.\end{aligned}}$

En cuanto a la linealización de los términos de la ecuación correspondiente a la conservación de masa (43), los términos resultantes son:

${\textstyle \operatorname {=} \underbrace {\left\langle {\dot {g}},\left(\sum _{a}N^{a}\delta \mu ^{a}\right)\right\rangle } _{1}-\underbrace {\left\langle {\boldsymbol {h}},\left(\sum _{a}\nabla N^{a}\delta \mu ^{a}\right)\right\rangle } _{2}}$.

La linealización del término 1 de la expresión anterior puede escribirse como

${\begin{aligned}\left\langle \Delta {\dot {g}},\sum _{a}N^{a}\delta \mu ^{a}\right\rangle =&\left\langle {\frac {\partial {\dot {g}}}{\partial {\boldsymbol {C}}}}\Delta {\boldsymbol {C}},\sum _{a}N^{a}\delta \mu ^{a}\right\rangle +\\&+\left\langle {\frac {\partial {\dot {g}}}{\partial \Theta }}\Delta \Theta ,\sum _{a}N^{a}\delta \mu ^{a}\right\rangle +\left\langle {\frac {\partial {\dot {g}}}{\partial \mu }}\Delta \mu ,\sum _{a}N^{a}\delta \mu ^{a}\right\rangle \end{aligned}}$

En cuanto al primer término del miembro de la derecha de la expresión (66), se tiene:

$\left\langle {\frac {\partial {\dot {g}}}{\partial {\boldsymbol {C}}}}\Delta {\boldsymbol {C}},\sum _{a}N^{a}\delta \mu ^{a}\right\rangle =\underbrace {\sum _{a,b}\delta \mu ^{a}\left\langle f{\frac {\partial g}{\partial {\boldsymbol {F}}}}\nabla N^{b},N^{a}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}} _{\mu -{\boldsymbol {\varphi }}}$

En cuanto al segundo término del miembro de la derecha de la expresión (66), se tiene:

$\left\langle {\frac {\partial {\dot {g}}}{\partial \Theta }}\Delta \Theta ,\sum _{a}N^{a}\delta \mu ^{a}\right\rangle =\underbrace {\sum _{a,b}\delta \mu ^{a}\left\langle f{\frac {\partial g}{\partial \Theta }}N^{b},N^{a}\right\rangle \Delta \Theta ^{b}} _{\mu {-\Theta }}.$

En cuanto al tercer término del miembro de la derecha de la expresión (66), se tiene:

$\left\langle {\frac {\partial {\dot {g}}}{\partial \mu }}\Delta \mu ,\sum _{a}N^{a}\delta \mu ^{a}\right\rangle =\underbrace {\sum _{a,b}\delta \mu ^{a}\left\langle f{\frac {\partial g}{\partial \mu }}N^{b},N^{a}\right\rangle \Delta \mu ^{b}} _{\mu {-\mu }}.$

La linealización del término 2 de la expresión (65) puede escribirse:

${\begin{aligned}&-\left\langle \Delta {\boldsymbol {h}},\left(\sum _{a}\nabla N^{a}\delta \mu ^{a}\right)\right\rangle =-\left\langle {\frac {\partial h}{\partial {\boldsymbol {C}}}}\Delta {\boldsymbol {C}}+{\frac {\partial h}{\partial \mu }}\Delta \mu +{\frac {\partial h}{\partial \Theta }}\Delta \Theta ,\left(\sum _{a}\nabla N^{a}\delta \mu ^{a}\right)\right\rangle =\\&=\underbrace {-\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {gD}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\operatorname {GRAD} (\mu ){\boldsymbol {F}}^{-1}\nabla N^{b},\nabla N^{a}\right\rangle \Delta \varphi ^{b}} _{\mu {-\varphi }}+\\&+\underbrace {\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {2D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {\partial g}{\partial {\boldsymbol {C}}}}\operatorname {GRAD} (\mu ){\boldsymbol {F}}^{T}\nabla N^{b},\nabla N^{a}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}} _{\mu -{\boldsymbol {\varphi }}}+\\&+\underbrace {\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {\partial g}{\partial \mu }}\operatorname {GRAD} (\mu )N^{b},\nabla N^{a}\right\rangle \Delta \mu ^{b}} _{\mu {-\mu }}+\\&+\underbrace {\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}g\nabla N^{b},\nabla N^{a}\right\rangle \Delta \mu ^{b}} _{\mu {-\mu }}+\\&+\underbrace {\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\operatorname {GRAD} (\mu ){\frac {\partial g}{\partial \Theta }}N^{b},\nabla N^{a}\right\rangle \Delta \Theta ^{b}} _{\mu {-\Theta }}-\\&-\underbrace {\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {Dg}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {1}{\Theta }}\operatorname {GRAD} (\mu )N^{b},\nabla N^{a}\right\rangle \Delta \Theta ^{b}} _{\mu {-\Theta }}.\end{aligned}}$

Las ecuaciones linealizadas correspondientes al balance de la cantidad de movimiento, balance de la energía y balance de masa, pueden escribirse matricialmente como un sistema de ecuaciones ${\textstyle {\mathcal {K}}{\boldsymbol {d}}={\boldsymbol {F}}}$, siendo ${\textstyle {\mathcal {K}}}$ la matriz tangente, ${\textstyle {\boldsymbol {d}}}$ el vector de desplazamientos y ${\textstyle {\boldsymbol {F}}}$ el vector de cargas, esto es:

$\left[{\begin{array}{ccc}K_{\varphi ,\varphi }&K_{\varphi ,\mu }&K_{\varphi ,\Theta }\\K_{\mu ,\varphi }&K_{\mu ,\mu }&K_{\mu ,\Theta }\\K_{\Theta ,\varphi }&K_{\Theta ,\mu }&K_{\Theta ,\Theta }\end{array}}\right]\left[{\begin{array}{l}d_{\varphi }\\d_{\mu }\\d_{\Theta }\end{array}}\right]=\left[{\begin{array}{l}{\boldsymbol {F}}_{\varphi }\\{\boldsymbol {F}}_{\mu }\\{\boldsymbol {F}}_{\Theta }\end{array}}\right],$

siendo:

${\begin{aligned}K_{\varphi ,\varphi }=&\sum _{a,b}\delta {\boldsymbol {\varphi }}^{a}\left\langle {\boldsymbol {F}}\mathbb {C} {\boldsymbol {F}}^{T}\nabla N^{a},\nabla N^{b}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}\\&+\sum _{a,b}\delta {\boldsymbol {\varphi }}^{a}\left\langle \nabla N^{a}{\boldsymbol {S}},\nabla N^{b}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}\\K_{\varphi ,\Theta }=&\sum _{a,b}\delta {\boldsymbol {\varphi }}^{a}\left\langle N^{b},{\frac {\boldsymbol {M}}{\Theta }}\nabla N^{a}\right\rangle \Delta \Theta ^{b}\end{aligned}}$

${\begin{aligned}&K_{\varphi ,\mu }=\sum _{a,b}\delta {\boldsymbol {\varphi }}^{a}\left\langle N^{b},{\boldsymbol {L}}\nabla N^{a}\right\rangle \Delta \mu ^{b}\\&K_{\mu ,\mu }=\sum _{a,b}\delta \mu ^{a}\left\langle f{\frac {\partial g}{\partial \mu }}N^{b},N^{a}\right\rangle \Delta \mu ^{b}\\&+\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {\partial g}{\partial \mu }}\operatorname {GRAD} (\mu )N^{b},\nabla N^{a}\right\rangle \Delta \mu ^{b}\\&+\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}g\nabla N^{b},\nabla N^{a}\right\rangle \Delta \mu ^{b}\\&K_{\mu ,\varphi }=\sum _{a,b}\delta \mu ^{a}\left\langle f{\frac {\partial g}{\partial {\boldsymbol {F}}}}\nabla N^{b},N^{a}\right\rangle \Delta \varphi ^{b}\\&-\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {gD}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\operatorname {GRAD} (\mu ){\boldsymbol {F}}^{-1}\nabla N^{b},\nabla N^{a}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}\\&+\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {2D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {\partial g}{\partial {\boldsymbol {C}}}}\operatorname {GRAD} (\mu ){\boldsymbol {F}}^{T}\nabla N^{b},\nabla N^{a}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}\\&K_{\mu ,\Theta }=\sum _{a,b}\delta \mu ^{a}\left\langle f{\frac {\partial g}{\partial \Theta }}N^{b},N^{a}\right\rangle \Delta \Theta ^{b}\\&+\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\operatorname {GRAD} (\mu ){\frac {\partial g}{\partial \Theta }}N^{b},\nabla N^{a}\right\rangle \Delta \Theta ^{b}\\&-\sum _{a,b}\delta \mu ^{a}\left\langle {\frac {Dg}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {1}{\Theta }}\operatorname {GRAD} (\mu )N^{b},\nabla N^{a}\right\rangle \Delta \Theta ^{b}\\&K_{\ominus ,\Theta }=-\sum _{a,b}\delta \Theta ^{a}\left\langle N^{a},N^{b}{\dot {s}}\right\rangle \Delta \Theta ^{b}\\&+\sum _{a,b}\delta \Theta ^{a}\left\langle N^{a},\Theta {\frac {1}{\Delta t}}{\frac {\partial s}{\partial \Theta }}N^{b}\right\rangle \Delta \Theta ^{b}\\&+\sum _{a,b}\delta \Theta ^{a}\left\langle \nabla N^{a},{\boldsymbol {F}}^{-1}K_{0}{\boldsymbol {F}}^{-T}\nabla N^{b}\right\rangle \Delta \Theta ^{b}-\\&-\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\operatorname {GRAD} (\mu ){\frac {\partial g}{\partial \Theta }}N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta \Theta ^{b}\\&+\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {Dg}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {1}{\Theta }}\operatorname {GRAD} (\mu )N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta \Theta ^{b}\end{aligned}}$

${\begin{aligned}K_{\Theta ,\varphi }=&\sum _{a,b}\delta \Theta ^{a}\left\langle N^{a},\Theta {\frac {1}{\Delta t}}{\frac {\partial s}{\partial {\boldsymbol {F}}}}\nabla N^{b}\right\rangle \delta {\boldsymbol {\varphi }}^{b}\\&-\sum _{a,b}\delta \Theta ^{a}\left\langle \nabla N^{a},K_{0}{\frac {\partial {\boldsymbol {C}}^{-1}}{\partial {\boldsymbol {C}}}}2{\boldsymbol {F}}^{T}\nabla N^{b}\operatorname {GRAD} \Theta \right\rangle \Delta {\boldsymbol {\varphi }}^{b}\\&+\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {gD}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\operatorname {GRAD} (\mu ){\boldsymbol {F}}^{-1}\nabla N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}\\&-\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {2D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {\partial g}{\partial {\boldsymbol {C}}}}\operatorname {GRAD} (\mu ){\boldsymbol {F}}^{T}\nabla N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta {\boldsymbol {\varphi }}^{b}\\K_{\Theta ,\mu }=&\sum _{a,b}\delta \Theta ^{a}\left\langle N^{a},\Theta {\frac {1}{\Delta t}}{\frac {\partial s}{\partial \mu }}N^{b}\right\rangle \Delta \mu ^{b}\\&-\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {D}{R\Theta \operatorname {det} {\boldsymbol {F}}}}{\frac {\partial g}{\partial \mu }}\operatorname {GRAD} (\mu )N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta \mu ^{b}\\&-\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {Dg}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\nabla N^{b}\operatorname {GRAD} (\mu ),N^{a}\right\rangle \Delta \mu ^{b}\\&-\sum _{a,b}\delta \Theta ^{a}\left\langle {\frac {Dg}{R\Theta \operatorname {det} {\boldsymbol {F}}}}\operatorname {GRAD} (\mu )\nabla N^{b},N^{a}\right\rangle \Delta \mu ^{b}\end{aligned}}$