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
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: Observation point coordinates along x2-axis.
59 If :math:`x_{obs} < -a` or :math:`x_{obs} > a`, output dislocations
60 are zero
61 :type x_obs: :py:class:`numpy.ndarray`, ``(N,)``
63 :return: dislocations at each observation point in strike, dip and
64 tensile direction.
65 :rtype: :py:class:`numpy.ndarray`, ``(N, 3)``
66 '''
68 if type(x_obs) is not num.ndarray:
69 x_obs = num.array(x_obs)
71 factor = num.array([2. / self.shearmod])
72 factor = num.append(
73 factor, num.tile(
74 2. * (1. - self.poisson) / self.shearmod, (1, 2)))
75 factor[1] *= -1.
77 crack_el = (x_obs > -self.a) | (x_obs < self.a)
79 disl = num.zeros((x_obs.shape[0], 3))
80 disl[crack_el, :] = \
81 self.stressdrop * num.sqrt(
82 self.a**2 - num.tile(x_obs[crack_el, num.newaxis], (1, 3))**2) * \
83 factor
85 return disl
87 def disloc_circular(self, x_obs):
88 '''
89 Calculation of dislocation at crack surface along x2 axis
91 Follows equations by Pollard and Segall (1987) to calculate
92 displacements for a circulat crack extended in x2 and x3 direction and
93 opening in x1 direction.
95 :param x_obs: Observation point coordinates along axis through crack
96 centre. If :math:`x_{obs} < -a` or :math:`x_{obs} > a`, output
97 dislocations are zero
98 :type x_obs: :py:class:`numpy.ndarray`, ``(N,)``
100 :return: dislocations at each observation point in strike, dip and
101 tensile direction.
102 :rtype: :py:class:`numpy.ndarray`, ``(N, 3)``
103 '''
105 if type(x_obs) is not num.ndarray:
106 x_obs = num.array(x_obs)
108 factor = num.array([4. / (self.shearmod * num.pi)])
109 factor = num.append(
110 factor, num.tile(
111 4. * (1. - self.poisson) / (self.shearmod * num.pi), (1, 2)))
112 factor[1] *= -1.
114 crack_el = (x_obs > -self.a) | (x_obs < self.a)
116 disl = num.zeros((x_obs.shape[0], 3))
117 disl[crack_el] = \
118 self.stressdrop * num.sqrt(
119 self.a**2 - num.tile(x_obs[crack_el, num.newaxis], (1, 3))**2) * \
120 factor
122 return disl
124 def _displ_infinite2d_along_x1(self, x1_obs):
125 '''
126 Calculation of displacement at crack surface along x1-axis
128 Follows equations by Pollard and Segall (1987) to calculate
129 displacements for an infinite 2D crack extended in x3 direction,
130 opening in x1 direction and the crack tip in x2 direction.
132 :param x1_obs: Observation point coordinates along x1-axis.
133 :type x1_obs: :py:class:`numpy.ndarray`, ``(N,)``
135 :return: displacements at each observation point in strike, dip and
136 tensile direction.
137 :rtype: :py:class:`numpy.ndarray`, ``(M, 3)``
138 '''
139 displ = num.zeros((x1_obs.shape[0], 3))
141 if self.stressdrop[0] != 0.:
142 sign = num.sign(x1_obs)
143 x1_ratio = x1_obs / self.a
145 displ[:, 0] = (
146 num.sqrt((x1_ratio)**2 + 1.) - num.abs(
147 x1_ratio)
148 ) * self.stressdrop[0] * self.a / self.shearmod * sign
150 return displ
152 def _displ_infinite2d_along_x2(self, x2_obs):
153 '''
154 Calculation of displacement at crack surface along x2-axis
156 Follows equations by Pollard and Segall (1987) to calculate
157 displacements for an infinite 2D crack extended in x3 direction,
158 opening in x1 direction and the crack tip in x2 direction.
160 :param x2_obs: Observation point coordinates along x2-axis.
161 :type x2_obs: :py:class:`numpy.ndarray`, ``(N,)``
163 :return: displacements at each observation point in strike, dip and
164 tensile direction.
165 :rtype: :py:class:`numpy.ndarray`, ``(N, 3)``
166 '''
168 crack_el = (x2_obs >= -self.a) & (x2_obs <= self.a)
170 displ = num.zeros((x2_obs.shape[0], 3))
172 if self.stressdrop[1] != 0.:
173 factor = (1. - 2. * self.poisson) / (2. * self.shearmod)
175 displ[crack_el, 2] = \
176 self.stressdrop[1] * factor * x2_obs[crack_el]
178 sign = num.sign(x2_obs)
180 displ[~crack_el, 2] = \
181 self.stressdrop[1] * factor * self.a * sign[~crack_el] * (
182 num.abs(x2_obs[~crack_el] / self.a) -
183 num.sqrt(x2_obs[~crack_el]**2 / self.a**2 - 1.))
185 return displ
187 def displ_infinite2d(self, x1_obs, x2_obs):
188 '''
189 Calculation of displacement at crack surface along different axis
191 Follows equations by Pollard and Segall (1987) to calculate
192 displacements for an infinite 2D crack extended in x3 direction,
193 opening in x1 direction and the crack tip in x2 direction.
195 :param x1_obs: Observation point coordinates along x1-axis.
196 If :math:`x1_obs = 0.`, displacment is calculated along x2-axis
197 :type x1_obs: :py:class:`numpy.ndarray`, ``(M,)``
198 :param x2_obs: Observation point coordinates along x2-axis.
199 If :math:`x2_obs = 0.`, displacment is calculated along x1-axis
200 :type x2_obs: :py:class:`numpy.ndarray`, ``(N,)``
202 :return: displacements at each observation point in strike, dip and
203 tensile direction.
204 :rtype: :py:class:`numpy.ndarray`, ``(M, 3)`` or ``(N, 3)``
205 '''
207 if type(x1_obs) is not num.ndarray:
208 x1_obs = num.array(x1_obs)
209 if type(x2_obs) is not num.ndarray:
210 x2_obs = num.array(x2_obs)
212 if (x1_obs == 0.).all():
213 return self._displ_infinite2d_along_x2(x2_obs)
214 elif (x2_obs == 0.).all():
215 return self._displ_infinite2d_along_x1(x1_obs)
218__all__ = [
219 'GriffithCrack']