1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5import numpy as num
6import logging
8from pyrocko.guts import Float, Object
9from pyrocko.guts_array import Array
11guts_prefix = 'modelling'
13logger = logging.getLogger(__name__)
16class GriffithCrack(Object):
17 '''
18 Analytical Griffith crack model.
19 '''
21 width = Float.T(
22 help='Width equals to :math:`2 \\cdot a`.',
23 default=1.)
25 poisson = Float.T(
26 help='Poisson ratio :math:`\\nu`.',
27 default=.25)
29 shearmod = Float.T(
30 help='Shear modulus :math:`\\mu` [Pa].',
31 default=1.e9)
33 stressdrop = Array.T(
34 help='Stress drop array:'
35 '[:math:`\\sigma_{3,r} - \\sigma_{3,c}`, '
36 ':math:`\\sigma_{2,r} - \\sigma_{2,c}`, '
37 ':math:`\\sigma_{1,r} - \\sigma_{1,c}`] :math:`=` '
38 '[:math:`\\Delta stress_{strike}, '
39 '\\Delta stress_{dip}, \\Delta stress_{tensile}`].',
40 default=num.array([0., 0., 0.]))
42 @property
43 def a(self):
44 '''
45 Half width of the crack in [m].
46 '''
48 return self.width / 2.
50 def disloc_infinite2d(self, x_obs):
51 '''
52 Calculation of dislocation at crack surface along x2 axis.
54 Follows equations by Pollard and Segall (1987) to calculate
55 dislocations for an infinite 2D crack extended in x3 direction,
56 opening in x1 direction and the crack extending in x2 direction.
58 :param x_obs:
59 Observation point coordinates along x2-axis.
60 If :math:`x_{obs} < -a` or :math:`x_{obs} > a`, output dislocations
61 are zero.
62 :type x_obs:
63 :py:class:`~numpy.ndarray`: ``(N,)``
65 :return:
66 Dislocations at each observation point in strike, dip and
67 tensile direction.
68 :rtype:
69 :py:class:`@numpy.ndarray`: ``(N, 3)``
70 '''
72 if type(x_obs) is not num.ndarray:
73 x_obs = num.array(x_obs)
75 factor = num.array([2. / self.shearmod])
76 factor = num.append(
77 factor, num.tile(
78 2. * (1. - self.poisson) / self.shearmod, (1, 2)))
79 factor[1] *= -1.
81 crack_el = (x_obs > -self.a) | (x_obs < self.a)
83 disl = num.zeros((x_obs.shape[0], 3))
84 disl[crack_el, :] = \
85 self.stressdrop * num.sqrt(
86 self.a**2 - num.tile(x_obs[crack_el, num.newaxis], (1, 3))**2) * \
87 factor
89 return disl
91 def disloc_circular(self, x_obs):
92 '''
93 Calculation of dislocation at crack surface along x2 axis.
95 Follows equations by Pollard and Segall (1987) to calculate
96 displacements for a circulat crack extended in x2 and x3 direction and
97 opening in x1 direction.
99 :param x_obs:
100 Observation point coordinates along axis through crack
101 centre. If :math:`x_{obs} < -a` or :math:`x_{obs} > a`, output
102 dislocations are zero.
103 :type x_obs:
104 :py:class:`~numpy.ndarray`: ``(N,)``
106 :return:
107 Dislocations at each observation point in strike, dip and
108 tensile direction.
109 :rtype:
110 :py:class:`~numpy.ndarray`: ``(N, 3)``
111 '''
113 if type(x_obs) is not num.ndarray:
114 x_obs = num.array(x_obs)
116 factor = num.array([4. / (self.shearmod * num.pi)])
117 factor = num.append(
118 factor, num.tile(
119 4. * (1. - self.poisson) / (self.shearmod * num.pi), (1, 2)))
120 factor[1] *= -1.
122 crack_el = (x_obs > -self.a) | (x_obs < self.a)
124 disl = num.zeros((x_obs.shape[0], 3))
125 disl[crack_el] = \
126 self.stressdrop * num.sqrt(
127 self.a**2 - num.tile(x_obs[crack_el, num.newaxis], (1, 3))**2) * \
128 factor
130 return disl
132 def _displ_infinite2d_along_x1(self, x1_obs):
133 '''
134 Calculation of displacement at crack surface along x1-axis.
136 Follows equations by Pollard and Segall (1987) to calculate
137 displacements for an infinite 2D crack extended in x3 direction,
138 opening in x1 direction and the crack tip in x2 direction.
140 :param x1_obs:
141 Observation point coordinates along x1-axis.
142 :type x1_obs:
143 :py:class:`~numpy.ndarray`: ``(N,)``
145 :return:
146 Displacements at each observation point in strike, dip and
147 tensile direction.
148 :rtype:
149 :py:class:`~numpy.ndarray`: ``(M, 3)``
150 '''
151 displ = num.zeros((x1_obs.shape[0], 3))
153 if self.stressdrop[0] != 0.:
154 sign = num.sign(x1_obs)
155 x1_ratio = x1_obs / self.a
157 displ[:, 0] = (
158 num.sqrt((x1_ratio)**2 + 1.) - num.abs(
159 x1_ratio)
160 ) * self.stressdrop[0] * self.a / self.shearmod * sign
162 return displ
164 def _displ_infinite2d_along_x2(self, x2_obs):
165 '''
166 Calculation of displacement at crack surface along x2-axis.
168 Follows equations by Pollard and Segall (1987) to calculate
169 displacements for an infinite 2D crack extended in x3 direction,
170 opening in x1 direction and the crack tip in x2 direction.
172 :param x2_obs:
173 Observation point coordinates along x2-axis.
174 :type x2_obs:
175 :py:class:`~numpy.ndarray`: ``(N,)``
177 :return:
178 Displacements at each observation point in strike, dip and
179 tensile direction.
180 :rtype:
181 :py:class:`~numpy.ndarray`: ``(N, 3)``
182 '''
184 crack_el = (x2_obs >= -self.a) & (x2_obs <= self.a)
186 displ = num.zeros((x2_obs.shape[0], 3))
188 if self.stressdrop[1] != 0.:
189 factor = (1. - 2. * self.poisson) / (2. * self.shearmod)
191 displ[crack_el, 2] = \
192 self.stressdrop[1] * factor * x2_obs[crack_el]
194 sign = num.sign(x2_obs)
196 displ[~crack_el, 2] = \
197 self.stressdrop[1] * factor * self.a * sign[~crack_el] * (
198 num.abs(x2_obs[~crack_el] / self.a) -
199 num.sqrt(x2_obs[~crack_el]**2 / self.a**2 - 1.))
201 return displ
203 def displ_infinite2d(self, x1_obs, x2_obs):
204 '''
205 Calculation of displacement at crack surface along different axis.
207 Follows equations by Pollard and Segall (1987) to calculate
208 displacements for an infinite 2D crack extended in x3 direction,
209 opening in x1 direction and the crack tip in x2 direction.
211 :param x1_obs:
212 Observation point coordinates along x1-axis.
213 If :math:`x1_obs = 0.`, displacment is calculated along x2-axis.
214 :type x1_obs:
215 :py:class:`~numpy.ndarray`: ``(M,)``
216 :param x2_obs:
217 Observation point coordinates along x2-axis.
218 If :math:`x2_obs = 0.`, displacment is calculated along x1-axis.
219 :type x2_obs:
220 :py:class:`~numpy.ndarray`: ``(N,)``
222 :return:
223 Displacements at each observation point in strike, dip and
224 tensile direction.
225 :rtype:
226 :py:class:`~numpy.ndarray`: ``(M, 3)`` or ``(N, 3)``
227 '''
229 if type(x1_obs) is not num.ndarray:
230 x1_obs = num.array(x1_obs)
231 if type(x2_obs) is not num.ndarray:
232 x2_obs = num.array(x2_obs)
234 if (x1_obs == 0.).all():
235 return self._displ_infinite2d_along_x2(x2_obs)
236 elif (x2_obs == 0.).all():
237 return self._displ_infinite2d_along_x1(x1_obs)
240__all__ = [
241 'GriffithCrack']