1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

# https://pyrocko.org - GPLv3 

# 

# The Pyrocko Developers, 21st Century 

# ---|P------/S----------~Lg---------- 

import numpy as num 

import logging 

 

from pyrocko.guts import Float, Object 

from pyrocko.guts_array import Array 

 

guts_prefix = 'modelling' 

 

logger = logging.getLogger('pyrocko.modelling.cracksol') 

 

 

class GriffithCrack(Object): 

''' 

Analytical Griffith crack model 

''' 

 

width = Float.T( 

help='Width equals to 2*a', 

default=1.) 

 

poisson = Float.T( 

help='Poisson ratio', 

default=.25) 

 

shearmod = Float.T( 

help='Shear modulus [Pa]', 

default=1.e9) 

 

stressdrop = Array.T( 

help='Stress drop array:' 

'[sigi3_r - sigi3_c, sigi2_r - sigi2_c, sigi1_r - sigi1_c]' 

'[dstress_Strike, dstress_Dip, dstress_Tensile]', 

default=num.array([0., 0., 0.])) 

 

@property 

def a(self): 

''' 

Half width of the crack 

''' 

 

return self.width / 2. 

 

@property 

def youngsmod(self): 

''' 

Shear (Youngs) modulus 

''' 

 

return 2. * self.shearmod * (1. + self.poisson) 

 

def disloc_infinite2d(self, x_obs): 

''' 

Calculation of dislocation at crack surface along x2 axis 

 

Follows equations by Pollard and Segall (1987) to calculate 

dislocations for an infinite 2D crack extended in x3 direction, 

opening in x1 direction and the crack extending in x2 direction. 

 

:param x_obs: Observation point coordinates along x2-axis. 

If :math:`x_{obs} < -a` or :math:`x_{obs} > a`, output dislocations 

are zero 

:type x_obs: :py:class:`numpy.ndarray`, ``(N,)`` 

 

:return: dislocations at each observation point in strike, dip and 

tensile direction. 

:rtype: :py:class:`numpy.ndarray`, ``(N, 3)`` 

''' 

 

if type(x_obs) is not num.ndarray: 

x_obs = num.array(x_obs) 

 

factor = num.array([2. / self.shearmod]) 

factor = num.append( 

factor, num.tile( 

2. * (1. - self.poisson) / self.shearmod, (1, 2))) 

factor[1] *= -1. 

 

crack_el = (x_obs > -self.a) | (x_obs < self.a) 

 

disl = num.zeros((x_obs.shape[0], 3)) 

disl[crack_el, :] = \ 

self.stressdrop * num.sqrt( 

self.a**2 - num.tile(x_obs[crack_el, num.newaxis], (1, 3))**2) * \ 

factor 

 

return disl 

 

def disloc_circular(self, x_obs): 

''' 

Calculation of dislocation at crack surface along x2 axis 

 

Follows equations by Pollard and Segall (1987) to calculate 

displacements for a circulat crack extended in x2 and x3 direction and 

opening in x1 direction. 

 

:param x_obs: Observation point coordinates along axis through crack 

centre. If :math:`x_{obs} < -a` or :math:`x_{obs} > a`, output 

dislocations are zero 

:type x_obs: :py:class:`numpy.ndarray`, ``(N,)`` 

 

:return: dislocations at each observation point in strike, dip and 

tensile direction. 

:rtype: :py:class:`numpy.ndarray`, ``(N, 3)`` 

''' 

 

if type(x_obs) is not num.ndarray: 

x_obs = num.array(x_obs) 

 

factor = num.array([4. / (self.shearmod * num.pi)]) 

factor = num.append( 

factor, num.tile( 

4. * (1. - self.poisson) / (self.shearmod * num.pi), (1, 2))) 

factor[1] *= -1. 

 

crack_el = (x_obs > -self.a) | (x_obs < self.a) 

 

disl = num.zeros((x_obs.shape[0], 3)) 

disl[crack_el] = \ 

self.stressdrop * num.sqrt( 

self.a**2 - num.tile(x_obs[crack_el, num.newaxis], (1, 3))**2) * \ 

factor 

 

return disl 

 

def _displ_infinite2d_along_x1(self, x1_obs): 

''' 

Calculation of displacement at crack surface along x1-axis 

 

Follows equations by Pollard and Segall (1987) to calculate 

displacements for an infinite 2D crack extended in x3 direction, 

opening in x1 direction and the crack tip in x2 direction. 

 

:param x1_obs: Observation point coordinates along x1-axis. 

:type x1_obs: :py:class:`numpy.ndarray`, ``(N,)`` 

 

:return: displacements at each observation point in strike, dip and 

tensile direction. 

:rtype: :py:class:`numpy.ndarray`, ``(M, 3)`` 

''' 

displ = num.zeros((x1_obs.shape[0], 3)) 

 

if self.stressdrop[0] != 0.: 

sign = num.sign(x1_obs) 

x1_ratio = x1_obs / self.a 

 

displ[:, 0] = ( 

num.sqrt((x1_ratio)**2 + 1.) - num.abs( 

x1_ratio) 

) * self.stressdrop[0] * self.a / self.shearmod * sign 

 

return displ 

 

def _displ_infinite2d_along_x2(self, x2_obs): 

''' 

Calculation of displacement at crack surface along x2-axis 

 

Follows equations by Pollard and Segall (1987) to calculate 

displacements for an infinite 2D crack extended in x3 direction, 

opening in x1 direction and the crack tip in x2 direction. 

 

:param x2_obs: Observation point coordinates along x2-axis. 

:type x2_obs: :py:class:`numpy.ndarray`, ``(N,)`` 

 

:return: displacements at each observation point in strike, dip and 

tensile direction. 

:rtype: :py:class:`numpy.ndarray`, ``(N, 3)`` 

''' 

 

crack_el = (x2_obs >= -self.a) & (x2_obs <= self.a) 

 

displ = num.zeros((x2_obs.shape[0], 3)) 

 

if self.stressdrop[1] != 0.: 

factor = (1. - 2. * self.poisson) / (2. * self.shearmod) 

 

displ[crack_el, 2] = \ 

self.stressdrop[1] * factor * x2_obs[crack_el] 

 

sign = num.sign(x2_obs) 

 

displ[~crack_el, 2] = \ 

self.stressdrop[1] * factor * self.a * sign[~crack_el] * ( 

num.abs(x2_obs[~crack_el] / self.a) - 

num.sqrt(x2_obs[~crack_el]**2 / self.a**2 - 1.)) 

 

return displ 

 

def displ_infinite2d(self, x1_obs, x2_obs): 

''' 

Calculation of displacement at crack surface along different axis 

 

Follows equations by Pollard and Segall (1987) to calculate 

displacements for an infinite 2D crack extended in x3 direction, 

opening in x1 direction and the crack tip in x2 direction. 

 

:param x1_obs: Observation point coordinates along x1-axis. 

If x1_obs = 0., displacment is calculated along x2-axis 

:type x1_obs: :py:class:`numpy.ndarray`, ``(M,)`` 

:param x2_obs: Observation point coordinates along x2-axis. 

If x2_obs = 0., displacment is calculated along x1-axis 

:type x2_obs: :py:class:`numpy.ndarray`, ``(N,)`` 

 

:return: displacements at each observation point in strike, dip and 

tensile direction. 

:rtype: :py:class:`numpy.ndarray`, ``(M, 3)`` or ``(N, 3)`` 

''' 

 

if type(x1_obs) is not num.ndarray: 

x1_obs = num.array(x1_obs) 

if type(x2_obs) is not num.ndarray: 

x2_obs = num.array(x2_obs) 

 

if (x1_obs == 0.).all(): 

return self._displ_infinite2d_along_x2(x2_obs) 

elif (x2_obs == 0.).all(): 

return self._displ_infinite2d_along_x1(x1_obs) 

 

 

__all__ = [ 

'GriffithCrack']