Introduction¶
fdfault is a finite difference code for numerical simulation of elastodynamic fracture and friction problems in 2 and 3 dimensions, principally those arising in study of dynamic earthquake rupture. The code solves the elastic or elastic-plastic wave equation in the bulk material, coupled to frictional failure on the fault and external boundary conditions.
Governing Equations¶
Material Governing Equations¶
The code solves the elastodynamic wave equation in 2 or 3 dimensions, with either an elastic or viscoplastic bulk material. For a 3D continuum, momentum balance requires that
where \({\rho}\) is density. Additionally, a constitutive law relates the stresses to the deformations. For a homogeneous, isotropic elastic-plastic material, these take the form
where the elastic tensor is given by Hooke’s Law for a homogeneous isotropic material
The elastic tensor includes two material parameters: \({\lambda}\) is the first Lamé parameter and \({G}\) is the shear modulus. Summation over repeated indices is implied. For an elastic material, \({\dot{\epsilon}^{pl}_{ij}=0}\), while for a viscoplastic material the plastic strains are determined by the Drucker-Prager viscoplastic flow rule described below. For 2D problems, the code can handle either in-plane (\({v_x}\), \({v_y}\), \({\sigma_{xx}}\), \({\sigma_{xy}}\), \({\sigma_{yy}}\) are nonzero) or anti-plane (\({v_z}\), \({\sigma_{xz}}\), \({\sigma_{yz}}\) are nonzero) problems, where all variables are functions of \({x}\) and \({y}\).
Initial conditions must be provided for the velocities and stresses. The code assumes that the velocities are initially zero, and provides numerous ways to set the initial stress depending on the complexity of the problem. More details are provided when describing the code specifics.
For Drucker-Prager Viscoplasticity, plastic flow occurs when stresses exceed the yield function \({F(\sigma_{ij})}\):
where \({\bar{\tau}=\sqrt{s_{ij}s_{ij}/2}}\) is the second invariant of the deviatoric stress tensor \({s_{ij}=\sigma_{ij}-(\sigma_{kk}/3)\delta_{ij}}\), \({c}\) is related to the cohesion, and \({\mu}\) is related to the internal coefficient of friction. For viscoplasticity, flow is allowed to exceed the yield criterion according to
where \({\Lambda=\sqrt{2\dot{e}^{pl}_{ij}\dot{e}^{pl}_{ij}}}\) is the equivalent plastic strain rate from the deviatoric plastic strain rate \({\dot{e}^{pl}_{ij}=\dot{\epsilon}^{pl}_{ij}-(\dot{\epsilon}^{pl}_{kk}/3)\delta_{ij}}\) and \({\eta}\) is a viscoplastic “viscosity” defining the time scale over which stresses can exceed the yield stress. If stresses are accumulated at a rate faster than the relaxation time of the viscoplastic material, the material behaves elastically, and the stress then decays towards the yield surface if no further stresses are applied.
The components of plastic flow are determined by
with \({P_{ij}(\sigma_{ij})=s_{ij}/(2\bar{\sigma})+(\beta/3)\delta_{ij}}\), where the \({\beta}\) parameter determines the ratio of volumetric to plastic strain. Thus, viscoplastic materials are determined by \({\rho}\), \({\lambda}\), \({G}\), \({c}\), \({\mu}\), \({\beta}\), and \({\eta}\), while elastic materials require specification of \({\rho}\), \({\lambda}\), and \({G}\). Rate independent plasticity arises in the limit that \({\eta \rightarrow 0}\), though the equations become increasingly stiff in this limit and plastic strain can exhibit localization that is not resolved by the spatial grid in many problems.
The code can handle variable material properties in two different ways. The first type is “block-like” structures, where the material properties are piecewise constant in an arbitrary number of blocks. Alternatively, one can specify the elastic properties in a point-by-point fashion (the plasticity material properties must be block-like even in the case of continuously varying elastic properties in the present version of the code).
Regardless of the method used, the domain is composed of a regular grid of blocks with conforming edges, though the block boundaries can have complex shapes, so this does not limit the ability of the code to handle complex earth structures. Blocks are coupled together through interfaces, which can either be locked or frictional. Locked interfaces require continuity of velocity and forces across block edges, while frictional interfaces allow for relative slip across the interface according to several possible friction laws:
Frictionless interfaces do not support shear stresses, so any stress applied to such interfaces result in shear slip.
Kinematic forcing allows for a forced rupture following a prescribed set of rupture times.
Slip weakening interfaces follow a slip-dependent friction law, where the friction coefficient has static and dynamic values and transitions between the two according to a function that decreases linearly with slip. Slip weakening laws also allow for cohesion. Slip weakening can also be combined with kinematic forcing to nucleate a rupture.
Shear Transformation Zone (STZ) Theory interfaces have a frictional strength that depends on the slip rate and a dynamic state variable representing the configurational disorder in the fault gouge.
Specific details on each of the friction models is provided below. In addition to these base friction laws, a future version of the code will allow you to specify arbitrary slip- or rate- and state-dependent friction laws through the Python module, which automatically generates the required source code for new friction laws.
External boundaries (i.e. boundaries not between two blocks) can have absorbing, free surface (traction-free), or rigid (velocity-free) boundary conditions.
Frictional Descriptions¶
General Formulation¶
Friction laws are used to set the boundary conditions on the edges connecting block interfaces. At each edge, the velocities and stresses for each block are rotated into normal and tangential components. The material governing equations above provide one set of relationships between the velocities and tractions; an additional condition is required to fully specify the boundary values.
For the normal components, the code requires continuity of velocities and tractions across the interface. In the event that tensile normal stresses occur, the code sets the normal traction to zero but does not allow for the interface to open.
For the shear components, the tractions are required to be continuous, but the velocities can be discontinuous across the fault. The velocity discontinuity, called the slip velocity \({V}\) with vector components \({V_1}\) and \({V_1}\), is the primary quantity of interest, as it describes the rate at which fault slip occurs. The wave equation provides one relationship between the shear traction and the slip velocity, and the other is given by another equation (along with the constraint that the traction and slip are parallel to one another). This relationship may simply specify the traction (so that the slip velocity can be solved for directly), or require a nonlinear solver if the relationship is nonlinear and cannot be solved in closed form.
The slip velocity is integrated in time, resulting in slip components \({U_1}\) and \({U_2}\), along with a scalar slip \({U}\) that is computed as a line integral. The slip is used by some friction laws to control the evolution of strength as a function of time, but is computed for all frictional descriptions.
Frictionless Interface¶
The frictionless interface is the simplest friction description – the interface cannot support a shear traction, so the slip velocity is set to be whatever value is needed to ensure that no stress accumulates. This applies to both components of traction for 3D problems. This model requires no additional parameter specifications beyond the material properties.
Kinematic Forcing¶
Rupture propagation can be prescribed using kinematic forcing. Kinematic forcing requires specification of 4 parameters: a static friction coefficient \({\mu_s}\), a dynamic friction coefficient \({\mu_d}\), a time scale \({t_c}\) which sets the time scale over which friction linearly weakens from static to dynamic, and a rupture time \({t_{rup}}\) that determines when the frictional weakening initiates at a given point. The friction coefficient can be determined directly from the time, and the friction coefficient combined with the normal traction determines the shear traction on the fault. If the value of the shear traction set by the wave equation is less than the value set by the friction law, the fault is locked, the slip rate is zero, and the shear traction takes the value from the wave equation.
Linear Slip-Weakening¶
In order for rupture propagation to be truly spontaneous, a friction law that does not prescribe rupture time is required. The most common form used in dynamic rupture modeling is the linear slip-weakening law. As with kinematic forcing, a static and dynamic friction coefficient must be specified. However, instead of a rupture time and a weakening time, the weakening process is characterized by a slip weakening distance \({d_c}\). Friction weakens linearly with slip from the static to dynamic value:
Once the friction coefficient is known, the shear traction is set in a similar fashion to the Kinematic Forcing law described above. The code also allows for frictional cohesion \({c_0}\), in which case the shear traction \({\tau}\) is:
where \({\sigma_n}\) is the normal traction (negative in compression).
For 3D problems with vector slip, each vector component of the velocity/traction is solved separately, but the total slip \({U}\) is used to determine the weakening behavior.
Additionally, the code allows for a combination kinematic/slip-weakening law, where the code uses the minimum friction coefficient that is calculated for the kinematic and slip-weakening laws. This is used in cases where the rupture is initiated with a kinematic procedure, but then is allowed to propagate spontaneously.
Shear Transformation Zone Theory¶
Shear Transformation Zone (STZ) Theory is a rate- and state-dependent constitutive law, which ties the fault strength to the dynamic evolution of a state variable representing the effective disorder temperature \({\chi}\). Frictional strength \({\mu}\) is determined by the slip rate and the effective temperature:
Note that this cannot be solved in closed form for the friction coefficient, so the code solves this equation simultaneously with the elastic wave equation for \({\mu}\) and \({V}\). Additionally, the effective temperature evolves in time according to
where the rate-dependent steady state effective temperature is \({\hat{\chi}(V)=\chi_w/\log(V_1/V)}\). The STZ model introduces several additional parameters: a reference slip rate \({V_0}\), an activation barrier for slip \({f_0}\), the frictional direct effect \({a}\), the friction coefficient at jamming \({\mu_y}\), the effective temperature specific heat \({c_0}\), the normalized effective temperature activation barrier \({\chi_w}\), the reference slip rate for STZ activation \({V_1}\), the STZ relaxation rate \({R}\), and the normalized effective relaxation barrier \({\beta}\).
More details on the STZ model and the parameters can be found in the papers listed below.
Numerical Details¶
The code solves the governing equations numerically using finite differences. The outer boundaries of each block is described by a series of 6 surfaces (4 curves in 2D), and each block is transformed from physical space to the unit cube (unit square in 2D). The governing equations are transformed as well, and the code solves the resulting problem on a structured grid using high order finite differences. The grid is generated using standard transfinite interpolation, and the required metric derivatives for solving the governing equations and applying boundary conditions are automatically calculated using finite differences. Grids between neighboring blocks must be conforming, though no other continuity condition is required across block interfaces. The grid must satisfy certain smoothness constraints (these are checked during the initialization steps in the code), though non-uniform grid spacing along interfaces is allowed, provided that the resulting grid meets the smoothness requirements. Boundary conditions at external boundaries and interfaces are applied in locally rotated normal/tangential coordinate systems using characteristic variables.
The specific finite difference operators used exhibit a summation-by-parts property that mimics the properties of integration by parts. This allows for estimates of the energy dissipation rate of the numerical scheme. Boundary conditions are imposed weakly using the Simultaneous Approximation Term approach, and this combined with the summation by parts difference operators allows for a provably stable numerical scheme that matches the energy dissipation rate of the continuous problem.
The code allows for central finite difference operators that are globally second, third, or fourth order accurate. Time integration is performed with a low memory Runge-Kutta method, with either first, second, third, or fourth order accuracy in time. Artificial dissipation can also be used for the finite difference operators, which will reduce the numerical artifact oscillations that can occur with large grid spacings.
The finite difference method is only applied to the elastic part of the problem. The plasticity equations are handled through an operator splitting procedure, where the elastic problem is solved first and then used as initial conditions for the plasticity problem. This is done using an implicit backward Euler method.
For more details on the numerical methods used, please consult the papers listed below.