Fully mixed finite element methods for the coupling of viscoelasticity and reaction-diffusion models

1 Introduction

• The fully mixed formulation treats elastic and viscous stresses, velocity, spin, and calcium flux as primary unknowns, giving direct access to coupling quantities and the diffusive flux.The symmetry of total stress is imposed weakly; the mixed diffusion formulation also provides local mass conservation.Source: §1 · p. 2
• Separate elastic and viscous stresses are retained, while only their sum requires H(div)-conformity, and the nonlinear coupling acts in both directions.Calcium enters the constitutive law through a Hill function, while the elastic-stress trace enters the calcium source term.Source: §1 · pp. 2–3
• A reduction to an integro-differential problem yields a unique global continuous solution, after which the remaining mixed variables are recovered uniquely.The local existence interval can be chosen independently of initial data, allowing the construction to be iterated.Source: §1 · p. 3
• The general semidiscrete framework is well posed for conforming finite-element spaces under a single inf-sup condition, with a discrete stability estimate.Eliminating algebraic variables from the differential-algebraic system gives an ordinary differential equation for the remaining variables.Source: §1 · p. 3
• An interpolation operator that commutes with divergence and preserves weak symmetry pairing supports optimal-order a priori error estimates.This structure also permits the spin variable to be eliminated from the main energy argument.Source: §1 · p. 3
• The Arnold–Falk–Winther, Farhloul–Fortin, and PEERS families satisfy the abstract assumptions, and the Farhloul–Fortin family admits superconvergent local displacement post-processing.The post-processing exploits the higher approximation order of the stress space.Source: §1 · p. 3
• Two- and three-dimensional experiments confirm predicted convergence rates, reproduce Kelvin–Voigt creep and cyclic loading, and illustrate mechanochemical effects on calcium wavetrains.The experiments also examine how mechanochemical feedback influences the wavetrains.Source: §1 · p. 3

2 Governing equations

• The coupling is bidirectional: calcium contributes active stress through a Hill-type function, while local dilation feeds into the calcium source.The dilation term can be rewritten using the trace of the elastic stress.Source: §2 · pp. 4, 6
• The Kelvin–Voigt mechanical equations couple momentum balance to elastic and viscous stresses, with zero velocity on Γu and zero total traction on Γσ.The total stress is the sum of the elastic and viscous stresses.Source: §2 · pp. 4–5
• Calcium and inactivated inositol triphosphate evolve through a reaction–diffusion system, with calcium diffusion represented by a flux in the mixed formulation.The calcium equation includes a source term proportional to the trace of the elastic stress.Source: §2 · pp. 5–6
• The analysis assumes globally Lipschitz reaction terms, although the stated physical reaction function may be singular and is not globally Lipschitz.The authors note that bounded solutions away from singular values permit Lipschitz behavior on the solution range; regularisation is another possible route.Source: §2 · pp. 5–6
• The weak formulation treats elastic and viscous stresses, velocity, spin, diffusive flux, calcium, and iIT as unknowns, with weak stress symmetry and mixed boundary conditions.The calcium flux is introduced as an additional unknown.Source: §2 · p. 6

3.1 The reduced integro-differential equations

• The reduction closes the dynamics in velocity and the two concentrations; displacement is recovered by time integration and stresses from the constitutive relations.The initial displacement enters the reduced system as a forcing offset.Source: §3 · p. 6, §3.1 · p. 7
• The weak formulation gives velocity and calcium time derivatives in dual spaces, making their initial conditions meaningful through the Lions–Magenes lemma.The reduced solution has velocity and calcium in L2-in-time energy spaces, while the inactivated inositol triphosphate variable is in H1-in-time with values in L2.Source: §3.1 · pp. 7–8

3.2 Well-posedness of the reduced problem

• On a sufficiently short interval, the coupled reduced problem is uniquely solvable by a contraction argument after freezing its coupling terms.The local existence time depends on coercivity and Lipschitz constants and model parameters, but not on initial data.Source: §3.2 · pp. 10–12
• A unique solution of the reduced coupled problem exists on the full interval (0, t_end], with the stated time-dependent regularity for velocity and the two concentrations.The local existence time is independent of the initial data, allowing the local result to be iterated over the full interval.Source: §3.2 · pp. 10, 12, 13
• For u0, v0 ∈ V and φ0 ∈ Q, the reduced solution satisfies ˙v, ˙φ ∈ L2(0, tend; L2(Ω)).This enhanced regularity supports reconstruction of the mixed variables.Source: §3.2 · pp. 13–14

3.3 Well-posedness of the mixed system

• Under the hypotheses of Theorem 3.6, the reduced-system solution and reconstructed mixed variables solve the mixed system (2.12).The mixed variables are reconstructed from the reduced solution in the distributional sense.Source: §3.3 · p. 15
• The mixed-system stability result implies at most one solution for fixed body force and initial data.The stability constant depends on the model parameters and final time, but not on the data.Source: §3.3 · pp. 16–17

• The discrete equations couple elastic and viscous stress relations with the weak stress–velocity–spin pairing and the calcium equation’s flux and elastic-stress terms.Source: §4 · pp. 17–18
• The discrete initial data are projections of the continuous data; no initial values are prescribed for σ2, spin, or flux because the algebraic relations determine them.Source: §4 · p. 18
• The differential variables are σ1, velocity, and the two concentrations, while σ2, spin, and flux are algebraic unknowns that can be eliminated to obtain an ODE.Source: §4 · pp. 19–21
• Under (H1), Assumption 3.1, and b ∈ C0([0, tend]; L2(Ω)), the semidiscrete problem has a unique solution on [0, tend].Source: §4 · p. 20
• The compatibility condition gives the stress–spin block full row rank, and positive definiteness makes the algebraic variables uniquely recoverable from the differential ones.Source: §4 · p. 20
• A discrete energy bound and Grönwall’s inequality prevent blow-up before tend, extending the ODE solution across the full interval.Source: §4 · p. 21
• With the weaker assumption b ∈ L2(0, tend; L2(Ω)), the solution remains unique and absolutely continuous, with its time-derivative regularity relaxed accordingly.Source: §4 · p. 21

5.1 Projection operators and error splitting

• The analysis assumes divergence-matching finite element spaces and compatible stress and flux projections that commute with divergence and preserve the weak-symmetry pairing.These properties support the error splitting and the orthogonality relations used to simplify the discrete error equations.Source: §5.1 · p. 22

5.2 The error equations

• Subtracting the semidiscrete equations from the continuous equations and inserting the projection-error splitting yields equations for the discrete errors.The discrete initial errors follow from the chosen L2-projected initial data.Source: §5.2 · p. 23

5.3 A priori error estimate

• Theorem 5.1 bounds the discrete errors by projection and initial-error contributions, with a constant independent of h.The initial stress-error contribution is of approximation order for the stated projected initial data; choosing σ1,h(0) = 0 removes that term.Source: §5.3 · pp. 24, 26
• Under Theorem 5.1’s hypotheses, Theorem 5.2 provides a mesh-independent a priori error estimate for the spin tensor.Its constant additionally depends on the inf-sup constant βd.Source: §5.3 · p. 26
• With additional regularity and bounded-derivative assumptions, Theorem 5.3 establishes mesh-independent error estimates in the divergence norms.The proof differentiates the error system in time; the resulting estimates use bounds on the time derivatives of the calcium and inactivated-IP3 variables.Source: §5.3 · pp. 27–28
• Estimate (5.27) follows by taking L2(0, t; •)-norms in (5.28) and (5.29) and inserting the bounds obtained.Source: §5.3 · p. 28
• The initial value θ(0) is obtained by evaluating (2.12) and (4.2) at t = 0, and is of approximation order provided discrete initial data are chosen compatibly with the equations at t = 0.Source: §5.3 · p. 28

6.1 Stable families of finite element spaces and convergence rates

• Arnold–Falk–Winther, Farhloul–Fortin, and PEERS spaces satisfy the abstract stability assumptions, including an h-independent inf-sup condition and a weak-symmetry-preserving interpolation property.The Farhloul–Fortin family approximates stress one order more accurately than PEERS, enabling displacement post-processing without full elliptic regularity.Source: §6.1 · p. 29
• The divergence errors for the combined stress and diffusive flux are O(h^(k+1)) under the additional regularity assumptions of the divergence-norm estimate.The estimate applies to div(eσ1 + eσ2) and div ez in L2-in-time norms.Source: §6.1 · p. 30
• For Farhloul–Fortin spaces, velocity and concentrations determine the rates, while stress and, with BDM diffusion spaces, flux gain one approximation order; the estimates also hold for Arnold–Falk–Winther and PEERS.For Arnold–Falk–Winther, the spin rate is hk+1 rather than hk+2; PEERS approximates stress at hk+1.Source: §6.1 · p. 30

6.2 Displacement superconvergence for specific finite element pairs

• The local displacement reconstruction is defined elementwise by matching the discrete displacement’s P0 moments, using the continuous relation ∇u = A1σ1 + ρ.The section studies this reconstruction alongside the projected-displacement superconvergence estimate.Source: §6.2 · p. 31
• For Farhloul–Fortin spaces paired with Raviart–Thomas diffusion spaces, the projected displacement converges one order faster than its h^(k+1) best-approximation order.The superconvergence result is established for Farhloul–Fortin viscoelasticity spaces with the stated Raviart–Thomas diffusion choice.Source: §6.2 · pp. 31–32

7.1 Corroboration of convergence against smooth solutions

• The predicted optimal convergence rates are observed for the Arnold–Falk–Winther, Farhloul–Fortin and PEERS discretisations, with component and total errors decreasing at rates set by polynomial approximation orders.The reported momentum and mass defects are small.Source: Table 7.1 · p. 33, §7.1 · p. 34, Table 7.2 · p. 34

7.2 Verification of the displacement superconvergence

• For the Farhloul–Fortin pair with k = 0, the post-processed displacement converges asymptotically at rate k + 2 = 2, while velocity and reconstructed displacement converge at rate k + 1 = 1.Post-processing into one degree above the usual velocity space also gives AFW and PEERS rates close to k + 2 for these smooth solutions, although the theorem does not cover those pairs.Source: §7.2 · p. 36, Table 7.4 · p. 36

7.3 Simple viscoelasticity benchmarks

• The creep test shows displacement rising toward its elastic limit under load, viscous stress decaying as elastic stress increases, and all quantities relaxing after unloading.Recovery after unloading is governed by a single time constant.Source: §7.3 · p. 37, Figure 7.2 · p. 37
• Under oscillatory loading, the uncoupled viscoelastic block reaches a periodic response with a phase lag between elastic and viscous branches and an open stress–strain hysteresis loop.The loop's enclosed area quantifies energy dissipated per cycle by the dashpot.Source: §7.3 · pp. 37–38, Figure 7.4 · p. 38

7.4 Mechanically driven pattern formation and calcium waves

• In the Atri case, target waves propagate at approximately 0.18 with period 6.23; one-way coupling leaves calcium unchanged while tissue area oscillates with the wavetrain.The area oscillation has about 0.1% peak-to-peak amplitude on a slow drift.Source: §7.4 · pp. 38–39, Figure 7.5 · p. 39, Figure 7.6 · p. 40
• For β = 1, 2, 5, the wavetrain period decreases from 6.23 to 5.78, 5.42 and 4.60, while peak concentration increases from 1.91 to 2.03, 2.05 and 2.21, respectively.For β = 10, positive feedback prevents a proper rest state.Source: §7.4 · p. 40
• The viscosity sweep shows area-change oscillation amplitude decreasing from 8% to 0.12%, while its lag approaches the quarter period of the purely viscous limit.For s ≥ 0.03, the lag follows the homogeneous Kelvin–Voigt estimate within 5%; the viscosity acts as a low-pass filter on mechanochemical feedback.Source: §7.4 · p. 40, Figure 7.7 · p. 41
• For β = 2, the two formulations agree within 0.001 periods in lag and 0.6% in area-oscillation amplitude.The comparison is between the mixed and primal formulations.Source: §7.4 · p. 40