Constitutive Models¶

Viscous Flow Models¶

ViscousFlowModel¶

class underworld3.constitutive_models.ViscousFlowModel[source]¶

Bases: Constitutive_Model

Viscous flow constitutive model for Stokes-type solvers.

Defines the relationship between deviatoric stress and strain rate:

\[\tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l} + \frac{\partial u_l}{\partial x_k} \right]\]

where \(\eta\) is the viscosity, which can be a scalar constant, SymPy function, Underworld mesh variable, or any valid combination. This results in an isotropic (but not necessarily homogeneous or linear) relationship between \(\tau\) and the velocity gradients.

Parameters:

unknowns (Unknowns) – The solver unknowns (typically velocity and pressure fields).
material_name (str, optional) – Name identifier for this material (used in multi-material setups).

Examples

>>> import underworld3 as uw
>>> stokes = uw.systems.Stokes(mesh)
>>> viscous = uw.constitutive_models.ViscousFlowModel(stokes.Unknowns)
>>> viscous.Parameters.shear_viscosity_0 = 1e21  # Pa.s
>>> stokes.constitutive_model = viscous

Diffusion Models¶

DiffusionModel¶

class underworld3.constitutive_models.DiffusionModel[source]¶

Bases: Constitutive_Model

Diffusion (Fourier/Fick) constitutive model for scalar transport.

Defines the flux-gradient relationship for scalar diffusion:

\[q_{i} = \kappa_{ij} \frac{\partial \phi}{\partial x_j}\]

For isotropic diffusion, \(\kappa_{ij} = \kappa \delta_{ij}\).

Parameters:

unknowns (Unknowns) – The solver unknowns (the scalar field being diffused).
material_name (str, optional) – Name identifier for this material.

Examples

>>> diffusion = uw.constitutive_models.DiffusionModel(poisson.Unknowns)
>>> diffusion.Parameters.diffusivity = 1e-6  # m^2/s
>>> poisson.constitutive_model = diffusion

TransverseIsotropicFlowModel¶

class underworld3.constitutive_models.TransverseIsotropicFlowModel[source]¶

Bases: ViscousFlowModel

Transversely isotropic (anisotropic) viscous flow model.

\[\tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l} + \frac{\partial u_l}{\partial x_k} \right]\]

where \(\eta\) is the viscosity tensor defined as:

\[\eta_{ijkl} = \eta_0 \cdot I_{ijkl} + (\eta_0-\eta_1) \left[ \frac{1}{2} \left[ n_i n_l \delta_{jk} + n_j n_k \delta_{il} + n_i n_l \delta_{jk} + n_j n_l \delta_{ik} \right] - 2 n_i n_j n_k n_l \right]\]

and \(\hat{\mathbf{n}} \equiv \{n_i\}\) is the unit vector defining the local orientation of the weak plane (a.k.a. the director).

The Mandel constitutive matrix is available in viscous_model.C and the rank-4 tensor form is in viscous_model.c.

Examples

>>> viscous_model = TransverseIsotropicFlowModel(dim)
>>> viscous_model.material_properties = viscous_model.Parameters(
...     eta_0=viscosity_fn,
...     eta_1=weak_viscosity_fn,
...     director=orientation_vector_fn
... )
>>> solver.constitutive_model = viscous_model
>>> tau = viscous_model.flux(gradient_matrix)
---

__init__(unknowns, material_name=None)[source]¶

Initialize a constitutive model.

Parameters:

unknowns (UnknownSet) – The solver’s unknowns (velocity, pressure, etc.)
material_name (str, optional) – A distinguishing name for this material’s symbols. If provided, symbols will be subscripted: η → η_{name} Useful when bundling multiple models in MultiMaterialModel.

property viscosity¶: Whatever the consistutive model defines as the effective value of viscosity in the form of an uw.expression

property K¶: Whatever the consistutive model defines as the effective value of viscosity in the form of an uw.expression

property grad_u¶: Symmetric strain rate tensor (with 1/2 factor).

\[\dot{\varepsilon}_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)\]

Viscoelastic Models¶

ViscoElasticPlasticFlowModel¶

class underworld3.constitutive_models.ViscoElasticPlasticFlowModel[source]¶

Bases: ViscousFlowModel

Viscoelastic-plastic flow constitutive model.

The stress (flux term) is given by:

\[\tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l} + \frac{\partial u_l}{\partial x_k} \right]\]

where \(\eta\) is the viscosity, a scalar constant, SymPy function, Underworld mesh variable, or any valid combination. This results in an isotropic (but not necessarily homogeneous or linear) relationship between \(\tau\) and the velocity gradients. You can also supply \(\eta_{IJ}\), the Mandel form of the constitutive tensor, or \(\eta_{ijkl}\), the rank-4 tensor.

The Mandel constitutive matrix is available in viscous_model.C and the rank 4 tensor form is in viscous_model.c. Apply the constitutive model using:

__init__(unknowns, order=1, integrator='bdf', material_name=None)[source]¶

Construct a viscoelastic-plastic flow model.

Parameters:

unknowns (Unknowns) – The solver unknowns (velocity, pressure).
order (int, default 1) –
Time-integration order. Combines with integrator:
- integrator='bdf', order=1: BDF-1 (backward Euler).
- integrator='bdf', order=2: BDF-2.
- integrator='etd', order=1: ETD-1 (single-step, fully L-stable, recommended default for VEP+yield).
- integrator='etd', order=2: ETD-2 (single-step with linear-quadrature forcing history; accurate on smooth VE but unstable in tight-yield VEP — produces global runaway, see EXPONENTIAL_VE_INTEGRATOR.md lessons #7, #9).
integrator (str, default "bdf") –
Time-integration scheme:
- "bdf": backward differentiation formula on the deviatoric-stress rate equation. Production default.
- "etd": exponential time-differencing — integrates the Maxwell relaxation operator analytically (α = exp(-Δt/τ)). order=1 is the recommended default for new code: same stability as BDF-1, exact handling of the relaxation factor at large Δt/τ. order=2 adds linear quadrature on the forcing history for higher accuracy on smooth VE (4.3× more accurate than BDF-2 on bench_ve_harmonic) but blows up under active yield in tight-yield TI faults. See docs/developer/design/EXPONENTIAL_VE_INTEGRATOR.md.
material_name (str, optional) – Name identifier for this material.

property order¶: Time integration order (1 or 2).

property effective_order¶

Effective order accounting for DDt history startup.

During the first few timesteps, the DDt may not have enough history to support the requested order. This property returns the lower of the requested order and the DDt’s effective order (which ramps from 1 to self.order as history accumulates).

property stress_history_ddt_kwargs¶

SemiLagrangian DDt kwargs based on integrator selection.

ETD-2 (order=2) needs the forcing-history slot; BDF and ETD-1 (order=1) do not.

property stress_star¶: Previous timestep stress \(\boldsymbol{\sigma}^*\) from history.

property stress_2star¶: Second-order stress history \(\boldsymbol{\sigma}^{**}\) (for 2nd order integration).

property E_eff¶

Effective strain rate including elastic-history coupling.

For BDF integration:

\[\dot{\varepsilon}_\mathrm{eff} = \dot{\varepsilon} - \sum_i c_i \frac{\sigma^{*(i)}}{2 \mu \Delta t}\]

For ETD-2 (exponential) integration:

\[\dot{\varepsilon}_\mathrm{eff} = (1-\varphi)\,\dot{\varepsilon} + \frac{\alpha}{2\eta}\,\sigma^* + (\varphi-\alpha)\,\dot{\varepsilon}^*\]

Both forms reduce to bare ε̇ when no elastic history is active. The yield criterion η_pl = τ_y/(2|E_eff|_II) is the same expression structure for both — it adapts naturally to the integrator’s E_eff.

property E_eff_inv_II¶

\(\dot{\varepsilon}_{II} = \sqrt{\frac{1}{2}\dot{\varepsilon}_{ij}\dot{\varepsilon}_{ij}}\).

Type:: Second invariant of effective strain rate

property K¶: Effective stiffness parameter (viscosity for visco-elastic-plastic flow).

property viscosity¶

Effective viscosity combining visco-elastic and plastic limits.

The yield mode controls how η_ve and η_pl are combined:

"smooth" (default): corrected harmonic η_ve·(1+f)/(1+f+f²) where f = η_ve/η_pl. Converges to η_pl at deep yielding, no Min/Max discontinuities.
"harmonic": 1/(1/η_ve + 1/η_pl). Smooth but undershoots τ_y when η_ve is small relative to η_pl.
"min": sharp Min(η_ve, η_pl). Exact yield stress but can cause SNES divergence with higher-order BDF time integration.

The unclipped η_ve depends on self._integrator — _unclipped_ve_viscosity returns the correct base for BDF or ETD-2 (raw η for ETD-2 since the time factor lives in E_eff; ve_effective_viscosity for BDF).

plastic_correction()[source]¶: Scaling factor to reduce stress to yield surface: \(f = \tau_y / \tau_{II}\).

property flux¶: Computes the effect of the constitutive tensor on the gradients of the unknowns. (always uses the c form of the tensor). In general cases, the history of the gradients may be required to evaluate the flux. For viscoelasticity, the

stress_projection()[source]¶: viscoelastic stress projection (no plastic response)

stress()[source]¶

Viscoelastic(-plastic) deviatoric stress for the weak form.

Both BDF and ETD-2 are written as σ = 2·viscosity·E_eff. E_eff carries the integrator-specific elastic-history coupling, and viscosity returns the appropriate yield- wrapped effective viscosity (ve_effective_viscosity for BDF, raw η for ETD-2 since the time factor is in E_eff).

property yield_mode¶

How to combine VE and plastic viscosities.

"softmin" (default): smooth approximation to Min —: η_ve / g(f) where g(f) ≈ max(1, f) with smoothing parameter δ (yield_softness, default 0.1). Approaches exact Min as δ → 0; smooth derivatives at the kink. Recommended default: gets within ~2 % of the true yield surface while avoiding the SNES kink penalties of "min".
"harmonic": parallel blending — 1/(1/η_ve + 1/η_pl).: Smooth but undershoots τ_y for soft materials.
"min": sharp cutoff — Min(η_ve, η_pl).: Exact yield but can cause SNES divergence with BDF-2 and BDF-2 phase-lag at BC discontinuities (see benchmarks).

Note: the previous "smooth" mode (corrected harmonic η_ve·(1+f)/(1+f+f²)) was retired — it under-clipped the yield surface by ~50 % under realistic forcing, with no compensating benefit over softmin. Recover from git history if needed (commit message keyword: smooth_yield).

property yield_softness¶

Regularisation parameter δ for "softmin" yield mode.

Controls how closely the soft minimum approximates the sharp Min. Smaller values → sharper yield (closer to Min, less robust). Larger values → smoother transition (more robust, lower stress).

Default 0.1. Only used when yield_mode == "softmin". Increase toward 0.5 if SNES convergence is difficult at yield onset.

property requires_stress_history¶: VEP models always require stress history tracking.

property plastic_fraction¶

1 - η_vep / η_ve.

Type:: Fraction of strain rate that is plastic

property is_elastic¶: True if elastic behavior is active (finite dt_elastic and shear_modulus).

property is_viscoplastic¶: True if plastic yielding is active (finite yield_stress).

Darcy Flow¶

DarcyFlowModel¶

class underworld3.constitutive_models.DarcyFlowModel[source]¶

Bases: Constitutive_Model

Darcy flow constitutive model for porous media flow.

Relates the Darcy flux to pressure gradients and body forces:

\[q_{i} = \kappa_{ij} \left( \frac{\partial p}{\partial x_j} - s_j \right)\]

where \(\kappa\) is the permeability (or hydraulic conductivity), \(p\) is the pressure (or hydraulic head), and \(s\) is the body force term (e.g., gravity: \(s = \rho g\)).

Parameters:

unknowns (Unknowns) – The solver unknowns (the pressure/head field).
material_name (str, optional) – Name identifier for this material.

Examples

>>> darcy = uw.constitutive_models.DarcyFlowModel(solver.Unknowns)
>>> darcy.Parameters.permeability = 1e-12  # m^2
>>> darcy.Parameters.s = [0, -rho * g]  # Gravity in y-direction
>>> solver.constitutive_model = darcy