1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5'''
8Analytical crack solutions for surface displacements and fault dislocations.
10'''
11import numpy as num
12import logging
14from pyrocko.guts import Float, Object
15from pyrocko.guts_array import Array
17guts_prefix = 'modelling'
19logger = logging.getLogger(__name__)
22class GriffithCrack(Object):
23 '''
24 Analytical Griffith crack model.
25 '''
27 width = Float.T(
28 help='Width equals to :math:`2 \\cdot a`.',
29 default=1.)
31 poisson = Float.T(
32 help='Poisson ratio :math:`\\nu`.',
33 default=.25)
35 shearmod = Float.T(
36 help='Shear modulus :math:`\\mu` [Pa].',
37 default=1.e9)
39 stressdrop = Array.T(
40 help='Stress drop array:'
41 '[:math:`\\sigma_{3,r} - \\sigma_{3,c}`, '
42 ':math:`\\sigma_{2,r} - \\sigma_{2,c}`, '
43 ':math:`\\sigma_{1,r} - \\sigma_{1,c}`] :math:`=` '
44 '[:math:`\\Delta \\sigma_{strike}, '
45 '\\Delta \\sigma_{dip}, \\Delta \\sigma_{tensile}`].',
46 default=num.array([0., 0., 0.]))
48 @property
49 def a(self):
50 '''
51 Half width of the crack in [m].
52 '''
54 return self.width / 2.
56 def disloc_infinite2d(self, x_obs):
57 '''
58 Calculation of dislocation at crack surface along x2 axis.
60 Follows equations by Pollard and Segall (1987) to calculate
61 dislocations for an infinite 2D crack extended in x3 direction,
62 opening in x1 direction and the crack extending in x2 direction.
64 :param x_obs:
65 Observation point coordinates along x2-axis.
66 If :math:`x_{obs} < -a` or :math:`x_{obs} > a`, output dislocations
67 are zero.
68 :type x_obs:
69 :py:class:`~numpy.ndarray`: ``(N,)``
71 :return:
72 Dislocations at each observation point in strike, dip and
73 tensile direction.
74 :rtype:
75 :py:class:`~numpy.ndarray`: ``(N, 3)``
76 '''
78 if type(x_obs) is not num.ndarray:
79 x_obs = num.array(x_obs)
81 factor = num.array([2. / self.shearmod])
82 factor = num.append(
83 factor, num.tile(
84 2. * (1. - self.poisson) / self.shearmod, (1, 2)))
85 factor[1] *= -1.
87 crack_el = (x_obs > -self.a) | (x_obs < self.a)
89 disl = num.zeros((x_obs.shape[0], 3))
90 disl[crack_el, :] = \
91 self.stressdrop * num.sqrt(
92 self.a**2 - num.tile(x_obs[crack_el, num.newaxis], (1, 3))**2) * \
93 factor
95 return disl
97 def disloc_circular(self, x_obs):
98 '''
99 Calculation of dislocation at crack surface along x2 axis.
101 Follows equations by Pollard and Segall (1987) to calculate
102 displacements for a circulat crack extended in x2 and x3 direction and
103 opening in x1 direction.
105 :param x_obs:
106 Observation point coordinates along axis through crack
107 centre. If :math:`x_{obs} < -a` or :math:`x_{obs} > a`, output
108 dislocations are zero.
109 :type x_obs:
110 :py:class:`~numpy.ndarray`: ``(N,)``
112 :return:
113 Dislocations at each observation point in strike, dip and
114 tensile direction.
115 :rtype:
116 :py:class:`~numpy.ndarray`: ``(N, 3)``
117 '''
119 if type(x_obs) is not num.ndarray:
120 x_obs = num.array(x_obs)
122 factor = num.array([4. / (self.shearmod * num.pi)])
123 factor = num.append(
124 factor, num.tile(
125 4. * (1. - self.poisson) / (self.shearmod * num.pi), (1, 2)))
126 factor[1] *= -1.
128 crack_el = (x_obs > -self.a) | (x_obs < self.a)
130 disl = num.zeros((x_obs.shape[0], 3))
131 disl[crack_el] = \
132 self.stressdrop * num.sqrt(
133 self.a**2 - num.tile(x_obs[crack_el, num.newaxis], (1, 3))**2) * \
134 factor
136 return disl
138 def _displ_infinite2d_along_x1(self, x1_obs):
139 '''
140 Calculation of displacement at crack surface along x1-axis.
142 Follows equations by Pollard and Segall (1987) to calculate
143 displacements for an infinite 2D crack extended in x3 direction,
144 opening in x1 direction and the crack tip in x2 direction.
146 :param x1_obs:
147 Observation point coordinates along x1-axis.
148 :type x1_obs:
149 :py:class:`~numpy.ndarray`: ``(N,)``
151 :return:
152 Displacements at each observation point in strike, dip and
153 tensile direction.
154 :rtype:
155 :py:class:`~numpy.ndarray`: ``(M, 3)``
156 '''
157 displ = num.zeros((x1_obs.shape[0], 3))
159 if self.stressdrop[0] != 0.:
160 sign = num.sign(x1_obs)
161 x1_ratio = x1_obs / self.a
163 displ[:, 0] = (
164 num.sqrt((x1_ratio)**2 + 1.) - num.abs(
165 x1_ratio)
166 ) * self.stressdrop[0] * self.a / self.shearmod * sign
168 return displ
170 def _displ_infinite2d_along_x2(self, x2_obs):
171 '''
172 Calculation of displacement at crack surface along x2-axis.
174 Follows equations by Pollard and Segall (1987) to calculate
175 displacements for an infinite 2D crack extended in x3 direction,
176 opening in x1 direction and the crack tip in x2 direction.
178 :param x2_obs:
179 Observation point coordinates along x2-axis.
180 :type x2_obs:
181 :py:class:`~numpy.ndarray`: ``(N,)``
183 :return:
184 Displacements at each observation point in strike, dip and
185 tensile direction.
186 :rtype:
187 :py:class:`~numpy.ndarray`: ``(N, 3)``
188 '''
190 crack_el = (x2_obs >= -self.a) & (x2_obs <= self.a)
192 displ = num.zeros((x2_obs.shape[0], 3))
194 if self.stressdrop[1] != 0.:
195 factor = (1. - 2. * self.poisson) / (2. * self.shearmod)
197 displ[crack_el, 2] = \
198 self.stressdrop[1] * factor * x2_obs[crack_el]
200 sign = num.sign(x2_obs)
202 displ[~crack_el, 2] = \
203 self.stressdrop[1] * factor * self.a * sign[~crack_el] * (
204 num.abs(x2_obs[~crack_el] / self.a) -
205 num.sqrt(x2_obs[~crack_el]**2 / self.a**2 - 1.))
207 return displ
209 def displ_infinite2d(self, x1_obs, x2_obs):
210 '''
211 Calculation of displacement at crack surface along different axis.
213 Follows equations by Pollard and Segall (1987) to calculate
214 displacements for an infinite 2D crack extended in x3 direction,
215 opening in x1 direction and the crack tip in x2 direction.
217 :param x1_obs:
218 Observation point coordinates along x1-axis.
219 If :math:`x1_obs = 0.`, displacment is calculated along x2-axis.
220 :type x1_obs:
221 :py:class:`~numpy.ndarray`: ``(M,)``
223 :param x2_obs:
224 Observation point coordinates along x2-axis.
225 If :math:`x2_obs = 0.`, displacment is calculated along x1-axis.
226 :type x2_obs:
227 :py:class:`~numpy.ndarray`: ``(N,)``
229 :return:
230 Displacements at each observation point in strike, dip and
231 tensile direction.
232 :rtype:
233 :py:class:`~numpy.ndarray`: ``(M, 3)`` or ``(N, 3)``
234 '''
236 if type(x1_obs) is not num.ndarray:
237 x1_obs = num.array(x1_obs)
238 if type(x2_obs) is not num.ndarray:
239 x2_obs = num.array(x2_obs)
241 if (x1_obs == 0.).all():
242 return self._displ_infinite2d_along_x2(x2_obs)
243 elif (x2_obs == 0.).all():
244 return self._displ_infinite2d_along_x1(x1_obs)
247__all__ = [
248 'GriffithCrack']