Composed friction partition

src/porepy/compositional/materials.py

@@ -427,11 +427,74 @@ class SolidConstants(Constants):
 
 
 @dataclass(kw_only=True, eq=False)
-class FractureDamageSolidConstants(SolidConstants):
-    """Solid parameters for fracture damage models."""
+class AsperityContactSolidConstants(SolidConstants):
+    """Solid parameters for friction composed from asperity contact.
+
+    See :class:`~porepy.constitutive_laws.AsperityStressPartition`. The basic
+    friction coefficient and the dilation angle the law also needs are inherited
+    from :class:`SolidConstants`.
+
+    """
 
     # NOTE this makes a deep copy of the solid constants dict.
     SI_units: ClassVar[dict[str, str]] = dict(**SolidConstants.SI_units)
+    SI_units.update(
+        {
+            "ploughing_friction_coefficient": "-",
+            "transitional_normal_traction": "Pa",
+            "stress_partition_exponent": "-",
+        }
+    )
+
+    ploughing_friction_coefficient: float = 0.0
+    r"""Ploughing friction coefficient :math:`\mu_p^0` [-].
+
+    The frictional resistance of asperities being sheared through, in the limit where
+    all contact is doing so. It enters the composed friction coefficient as
+    :math:`a_s \mu_p^0 d^f`, so it is the value reached at
+    :math:`\sigma_n = \sigma_T` with the asperities intact; see
+    :meth:`~porepy.constitutive_laws.FractureDamage.friction_coefficient`.
+
+    The default of zero leaves the ploughing term absent, reducing the composed
+    coefficient to the sliding law alone.
+
+    """
+
+    transitional_normal_traction: float = float("inf")
+    r"""Transitional normal traction :math:`\sigma_T` [Pa].
+
+    The normal traction at which the asperities are fully sheared through rather than
+    slid over, so that the stress partition
+    :meth:`~porepy.constitutive_laws.AsperityStressPartition.stress_partition`
+    reaches one. Ladanyi and Archambault (1970) identify it with the strength of the
+    asperity material, hence a value of the order of the uniaxial compressive strength.
+
+    The default is infinite, which holds :math:`a_s \equiv 0` and thereby leaves the
+    partition inert: all contact is then sliding contact, which is the low-traction
+    limit of the law.
+
+    """
+
+    stress_partition_exponent: float = 1.5
+    r"""Exponent :math:`K` of the stress partition [-].
+
+    Sets how abruptly contact transfers from sliding to shearing as the normal traction
+    approaches :math:`\sigma_T`. The default of 1.5 is the value of Ladanyi and
+    Archambault (1970).
+
+    Values in :math:`(1, 2)` are non-integer, so the partition's base must be clipped to
+    the non-negative branch before it is raised; see
+    :meth:`~porepy.constitutive_laws.AsperityStressPartition.stress_partition`.
+
+    """
+
+
+@dataclass(kw_only=True, eq=False)
+class FractureDamageSolidConstants(AsperityContactSolidConstants):
+    """Solid parameters for fracture damage models."""
+
+    # NOTE this makes a deep copy of the solid constants dict.
+    SI_units: ClassVar[dict[str, str]] = dict(**AsperityContactSolidConstants.SI_units)
     SI_units.update(
         {
             "residual_dilation_damage": "-",
@@ -448,8 +511,9 @@ class FractureDamageSolidConstants(SolidConstants):
 
     :math:`d_0^f = 1` holds :math:`d^f \equiv 1` and thereby switches friction damage
     off, leaving the intact coefficient; this is how an undamaged reference run is
-    configured. :math:`d_0^f` is a floor on the friction coefficient, since the damage
-    state multiplies it in full.
+    configured. :math:`d_0^f` floors the ploughing term alone, not the whole friction
+    coefficient, so :math:`d_0^f = 0` is admissible: it says that asperities worn flat
+    cease to ploughing-resist, while the basic friction of the rock surfaces remains.
 
     """
 
@@ -471,21 +535,25 @@ class FractureDamageSolidConstants(SolidConstants):
 
     The frictional work per unit area that must be dissipated against the asperities for
     the dilation damage state to decay by a factor :math:`e` towards its residual; see
-    :class:`~porepy.constitutive_laws.DilationDamage`. It absorbs Archard's
+    :class:`~porepy.constitutive_laws.FractureDamage`. It absorbs Archard's
     dimensionless wear coefficient and is therefore an effective calibrated quantity
     rather than a purely geometric one.
 
-    Associated with the longer-wavelength waviness of the fracture surface, which is the
-    scale that produces resolvable aperture change.
+    Together with :attr:`friction_wear_energy_scale` it sets how fast the two roles of
+    one asperity population degrade: this one governs the loss of the surfaces' ability
+    to force climbing, and hence of resolvable aperture change.
 
     """
 
     friction_wear_energy_scale: float = 1.0
     r"""Wear energy scale for friction damage :math:`\Lambda_c^f` [J * m^-2].
 
-    As :attr:`dilation_wear_energy_scale`, but for the friction channel, and associated
-    with the shorter-wavelength unevenness, which contributes frictional resistance
-    without geometrically resolvable dilation.
+    As :attr:`dilation_wear_energy_scale`, but governing the loss of the asperities'
+    resistance to being sheared through, i.e. of the ploughing term. Which of the two
+    roles an asperity plays at a given moment is decided by the normal traction through
+    :meth:`~porepy.constitutive_laws.AsperityStressPartition.stress_partition`,
+    not by the asperity's size, so the two scales are not a decomposition of the
+    roughness into wavelength bands.
 
     """
 

src/porepy/examples/fracture_damage.py

@@ -153,6 +153,9 @@ def collect_data(self) -> Any:
 class FractureDamageMomentumBalance(  # type: ignore[misc]
     pp.models.solution_strategy.ContactIndicators,
     DamageDataSaving,
+    pp.constitutive_laws.FractureDamage,
+    pp.constitutive_laws.AsperityStressPartition,
+    pp.constitutive_laws.DilationRotatedFriction,
     pp.constitutive_laws.FractureDamageEvolutionCoefficients,
     TimeDependentDamageBCs,
     pp.MomentumBalance,
@@ -180,23 +183,14 @@ class FractureDamageHistoryMixin(
 ):
     """The damage history variable and its convolution equation.
 
-    Mixed in whenever any damage channel is active, since the history is common to
-    them. The channels themselves are added by the constitutive mixins in
-    :data:`damage_types`.
+    Separate from the constitutive laws in
+    :class:`~porepy.constitutive_laws.FractureDamage` so that the history, which is
+    common to both channels, is supplied exactly once.
     """
 
     pass
 
 
-# Collect the damage types in a dictionary for easy access when building models with
-# different regimes. These supply the constitutive laws only; the history they read is
-# provided once by :class:`FractureDamageHistoryMixin`.
-damage_types = {
-    "dilation": pp.constitutive_laws.DilationDamage,
-    "friction": pp.constitutive_laws.FrictionDamage,
-}
-
-
 class ExactSolution:
     """Exact solution for the damage model.
 
@@ -461,6 +455,7 @@ def oriented_length(j):
         # dilation damage without incurring too much normal opening and stress.
         "maximum_elastic_fracture_opening": 0.0,  # [m] Simplify by assuming no elastic
         # opening.
+        "fracture_gap": 0.0,  # [m] Mated aperture.
     }
 )
 # Increase shear modulus to suppress shear displacements relative to normal ones.
@@ -526,7 +521,9 @@ def create_displacement_controlled_setup(
     Parameters:
         isotropic: If True, use isotropic damage length; otherwise anisotropic.
         dim: Spatial dimension of the bulk domain (2 or 3).
-        damages: Damage types to include (subset of {"dilation", "friction"}).
+        damages: Damage channels to activate (subset of {"dilation", "friction"}).
+            Channels left out have their residual damage state pinned at one, which
+            holds them intact.
 
     Returns:
         Tuple ``(model_class, model_params, solver_params)`` with caller-owned mutable
@@ -542,8 +539,6 @@ def create_displacement_controlled_setup(
         params["exact_solution"] = ExactSolutionAnisotropic
         model_class = add_mixin(damage.AnisotropicFractureDamageLength, model_class)
 
-    for name in damages:
-        model_class = add_mixin(damage_types[name], model_class)
     model_class = add_mixin(FractureDamageHistoryMixin, model_class)
 
     geom = (
@@ -553,18 +548,22 @@ def create_displacement_controlled_setup(
 
     displacements = north_displacements_3d.copy()
     displacements = displacements[:dim]
-    # Keep compression for the first steps and open the fracture in the last one.
-    displacements[1] = 0.98e-3
-    displacements[1, 4] = 3e-3
+    # Keep the fracture closed for the first steps and open it in the last one.
+    displacements[1] = -2e-5
+    displacements[1, 4] = 2e-3
 
     params.update(
         {
             "time_manager": pp.TimeManager(np.arange(0.0, 5.0), 1.0, True),
             "north_displacements": displacements,
         }
     )
+    solid = solid_params.copy()
+    for name in ("dilation", "friction"):
+        if name not in damages:
+            solid[f"residual_{name}_damage"] = 1.0
     params["material_constants"] = {
-        "solid": FractureDamageSolidConstants(**solid_params.copy()),  # type: ignore[arg-type]
+        "solid": FractureDamageSolidConstants(**solid),  # type: ignore[arg-type]
     }
 
     solver_params = {