Employment of Fracture Mechanics Criteria for Accurate Assessment of the Full Set of Elastic Constants of Orthorhombic/Tetragonal Mono-Crystalline YBCO

The effect of elastic constants, cij, on the nature (easy or difficult) of a cleavage system in mono-crystalline YBa2Cu3O7−δ is investigated by employing a novel three-dimensional eigenfunction expansion technique, based in part on the separation of the thickness variable and partly on a modified Frobenius-type series expansion technique in conjunction with Eshelby–Stroh formalism. Out of the three available, complete sets of elastic constants, only the experimental measurements using resonant ultrasound spectroscopy merit serious attention, despite reported values of c12 and, to a lesser extent, c66 being excessively high. The present investigation considers six through-thickness crack systems weakening orthorhombic mono-crystalline Yttrium barium copper oxide (YBCO) plates. More importantly, the present investigation establishes sufficient conditions for crack path stability/instability, which entail a cleavage system being easy or difficult, i.e., whether a crack would propagate in its original plane/direction or deflect to a different one. This criterion of fracture mechanics is then employed for accurate determination of the full set of elastic constants of superconducting mono-crystalline YBCO. Finally, heretofore unavailable results pertaining to the through-thickness variations of stress intensity factors and energy release rates for a crack corresponding to symmetric and skew-symmetric hyperbolic cosine loads, which also satisfy the boundary conditions on the plate surfaces, bridge a longstanding gap.


Introduction
The elastic constants of engineering materials are crucial for understanding the deformation and failure behaviors of structural components. From a microscopic perspective, their importance arises from their intimate relationships to such solid-state phenomena as specific heat, Debye temperature (Θ D ), and the Grünelsen parameter [1].
Modern applications of mono-crystalline high-Tc (critical temperature) superconductors (HTS), discovered by Bednorz and Muller [2], include Josephson junctions, which can act as a switch for magnetic fields, or alternatively, perform the function of a magnetic device, such as superconducting quantum interference devices (SQUID) [3]. The Θ D of a superconductor such as mono-crystalline YBCO (YBa 2 Cu 3 O 7−δ ) can be determined from the knowledge of the elastic constants, c ij , in a manner described by Equations (2) and (3) of Lei et al. [1], which, in combination with the electron-phonon coupling parameter, λ * can be used to compute the superconducting transition temperature, T c [4,5].
The above review of the literature reveals an absence of reliable and accurate experimentally measured complete sets of the nine elastic constants needed for characterization of the deformation/fracture, as well as other solid-state (e.g., Θ D , T c , etc.) behaviors of superconducting (orthorhombic) YBCO single crystals. This calls for a reliable criterion for assessment of the measured data that would allow us to come up with a reasonably accurate complete set of nine elastic constants. This is the primary objective of the present investigation. One effective way to address this important issue is to analytically examine the effects of elastic constants on crack path stability/instability in mono-crystalline YBCO and compare them with the experimental results for easy cleavage planes, reported by Cook et al. [6], Raynes [9], and Granozio and di Uccio [14], among others. In what follows, the above-mentioned modified eigenfunction expansion technique, based in part on the separation of the thickness variable and partly on the Eshelby-Stroh type affine transformation, is developed to derive a three-dimensional asymptotic stress field in the vicinity of the front of a semi-infinite through-thickness crack weakening an infinite orthorhombic mono-crystalline plate of finite thickness and subjected to far-field mode I/II loadings. Crack-face boundary conditions and those that are prescribed on the top and bottom (free or fixed) surfaces of the plate are exactly satisfied. The present investigation considers six through-crack systems-(010) [ length direction-weakening orthorhombic mono-crystalline plates. Explicit expressions for the singular stresses in the vicinity of the front of a through-thickness crack weakening an orthorhombic mono-crystalline plate, subjected to far-field mode I/II loadings, are presented. In addition, through-thickness distribution of the stress intensity factors, and energy release rates are also presented.
Next, the important issue of easy or difficult cleavage plane and the related question of crack deflection criterion is discussed. The latter is based on the relative fracture energy (or the energy release rate) available for possible "fracture paths" [18]. This said, it is noteworthy that the Griffith energy balance-based criterion cannot be regarded as a sufficient condition [80]. This is because Griffith's criterion is "Not really a fracture criterion but only a necessary condition for fracture" [80]. This calls for establishment of a sufficient condition for determination of an easy or difficult cleavage plane, and the related question of the crack path stability/instability criterion. This is the second and somewhat more important objective of the present study.
One Holy Grail issue in fracture mechanics of anisotropic solids is to find a dimensionless parameter akin to Reynold's number in fluid flow problems, crossing a critical value which signifies transition from one regime to another, such as the critical value of Reynold's number, above which the flow is turbulent and below which it is laminar. It is an attendant issue relating to crack deflection that remains heretofore unaddressed, which is the third objective of the present investigation. In a similar vein, just as the introduction of Reynold's number facilitated the design and setting up of experiments in addition to experimental verification of analytical and computational solutions in fluid dynamics, a similar expectation-cum-need exists in the important field of fracture mechanics of anisotropic solids, which is the final objective of the present study. Here, the accuracy and efficacy of the available experimental results on the elastic constants of YBCO single crystals, measured by modern experimental techniques with resolutions at the atomic scale (or nearly so), such as X-ray diffraction [65], the ultrasound technique [66,70,71], neutron diffraction [73]/scattering [67], Brillouin spectroscopy [68,69] /scattering [72], resonant ultrasound spectroscopy [1,74,75], and the like, is assessed with a powerful theoretical analysis on crack path stability/instability, in part based on a dimensionless parameter, such as the planar anisotropic ratio.
The present study, although to a smaller extent a review of earlier work on this topic, is largely based on original research on this subject. The topic, which covers mathematics (e.g., branch point/branch cut, asymptotic, solution to 3D mixed boundary-value problem, and necessary and sufficient condition for fracture), solid-state physics/chemistry (e.g., stoichiometry, superconductivity, crystal structures, and single crystal cleavage) and engineering (e.g., 3D fracture mechanics), has, so far, remained largely unexplored in the Literature.

Formulation of the Problem
One of the most important cleavage systems for mono-crystalline orthorhombic YBCO is (001) [100] × [010] ({crack plane}<crack front>x<initial propagation direction>). In what follows, the deformation behavior in the vicinity of the front of a semi-infinite through-thickness crack, weakening an infinite orthorhombic YBCO plate of thickness, 2 h (Figures 1 and 2) is analyzed in detail. Here, the z-axis is placed along the straight crack front, [100], while the coordinates x [010], y [001] are used to define the directions along the length of the crack (propagation direction) and the direction transverse to it, respectively, in the middle plane of the single crystal plate. u, v and w represent the components of the displacements in x [010], y [001] and z [100] directions, respectively. The corresponding stress-strain relations are given as follows: where c ij , i, j = 1,...,6, denotes the elastic stiffness constants of an orthorhombic monocrystalline plate. σ x , σ y , σ z represent the normal stresses, and τ x y , τ x z , τ y z denote the shear stresses, while ε x , ε y , ε z denote normal strains, and γ x y γ x z , γ y z represent the shear strains.    The assumed displacement functions for the three-dimensional crack problem under

Singular Stress Fields in the Vicinity of a Crack Front Weakening an Orthotropic/Orthorhombic Lamina/Single Crystal under General Loading
The assumed displacement functions for the three-dimensional crack problem under consideration are selected on the basis of the separation of z variables. These are as given below [22,56-58,64]: u( x, y, z) = e i k z U( x, y), (9) v( x, y, z) = e i k z V( x, y), (10) w( x, y, z) = e i k z W( x, y). (11) It may be noted that since the z-dependent term and its first partial derivative can either be bounded and integrable at most admitting ordinary discontinuities, or the first partial derivative at worst be square integrable (in the sense of Lebesgue integration) in its interval z∈ [−h, h], i.e., admitting singularities weaker than square root (i.e., z (−1/2+ε) , ε), it can be best represented by a Fourier series [22,25,58]. The latter case is justified by Parseval's theorem [81], and its physical implication is that of satisfying the criterion of finiteness of local strain energy and path independence [82]. Substitution of Equations (9)- (11) into Equations (2)-(5) yields the following system of coupled partial differential equations (PDE's): where Appl. Mech. 2023, 4

592
The solution to the system of coupled partial differential Equations (12)- (14) subjected to the most general loading can now be sought in the form of the following modified Frobenius type series in terms of the variable x 1 + p y 1 as follows: Out of the various combinations, such as (a , b , c ), (a, b, c), (a , b, c), (a, b , c), (a, b, c ), (a , b , c), (a , b, c ), and (a, b , c ), only the first two groupings can produce meaningful solutions, for the mode I/II and mode III loading cases, respectively. This step permits the separation of mode III from modes I/II. The first grouping is described below, while the second one has already been employed for the anti-plane shear case [59][60][61].

Singular Stress Fields in the Vicinity of a (001)[100] Through-Crack Front Propagating under Mode I (Extension/Bending) and Mode II (Sliding Shear/Twisting) in [010] Direction
While the separability of the x 1 and y 1 variables is the hallmark of the derivation of the singular stress fields weakening an isotropic plate 55] subjected to the far-field mode I/II loading, the same cannot be true for its orthorhombic counterpart, wherein the solution to the system of coupled partial differential Equation (5) must be sought in the form of a modified Frobenius type series in terms of a combined variable x 1 + p y 1 , which is a mathematical representation of an affine transformation employed earlier by Eshelby et al. [62], Stroh [16], and Shih et al. [17], as shown below [22,53,54,60]. It may be remarked here that the solution methodology offered by these authors, especially Refs. [17,62], are entirely different from what follows.
On substitution of Equations (20)- (22) into Equations (12)- (14), and equating the coefficients of ( x 1 + p y 1 ) s+2n−2 , a set of recurrent relationships can be derived, which, for n = 0, results in the following: which, in turn, yields, for the non-trivial case, the characteristic equation for the coupled partial differential Equations (3)-(5) or (12)-(14) as given below: where the normalized elastic parameter, κ = 1/ χ, is given by in which E 3 is z-direction Young's modulus, G 23 is the shear modulus in the y-z plane, while ν 12 is the major Poisson's ratio in the x-y plane. ν 32 denotes the minor Poisson's ratio in the y-z plane, while ν 13 and ν 31 represent the major and minor Poisson's ratios, respectively, in the x-z plane. It also is sometimes convenient to relate κ = 1/ χ, as will be seen later, to A = 1/ λ, with A being the planar anisotropic ratio (in the x [010]-y [001] plane), as shown below.
in which It can easily be seen from Equation (27) that A is higher when the shear stiffness (modulus) and major Poisson's ratio in the x [010], y [001] plane assume larger magnitudes. This simple fact assumes great importance as this investigation aims to solve one Holy Grail issue in the fracture mechanics of anisotropic media, of coming up with a dimensionless parameter akin to Reynold's number in fluid flow problems, crossing a critical value which signifies transition from one regime to another, such as the critical value of Reynold's number above which the flow is turbulent and below which it is laminar. It is an attendant issue relating to crack deflection in mono-crystalline orthorhombic YBCO.

Case (a): Complex Roots
where ξ = 1 √ 2 valid for A > 1 or, equivalently, κ > √ c 33 /c 22 . It may be noted here that a crack front and the corresponding semi-infinite crack represent a branch point and a branch cut, respectively. Therefore, the general asymptotic form for the displacement and stress fields in the vicinity of the crack front can most conveniently be obtained by employing a set of two polar coordinate systems, ( ρ, ψ) and ( ρ , ψ ), derived from the aforementioned affine transformation and expressed in terms of the cylindrical polar coordinate system ( r, θ, z), as follows: ρ cos ( ψ) = r cos( θ) + ξ sin( θ) , ρ sin ( ψ) = r η sin( θ) , (34) ρ cos ( ψ ) = r cos( θ) − ξ sin( θ) , ρ sin ( ψ ) = r − η sin( θ) , (35) where ρ = r cos( θ) + ξ sin( θ) 2 + η 2 sin 2 ( θ) and cos ψ( θ) = cos( θ)+ ξ sin( θ) (cos( θ)+ ξ sin( θ)) 2 + η 2 sin 2 ( θ) 1/2 , (39) cos The components of the displacement vector and the stress tensor in the immediate neighborhood of the crack from can now be written as shown below: and where and It may be noted that since s or Re s (when s is complex) is positive, all the higher order terms in Equations (45)-(48) vanish as r → 0. The non-vanishing components of asymptotic displacement vector and stress tensor in the cylindrical polar coordinate system ( r, θ, z) can now be obtained by using the standard vector and the second-rank tensor transformation rule: The stress component, σ z , is as given in Equation (48). The substitution of Equation (55), in conjunction with Equations (45)- (48), into the stress-free condition on the crack-side surfaces, given by Equation (7), gives rise to four homogeneous equations, which finally yields: Either or sin( s − 1)π = 0.

Symmetric (Mode I) Loading (Extension/Bending)
For the case of far-field symmetric (Mode I) loading, the following boundary conditions can be applied: when s = 1/2, substitution of Equation (55) in conjunction with Equations (45)-(47) into Equations (58) and (59), yields the following: Finally, on substitution of Equations (60)-(62) into Equations (42), (43) and (45)- (48), the relevant components of the displacement vector and the stress tensor in the immediate neighborhood of a semi-infinite crack front, under Mode I loading, can be written as given below: where the mode I stress intensity factor, K I ( z), is defined as follows: For the case of far-field skew-symmetric (Mode II) loading, the following boundary conditions can be applied: When s = 1/2, substitution of Equation (55) in conjunction with Equations (45)-(47) into Equations (70) and (71) yields the following: Finally, on substitution of Equations (72)- (74) into Equations (42), (43) and (45) in which the mode II stress intensity factor, K I I ( z), is defined as follows: σ z ( r, θ, z) is given by Equation (69). A critical examination of Equations (66) and (79) reveals the following interesting relationship: Since K I ( z) and K I I ( z) are directly proportional to far-field loadings, σ ∞ y ( z) and τ ∞ x y ( z), respectively, this interesting result implies that both these loadings produce identical (except the negative sign for τ x y ) θ-dependence of these two near-field driving stresses, σ y ( r, θ, z) Mode I / τ x y ( r, θ, z) Mode I I : This is unlike what happens in the case of an isotropic material, where no such θinvariance (i.e., identical θ-dependence) exists, as the following relationship would clearly demonstrate [25,55]: Comparison of Equations (81)- (84) shows that while the ratio of mode I to mode II stress intensity factors or its far-field loading counterpart is equal to negative times the ratio of the corresponding driving stresses, σ y ( r, θ, z) Mode I / τ x y ( r, θ, z) Mode I I , for the case of an orthorhombic crystal with complex roots (valid for A > 1 or, equivalently, κ > √ c 33 /c 22 ), and this relationship is invariant with respect to coordinate transformation (here specifically with respect to θ variation), the same cannot be said about the corresponding relationship for an isotropic material. The coupling between cos ψ/2 and sin ψ/2 , and similar coupling between cos ψ /2 and sin ψ /2 in the two driving stresses, σ y ( r, θ, z) Mode I and τ x y ( r, θ, z) Mode I I , which has also been reported in earlier studies [22,56], is instrumental in causing the above remarkable relationship given by Equations (81) and (82), which signifies mode mixity. It is also somewhat counter-intuitive that an invariant (more specifically, with respect to θ) relationship, such as Equations (81) and (82), resulting from complex roots (valid for A > 1 or, equivalently, κ > √ c 33 /c 22 ), guarantees the fact of the cleavage system under consideration being difficult (i.e., not easy), thus further concomitant in ensuring crack deflection or turning from this difficult cleavage system onto a nearby available easy one, violating the self-similarity of crack growth or propagation in the process. In contrast, again counter-intuitively, a non-invariant (more specifically, with respect to θ) relationship, such as Equations (83) and (84) resulting from the degenerate case ( A = 1 or, equivalently, κ = 1), guarantees the cleavage system under consideration of being supereasy, thus ensuring the self-similarity (i.e., θ-invariance) of crack growth or propagation in the process.
The dimensionless parameters, such as the anisotropic ratio, A, or, equivalently, the normalized elastic parameter, κ, can serve as the Holy Grail quantity for an a priori determination of the status of a cleavage system to be easy or difficult, very much akin to Reynold's number for fluid flow problems, crossing a critical value of which signifies transition from one regime to another. Here, the anisotropic ratio, A, or, equivalently, normalized elastic parameter, κ, for a (001) [100] × [010] cleavage system, crossing the critical value of 1 or √ c 33 /c 22 , respectively, signifies transitioning from self-similar crack growth or propagation to crack deflection or turning from a difficult cleavage system onto a nearby easy one.
As has been suggested by Sedov [83], similarity analysis is an effective tool to solve complex problems in mechanics. This can be employed to establish a sorely needed sufficient condition for the problem at hand and is further elaborated in Section 7.2 below.

Case (b): Imaginary Roots
Equation (24) can also have four imaginary roots given by where is valid for which corresponds to a candidate plane of minimum surface energy. As has been explained earlier, the general asymptotic form for the displacement and stress fields in the vicinity of the crack front (a branch point in 2D form) can most conveniently be obtained by employing a set of two polar coordinate systems, ( ρ, ψ) and ( ρ , ψ ), derived from the affine transformation discussed earlier, and expressed in terms of the cylindrical polar coordinate system ( r, θ, z), as follows: in which and cos ψ 1 ( θ) = cos( θ) The components of the displacement vector and the stress tensor in the immediate neighborhood of the crack front can now be written as shown below: and where and D b ( z) is same as given earlier in Equation (53). As before, since s or Re s (when s is complex) is positive, all the higher order terms in Equations (101)-(104) vanish as r→ 0. The non-vanishing components of the asymptotic displacement vector and the stress tensor in the cylindrical polar coordinate system ( r, θ, z), can now be obtained by using the standard vector and the second-rank tensor transformation rule, given earlier by Equations (54) and (55). The stress component, σz, is as given in Equation (104).
Substitution of Equation (55) in conjunction with Equations (101)-(103) into the stressfree condition on the crack-side surfaces, given by Equation (7), gives rise to four homogeneous equations, which finally yields or Equation (109) yields the lowest nonvanishing eigenvalue, s = 1/2, 0 < s < 1, thus meeting the criterion of locally finite energy. Equation (110), in contrast, gives rise to s = 0, 1, which can take care of rigid body translation and rotation, respectively. Interestingly, s = 1 also accounts for the T-stress.

Symmetric (Mode I) Loading (Extension/Bending)
For s = 1/2, and Finally, on substitution of Equations (111) and (112) into Equations (98), (99) and (101)-(104), the relevant components of the displacement vector and the stress tensor in the immediate neighborhood of a semi-infinite crack front, under Mode I far-field loading, can be written as given below: where the mode I stress intensity factor, K I ( z), is defined as follows: σ z ( r, θ, z) is given by Equation (69).

Skew-Symmetric (Mode II) Loading (Sliding Shear/Twisting)
For s = 1/2, and Finally, on substitution of Equations (119) and (120) into Equations (98), (99) and (101)-(104), the relevant components of the displacement vector and the stress tensor in the immediate neighborhood of a semi-infinite crack front, under Mode II loading, can be written as given below: in which the mode II stress intensity factor, K I I ( z), is defined as follows: is given by Equation (69). A critical examination of Equations (116) and (126) reveals the following relationship involving the two near-field crack driving stresses, and τ x y ( r, θ, z) Mode I I , and corresponding far-field applied stresses, σ ∞ y ( z) and τ ∞ x y ( z), proportional to K I ( z) and K I I ( z): Equation (127) shows an absence of θ-invariance in the relationship involving the two far-field and two near-field stresses mentioned above in the case of an orthorhombic crystal with imaginary roots (valid for A < 1 or, equivalently, κ < √ c 33 /c 22 ) as a result of the absence of coupling between cos ψ 1 /2 and sin ψ 1 /2 (and similar coupling between cos ψ 1 /2 and sin ψ 1 /2 ) in the two driving stresses, σ y ( r, θ, z) Mode I and τ x y ( r, θ, z) Mode I I . The lack of coupling between cos ψ 1 /2 and sin ψ 1 /2 (and similar coupling between cos ψ 1 /2 and sin ψ 1 /2 ) in the two driving stresses, σ y ( r, θ, z) Mode I and τ x y ( r, θ, z) Mode I I , which has also been reported in earlier studies [22,56], is instrumental in guaranteeing the fact of the cleavage system under consideration being easy, thus ensuring self-similar crack growth or propagation in the process, in a manner similar to what happens in an isotropic solid. Comparison of Equations (83), (84) and (127) also shows a measure of similarity in terms of the θ-dependence discussed above; see further discussion in Section 7.2 below.

Satisfaction of Traction-Free Boundary Conditions
The stress field in the vicinity of the front of a semi-infinite crack under inplane extension and out-of-plane bending, respectively, can be recovered if in Equations (68), (118) or (80), (126), even and odd functions are selected from D b ( z * ): where z * = z/h. By using the boundary condition on the free plate surfaces, z * = z/h = ±1, the general form of D bs ( z * ) can be obtained as Hence, K I ( z * ) = K Is ( z * ) and K I I ( z * ) = K I Is ( z * ) would represent symmetric stress intensity factors; see Section 8 and the penultimate figure below.
D ba ( z * ) that satisfies the traction-free condition on the plate surfaces is, in the absence of discontinuity of the function at z * = z/h = 0, generally given by In the presence of discontinuity of the function at z * = z/h = 0, D ba ( z * ) must be written as follows: As a consequence, K I ( z * ) = K Ia ( z * ) and K I I ( z * ) = K I Ia ( z * ) would represent antisymmetric stress intensity factors; see Section 8 below. If the boundary conditions on the traction-free plate surfaces are satisfied, all the displacements and singular stresses vanish on the plate surfaces in the vicinity of the crack front.

Hyperbolic Cosine Distributed Far-Field Loading
Hyperbolic cosine distributed far-field loading, which is proportional to cosh( z * ), | z * | < 1, is applied. The applied symmetric loading function and the corresponding "stress intensity factors" (valid for | z * | ≤ 1) are proportional to The corresponding Fourier series can be derived as follows: The applied antisymmetric loading function (valid for | z * | < 1) and the corresponding "stress intensity factors" (valid for | z * | ≤ 1) are proportional to Since D ba ( z * ) has a discontinuity at z * = 0, the corresponding Fourier series can be obtained as given below:

Through-Thickness Distribution of Stress Intensity Factors (Modes I and II)
The stress intensity factors, K I ( z * ) and K I I ( z), cannot be determined unless the farfield loading and a characteristic length (e.g., crack geometry) are specified. Sih et al. [17] have shown the applicability of the complex variable approach in conjunction with the eigenfunction expansion approach in the derivation of the two-dimensional stress intensity factors for anisotropic plates. The stress intensity factor for an infinite orthorhombic/tetragonal mono-crystalline plate with a central crack of length, 2a, and subjected to far-field mode I/II loading is available for the two-dimensional case [17], and can easily be extended to the present three-dimensional case as follows: where with C 1 being a constant. On substitution of Equations (138)-(149) into it, Equation (137) can be expressed in cylindrical polar coordinates as follows: which finally gives for both complex and imaginary roots. Equations (142) and (143) reduce to their twodimensional counterparts [53], by taking D b ( z * ) = 1. It may further be noted that the normalization factor, K i ( z)/ K i,2D , i = I, II, is equal to D b ( z * ) for a given far-field loading.

Through-Thickness Distribution of Energy Release Rates (Modes I and II)
The through-thickness distributions of the energy release rates due far-field loadings, σ ∞ y and τ ∞ x y , for a center-crack of length 2a, weakening an infinite plate of finite thickness, 2 h, can be derived by introducing the thickness-wise partial crack closure method as follows: which, on substitution of σ y ( x, 0, z * ) = σ y θ=0 and v(∆a − x, π, z * ) = v| θ= π , obtained from Equations (66) and (64), respectively, for complex roots, and Equations (116) and (114), respectively, for imaginary roots, yields Interestingly, solutions for both complex and imaginary roots reduce to the following: Similarly, which, on substitution of τ x y ( x, 0, z * ) = τ x y θ=0 and u(∆a − x, π, z * ) = u| θ=π , obtained from Equations (79) and (75), respectively, for complex roots, and Equations (125) and (121), respectively, for imaginary roots, yields the following: , f or imaginary roots (148) As before, solutions for both complex and imaginary roots reduce to the following: Equations (146) and (149) reduce to their two-dimensional (plane strain) counterparts [53], by taking D b ( z * ) = 1. It may further be noted that the normalization factor is equal to D b ( z * ) 2 for a given far-field loading. For the special case of a tetragonal single crystal, the above energy release rates reduce to (150)

Crack Deflection Criterion, Based on the Relative Fracture Energy
The important issue of a cleavage plane being deemed easy or difficult can be related to a crack deflection criterion, which is based on the relative fracture energy (or the energy release rate) available for possible "fracture paths" [17]. The deflection or kinking of a crack from the cleavage system 1 to the cleavage system 2 is favored if (however, not iff, i.e., if and only if): in which G i and i , i = 1, 2, are energy release rate and surface energy, respectively, of the ith cleavage system. The above Griffith-Irwin theory-based crack deflection criterion is not accepted as a sufficient condition (albeit being still extremely useful and widely employed, including in the present work) for a cleavage system deemed to be easy or difficult for crack propagation in single crystals, which calls for development of a new conceptual-cum-analytical tool, which would take into account short-range interactions.
Atomistic scale modeling of cracks, however, requires consideration of both the longrange elastic interactions and the short-range chemical reactions. The Griffith theory does not take the latter into account [22]. Secondly and more importantly, fracture criteria derived from equilibrium theories such as the Griffith (thermodynamics-based) energy balance criterion are not equipped to meet the sufficiency condition, because of the prevailing non-equilibrium conditions such as physico-chemical reactions during crack propagation. Hence, such criteria can only be regarded as necessary conditions for fracture [79,83], but not as sufficient [80,84]. The effect of short-range chemical reactions can obviously be encapsulated by atomic scale simulations, such as the investigation of lowspeed propagation instabilities in silicon using quantum-mechanical hybrid, multi-scale modelling due to Kermode et al. [85], which, however, entails extensive computational and other resources. Alternatively, and more importantly, such short-range interactions can also be captured by elastic properties-based parameters (with a few exceptions), such as the anisotropic ratio, A, or, equivalently, the normalized elastic parameter, κ. This is because elastic properties are controlled by various aspects of the underlying structural chemistry of single crystals, such as the Bravais lattice type, bonding (covalent, ionic, and metallic), bonding (including hybridized) orbitals, electro-negativity of constituent atoms in a compound, polarity, etc. [22]. The general theory behind these characteristics pertaining to the structural chemistry of crystals are available in well-known treatises (see, e.g., [86][87][88]). More specifically, the elastic properties of superconducting YBa 2 Cu 3 O 7-δ are strongly influenced by oxygen non-stoichiometry (as well as various structural defects). It is known to crystallize in a defect perovskite structure consisting of layers. When δ = 1, the O(1) sites in the Cu (1) [89]. Furthermore, Ledbetter [90] and Lin et al. [73] measured the elastic constants of polycrystalline YBCO using ultrasonic methods and found that while the "Elastic moduli corresponding mainly to shear modes increase monotonically with oxygen concentration", their counterparts due to "Dilation modes increase up to the values of 6.7 of the oxygen index, after which they begin to decrease"; see also Lubenets et al. [91].
In this connection, it must be noted that the invariant relationship (82), derived earlier in Section 4.1, equating the ratio of mode I to mode II stress intensity factors or its far-field loading counterpart (long range order) to negative times the ratio of the corresponding driving stresses, σ y ( r, θ, z) Mode I / τ x y ( r, θ, z) Mode I I (short range interaction) for the case of an orthorhombic crystal with complex roots (valid for A > 1 or, equivalently, κ > √ c 33 /c 22 ), guarantees the fact of the cleavage system under consideration being difficult, thus further concomitant in ensuring crack deflection or turning from this difficult cleavage system onto a nearby available easy one, violating the self-similarity of crack growth or propagation in the process. This type of behavior of A(or, equivalently, κ), crossing a threshold or critical value that signifies transition from one regime to another, very much establishes its Reynold's number-like character.
It may also be interesting to observe that experimental determination of surface energy, G i , can sometimes be notoriously challenging due to the presence of micro-to-nano scale defects, such as porosity, dislocation, twin boundaries, misalignment of bonds [22,86] with respect to the loading axis, and the like. In contrast, the above-derived invariant relationship (38) demands only measurement of strains or stresses at a point for a given farfield loading, which are, relatively speaking, much easier in comparison to determination of surface energies.

Similarity/Dissimilarity of Present Solutions with Their Isotropic Counterpart
As has been mentioned earlier in Section 4.1, similarity analysis is an effective tool to solve complex problems in the fracture mechanics of single crystals [22,83]. In what follows, such similarity or lack thereof to the present solutions for singular stress fields at crack fronts weakening orthorhombic single crystals with their isotropic counterpart is investigated. Such analyses can lead to a sufficient condition for determination of a cleavage system being easy or difficult for crack propagation.

Isotropic Materials
For an isotropic material, the in-plane displacements (for n = 0) can be given as follows [25,55]: V(x, y, z) = (ik) s b s (z)ρ s e ipψ , in which and giving rise to ψ = π/2 for all positive values of y, for x = 0.

Present Solution Involving Complex Roots
Going back to Equations (20), (21), (30) and (31), the in-plane displacements can be rewritten in the form (for n = 0), by dropping the overhead~over certain relevant variables and constants (in the interest of generality): U(x, y, z) = (ik) s a s (z)(x + py) s = (ik) s a s (z)ρ s e isψ , in which ρ and ψ can be rewritten as follows: Therefore, for an orthorhombic (tetragonal and cubic being special cases) crystal with complex roots when x = 0, for all positive values of y, which completely differs from its isotropic counterpart.

Present Solution Involving Imaginary Roots
Going back to Equations (20), (21), (85) and (86), the in-plane displacements can be rewritten in the form (for n = 0), by dropping the overhead~over certain relevant variables and constants: U(x, y, z) = (ik) s a s (z)(x + py) s = (ik) s a s (z)ρ s e isψ , V(x, y, z) = (ik) s b s (z)(x + py) s = (ik) s b s (z)ρ s e isψ , in which ρ and ψ can be rewritten as follows: again yielding ψ = π/2 for all positive values of y, when x = 0, thus affirming the similarity of crack propagation in such a cleavage system with its isotropic counterpart. Table 1 displays the elastic stiffness constants of orthorhombic (superconductor) and tetragonal (insulator) YBCO single crystals. If otherwise not specified, the elastic stiffness constants are measured at room temperature (Approx. 300 • K). Table 2 shows the elastic stiffness ratio square root, √ c 22 /c 11 , the normalized elastic parameter, κ, nature of the four roots of characteristic equation (complex or imaginary), and the character of the cleavageplane (easy or not) for a (010)[001] through-thickness crack with [100] being the initial propagation direction (see Figure 3), while Table 3 [79]. *** Same as *, except c 12 and c 66 measured by ultrasound by Saint-Paul and Henry [71].     Next, the effect of elastic constants, c ij (especially c 12 and, to a lesser extent, c 66 ), on the nature (i.e., easy or difficult) of a cleavage system in YBCO (YBa 2 Cu 3 O 7−δ ) is discussed. Only three complete sets of elastic constants are available in the Literature accessible to the present author, out of which those due to Ledbetter and Lei [79] are just estimates (marked ** in Table 1), while their experimental counterparts due to Reichard et al. [67] are based on the assumption of tetragonal symmetry; see Table 1. This only leaves the experimental measurements (by resonant ultrasound spectroscopy) due to Lei et al. [1], marked * in Table 1. However, their c 12 value appears to be excessively high. This is because, according to these authors themselves, "no wave speed in the crystal depends only on c 12 , it is no way to estimate it directly." It also is well-known that while c 12 and c 66 can be measured independently by static tests [76], these constants are always coupled in vibrations-based measurements [77,78]. Therefore, in Table 1 of the present investigation, both c 12 and c 66 , measured by ultrasound by Saint-Paul and Henry [71], have been utilized (marked ***) in replacement of their counterparts due to Lei et al. [1] in order to assess the fracture characteristics of YBCO, and to compare them with experiments by Cook et al.       Figure 3.   Table 4, which shows that A (respectively, κ ) value of 1.6266 (resp. 2.2619) have crossed the critical magnitude of 1 (resp. 0.9284). These results predict that (010) , to be easy, which is in agreement with the experimentally observed fracture characteristics of YBCO due to Cook et al. [6], Raynes et al. [9], and Goyal et al. [10], among others; see also Granozio and di Uccio [14] for a summary of the available experimental results. They [14] have also presented approximate theoretical results of fully oxidized YBCOs (δ = 0, 1), and concluded that the three lowest surface energies follow the inequality: γ (001) < γ (100) < γ (010). Furthermore, based on experimental results from transmission electron microscopy [92], X-ray photo-emission microscopy [93], low-energy ion scattering spectroscopy [94], and surface polarity [95] analyses performed on fully oxidized YBCO crystals, these authors [14] have shown that the low energy cut is between the Ba = O and Cu = O planes.

Numerical Results and Discussion
The efficacy of the indentation test has been extensively studied in the brittle fracture Literature [96][97][98]. Lawn [96] and Anstis et al. [97] have presented the following relationship between fracture toughness and the size of a radial crack produced by a Vickers-type sharp indenter: , finally giving rise to the following: in which P, c 0 , E, and H represent the indentation load, equilibrium half-crack length, Young's modulus, and hardness (of an isotropic material), respectively, and § R V denotes a material-independent constant for the Vickers-produced radial crack. Raynes et al. [9], following the lead of Anstis et al. [97], have determined the fracture toughness of monocrystalline YBa 2 Cu 3 O 7−δ , taking into account its anisotropy. Table 8 presents the critical stress intensity factor or fracture toughness (K c = K Ic ) and the critical energy release rate or fracture energy (G c = G Ic ) of the six easy cleavage systems of mono-crystalline superconducting YBCO. It is worthwhile to note here that there is some misconception about the computation of fracture energy, G c , from the corresponding measured value of K c of an anisotropic (e.g., orthorhombic) single crystal in the Literature; see, e.g., Granozio and di Uccio [14]. The factor (c 22 /c 33 ) + χ (see Equation (146) above), and/or its counterparts for other cleavage systems treated in Appendices A-E, are not accounted for in these authors' computations. The energy release rate in an anisotropic (e.g., orthorhombic) single crystal not only varies from one cleavage plane to another, but also varies according to propagation direction. Figure 6a, Figure 9.  Figure 10a,b shows variation of the normalized stress intensity factors, K * (z) = K(z)/K PlaneStrain , through the thickness of an orthorhombic mono-crystalline plate weakened by a through-thickness crack. Variation of the normalized stress intensity factor, K * (z), through the thickness of the same plate weakened by any of the six throughcracks investigated here is identical. Figure 10a shows the through-thickness variation of K * S (z) = K S (z)/K PlaneStrain for a far-field symmetrically distributed hyperbolic cosine load for mode I (stretching) or mode II (in-plane shear), while its skew-symmetric counterpart K * A (z) = K A (z)/K PlaneStrain for mode I (bending) or mode II (twisting) is displayed in Figure 10b. Of special significance is the discontinuity in the stress intensity factor at z* = 0 in the skew-symmetric loading case, shown in Figure 9. Figure 11 shows the corresponding variation of energy release rate, G * , through the top half of the plate thickness. For throughthickness symmetric far-field loading, the crack is expected to grow through thickness in a stable manner until the stress intensity factor or the energy release rate reaches its critical value at mid-thickness. With further increase in the magnitude of far-field loading, unstable crack growth is expected to progressively spread throughout the plate thickness. For skew-symmetric loading, as reported on earlier occasions [55], the bottom half will experience crack closure. Such types of results describing the three-dimensional distribution of stress intensity factors and energy release rates are generally unavailable in the fracture mechanics Literature. release rate reaches its critical value at mid-thickness. With further increase in the magnitude of far-field loading, unstable crack growth is expected to progressively spread throughout the plate thickness. For skew-symmetric loading, as reported on earlier occasions [55], the bottom half will experience crack closure. Such types of results describing the three-dimensional distribution of stress intensity factors and energy release rates are generally unavailable in the fracture mechanics Literature.        ER REVIEW Figure 11. The variation of (mode I, II) energy release rate through thickness due to f hyperbolic load.

Summary and Conclusions
A modified eigenfunction expansion technique, based partly on the sepa z-variable and, in part, on the Eshelby [60]-Stroh [15] type affine transform ployed to derive three-dimensional asymptotic displacement and stress fields ity of the front of a semi-infinite through-thickness crack weakening an in rhombic single crystal plate. Crack-face boundary conditions and those that a on the top and bottom (free) surfaces of the orthorhombic plate are exactly plicit expressions for the singular stresses in the vicinity of the front of the th ness crack, subjected to far-field in-plane mode I and II loadings, are present The present investigation considers six through-crack systems-(010)[0 More importantly, the present approach predicts whe would propagate in its original plane/direction or deflect to a different one. Figure 11. The variation of (mode I, II) energy release rate through thickness due to far-field cosine hyperbolic load.

Summary and Conclusions
A modified eigenfunction expansion technique, based partly on the separation of the zvariable and, in part, on the Eshelby [60]-Stroh [15] type affine transformation, is employed to derive three-dimensional asymptotic displacement and stress fields in the vicinity of the front of a semi-infinite through-thickness crack weakening an infinite orthorhombic single crystal plate. Crack-face boundary conditions and those that are prescribed on the top and bottom (free) surfaces of the orthorhombic plate are exactly satisfied. Explicit expressions for the singular stresses in the vicinity of the front of the through-thickness crack, subjected to far-field in-plane mode I and II loadings, are presented.
The present investigation considers six through-crack systems- (010) (001)[100] with the [010] length direction-weakening orthorhombic YBCO single crystal plates. More importantly, the present approach predicts whether a crack would propagate in its original plane/direction or deflect to a different one. The present study is unique in the sense that such a fracture mechanic criterion is employed for accurate determination of the full set of elastic constants of mono-crystalline YBCO.
The following interesting conclusions can be drawn from the present investigation: (i) Atomistic scale modeling of cracks requires consideration of both the long-range elastic interactions and the short-range chemical reactions. The Griffith thermodynamicbased theory does not take the latter into account, and hence must be regarded as a necessary condition (albeit still extremely useful and widely employed) but not as sufficient.
(ii) The effect of short-range chemical reactions can be adequately captured by the elastic properties-based parameters, such as the anisotropic ratio, A, or, equivalently, the normalized elastic parameter, κ. This is because the elastic properties are controlled by various aspects of the underlying structural chemistry of single crystals, such as the Bravais lattice type, bonding (covalent, ionic, and metallic), bonding (including hybridized) orbitals, electro-negativity of constituent atoms in a compound, polarity, etc.
(iii) More specifically, the elastic properties of superconducting YBa 2 Cu 3 O 7-δ are strongly influenced by oxygen non-stoichiometry (as well as various structural defects).
(iv) It is somewhat paradoxical (or counter-intuitive) that an invariant (more specifically, with respect to θ) relationship, such as Equations (81) and (82) equating the ratio of mode I to mode II stress intensity factors or its far-field loadings counterpart (long range order) to negative times the ratio of the corresponding driving stresses, σ y (r, θ, z) Mode I /τ xy (r, θ, z) Mode I I (short range interaction) for the case of an orthorhombic crystal with complex roots (valid for A > 1 (or, equivalently, κ > √ c 22 /c 11 )), guarantees the fact of the (010)[001] × [100] cleavage system being difficult, thus further concomitant in ensuring crack deflection or turning from this difficult cleavage system onto a nearby available easy one, violating the self-similarity of crack growth or propagation in the process.
(vi) Similarity or dissimilarity of the present asymptotic solutions for a (010)[001] × [100] cleavage system involving complex (A > 1 or, equivalently, κ > √ c 22 /c 11 ) and imaginary roots (A < 1 or, equivalently, κ < √ c 22 /c 11 ) with their isotropic (A = κ = 1) counterparts can lead to a sufficient condition for determination of a cleavage system being easy or difficult for crack propagation.
(vii) The dimensionless parameters, such as the anisotropic ratio, A, or, equivalently, the normalized elastic parameter, κ, can serve as the Holy Grail quantity for an a priori determination of the status of a cleavage system to be easy or difficult, very much akin to Reynold's number for fluid flow problems, crossing a critical value which signifies transition from one regime to another. Here, the anisotropic ratio, A, or, equivalently, normalized elastic parameter, κ, for a (010)[001] × [100] cleavage system, crossing the critical value of 1 or √ c 22 /c 11 , respectively, signifies transitioning from self-similar crack growth or propagation to crack deflection or turning from a difficult cleavage system onto a nearby easy one.
(viii) Just as the introduction of Reynold's number facilitated the design and setting up of experiments in addition to experimental verification of analytical and computational solutions in fluid dynamics, the accuracy and efficacy of the available test results on elastic constants of YBCO single crystals, measured by modern experimental techniques with resolutions at the atomic scale, or nearly so, such as X-Ray diffraction, the ultrasound technique, neutron diffraction/scattering, Brillouin spectroscopy/scattering, resonant ultrasound spectroscopy, and the like, is assessed with a powerful theoretical analysis on crack path stability/instability, in part based on a dimensionless parameter, such as the planar anisotropic ratio, A.
(ix) Experimental determination of surface energy, G i , can sometimes be notoriously challenging, due to the presence of micro-to-nano scale defects, such as porosity, dislocation, twin boundaries, misalignment of bonds with respect to the loading axis, and the like. In contrast, the above-derived invariant relationship, (38), requires only measurement of strains or stresses at a point for a given far-field loading, which are, relatively speaking, much easier in comparison to the determination of surface energies.
(x) The planar anisotropic ratio, A, or, equivalently, the normalized elastic parameter, κ, for YBCO * is larger than 1 or √ c 22 /c 11 , respectively, giving rise to complex roots (of the characteristic equation) for a (010)[001] × [100] through-crack, weakening a YBCO mono-crystalline orthorhombic plate. The same is true for a (100)[001] × [010] crack. These results predict that (010) and (100) are difficult cleavage planes, which are in contradiction with the experimental observations.
(xi) Only for YBCO ***, all the cleavage systems are predicted to be easy, which is in agreement with the experimentally observed fracture characteristics, thus ensuring that a reasonably accurate complete set of nine experimentally determined elastic constants has been arrived at, by employing the present theoretical fracture mechanic approach.
(xii) For tetragonal YBCO T , all six cleavage systems investigated here are found to be difficult, thus completely invalidating the values of the corresponding experimentally determined elastic constants reported by Reichard et al. [67].
(xiii) Finally, generally unavailable results, pertaining to the through-thickness variations of stress intensity factors and energy release rates for a crack corresponding to symmetric and skew-symmetric hyperbolic cosine loads that also satisfy the boundary conditions on the top and bottom surfaces of an orthorhombic mono-crystalline plate under investigation, bridge a longstanding gap in the stress singularity/fracture mechanics Literature.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Singular Stress Fields in the Vicinity of a (010)[001] Through-Crack Front Propagating under Mode I (Extension/Bending) and Mode II (Sliding Shear/Twisting) in [100] Direction
The cleavage plane considered is (010) (Figure 3). Here, the z-axis is placed along the straight crack front, while the coordinates x, y, are used to define the directions along the length of the crack and transverse to it, respectively, in the plane of the plate. u, v, and w represent the components of the displacements in the x, y, and z directions, respectively.
The stress-strain relations for an orthorhombic single crystal are given as follows: where c ij , i, j = 1, ..., 6, denotes the elastic stiffness constants of an orthorhombic monocrystalline plate. σ x , σ y , σ z represent the normal stresses, and τ xy , τ xz , τ yz denote the shear stresses, while ε x , ε y , ε z denote normal strains, and γ xy , γ xz , γ yz represent the shear strains. The three equilibrium equations for a linear elastic solid, made of an orthotropic/ orthorhombic material, can be expressed in terms of the displacement components, u, v, and w, as follows: The characteristic equations for the coupled partial differential Equations (A2)-(A4) can be written as follows: in which the normalized elastic parameter, κ = 1/χ, is given by in which E 1 is Young's modulus in the x direction, G 12 is the shear modulus in the x-y plane, while ν 12 is the major Poisson's ratio in the x-y plane. ν 13 and ν 31 denote the major and minor Poisson's ratios in the x-z plane, while ν 32 represents the minor Poisson's ratio in the y-z plane. κ can also be expressed in terms of the planar anisotropic ratio (in the x-y plane), A, as follows: where A is defined as Equation (A5) has either (a) four complex or (b) four imaginary roots, depending on whether A = 1, κ = 1 represents the degenerate isotropic material case, for which the solution is available in Chaudhuri and Xie [25].
For the extension-bending (mode I) and in-plane shear-twisting (mode II) loadings, it can easily be seen that for orthorhombic single crystals with A < 1 or, equivalently, κ < √ c 22 /c 11 , the (010) plane is the easy cleavage plane, and [100] is the easy propagation direction. Conversely, A > 1 or, equivalently, κ > √ c 22 /c 11 yields complex roots, implying that neither (010) is the easy cleavage plane nor is [100] the easy propagation direction, and the crack will likely deviate from this plane and this direction under mode I/II loadings.

Appendix B. Singular Stress Fields in the Vicinity of a (010)[100] Through-Crack Front Weakening an Orthorhombic Single Crystal under Mode I (Extension/Bending) and Mode II (Sliding Shear/Twisting)
The cleavage plane considered is (010) (Figure 4). Here, the z -axis is placed along the straight crack front, [100], while the coordinates x [001], y [010] are used to define the directions along the length of the crack (propagation direction) and the direction transverse to it, respectively, in the middle plane of the plate. u , v and w represent the components of the displacements in x , y , and z directions, respectively. The stress-strain relations for an orthorhombic single crystal are given by where c ij , i, j = 1, 2, 6, denote the elastic stiffness constants with respect to the rotated coordinate system, x , y (obtained by rotation of 90o about the z-axis): The three equilibrium equations for a linear elastic orthotropic/orthorhombic solid can now be expressed in terms of the displacement functions, u , v , and w, as follows: The characteristic equations for the coupled partial differential Equations (A17)-(A19) can be written as follows: in which the normalized elastic parameter, κ = 1/χ , is given by In which E 2 is y-direction Young's modulus, and G 12 is the shear modulus in the x-y plane, while ν 21 is the minor Poisson's ratio in the x-y plane. ν 31 denotes the minor Poisson's ratio in the x-z plane, while ν 23 and ν 32 represent the major and minor Poisson's ratios, respectively, in the y-z plane. κ . can also be expressed in terms of the planar anisotropic ratio (in the x [001]-y [010] plane), A , as follows: where A is defined as Equation (A20) has either (a) four complex or (b) four imaginary roots, depending on whether A = 1, κ = 1 represents the degenerate isotropic material case, for which the solution is available in Chaudhuri and Xie [25].
For the extension-bending (mode I) and in-plane shear-twisting (mode II) loadings, it can easily be seen that for orthorhombic single crystals with A < 1, or, equivalently, κ < √ c 22 /c 11 = √ c 22 /c 33 , the (010) plane is the easy cleavage plane, and [001] is the easy propagation direction. Conversely, A > 1, or, equivalently, κ > √ c 22 /c 33 yields complex roots, implying that neither (010) is the easy cleavage plane nor is [001] the easy propagation direction, and the crack will likely deviate from this plane and this direction under mode I/II loadings.

Appendix C. Singular Stress Fields in the Vicinity of a (100)[001] Through-Crack Front Propagating under Mode I (Extension/Bending) and Mode II (Sliding Shear/Twisting) in [010] Direction
The cleavage plane considered is (100) ( Figure 5). Here, the z-axis is placed along the straight crack front, [001], while the coordinates x [010], y [100] are used to define the directions along the length of the crack (propagation direction) and the direction transverse to it, respectively, in the middle plane of the plate. u , v and w represent the components of the displacements in x , y , and z directions, respectively. The stress-strain relations for an orthorhombic single crystal are given by where c ij , i, j = 1, 2, 6, denote the elastic stiffness constants with respect to the rotated coordinate system, x , y , z (obtained by rotation of 90 • about the -y-axis): The three equilibrium equations for a linear elastic orthotropic/orthorhombic solid can now be expressed in terms of the displacement functions, u , v , and w , as follows: The characteristic equations for the coupled partial differential Equations (A32)-(A34) can be written as in which the normalized elastic parameter, κ = 1/χ , is given by in which E 3 is z-direction Young's modulus, and G 23 is the shear modulus in the y-z plane, while ν 32 denotes the minor Poisson's ratio in the y-z plane. ν 12 is the major Poisson's ratio in the x-y plane, while ν 13 and ν 31 represent the major and minor Poisson's ratios, respectively, in the x-z plane. κ can also be expressed in terms of the planar anisotropic ratio (in the x [010]-y [100] plane), A as follows: which finally yields in which A is defined as follows: Equation (A35) has either (a) four complex or (b) four imaginary roots, depending on whether (a) A > 1, or, equivalently, κ > c 22 A = 1, κ = 1 represents the degenerate isotropic material case, for which the solution is available in Chaudhuri and Xie [25].
For the extension-bending (mode I) and in-plane shear-twisting (mode II) loadings, it can easily be seen that for orthorhombic single crystals with A < 1, or, equivalently, κ < √ c 22 /c 11 = √ c 11 /c 22 , the (100) plane is the easy cleavage plane (and y [010]-direction is the easy propagation direction). Conversely, A > 1 or, equivalently, κ > √ c 11 /c 22 yields complex roots, implying that neither (100) is the easy cleavage plane nor is [010] the easy propagation direction, and the crack will likely deviate from this plane and this direction under mode I/II loadings.

Appendix D. Singular Stress Fields in the Vicinity of a (100)[010] Through-Crack Front Propagating under Mode I (Extension/Bending) and Mode II (Sliding Shear/Twisting) in [001] Direction
The cleavage plane considered is (100). Here, the z-axis is placed along the straight crack front, [010], while the coordinates x [001], y [100] are used to define the directions along the length of the crack (propagation direction) and the direction transverse to it, respectively, in the middle plane of the plate. u, v and w represent the components of the displacements in x [001], y [100], and z [010] directions, respectively. The stress-strain relations for an orthorhombic single crystal are given by where c ij , i, j = 1, 2, 6, denote the elastic stiffness constants with respect to the transformed coordinate system, x The three equilibrium equations for a linear elastic orthotropic/orthorhombic solid can now be expressed in terms of the displacement functions, u, v and w, as follows: The characteristic equations for the coupled partial differential Equation (121) can be written as follows: in which the normalized elastic parameter, κ = 1/χ, is given by in which E 3 is z-direction Young's modulus, and G 13 is the shear modulus in the x-z plane, while ν 12 and ν 21 denote the major and minor Poisson's ratios, respectively in the x-y plane. ν 13 is the major Poisson's ratio in the x-z plane, while ν 23 represents the major and minor Poisson's ratio in the y-z plane. κ can also be expressed in terms of the planar anisotropic ratio (in the x [001]-y [100] plane), A, as follows: or (b) A < 1, or, equivalently, κ < c 11 A = 1, κ = 1 represents the degenerate isotropic material case, for which the solution is available in Chaudhuri and Xie [25].
For the extension-bending (mode I) and in-plane shear-twisting (mode II) loadings, it can easily be seen that for orthorhombic single crystals with A < 1, or, equivalently, κ < √ c 11 /c 33 , the (010) plane is the easy cleavage plane (and z [001] direction is the easy propagation direction). Conversely, A > 1, or, equivalently, κ > √ c 11 /c 33 yields complex roots, implying that neither (010) is the easy cleavage plane nor is [001] the easy propagation direction, and the crack will likely deviate from this plane and this direction under mode I/II loadings.

Appendix E. Singular Stress Fields in the Vicinity of a (001)[010]. Through-Crack Front Propagating under Mode I (Extension/Bending) and Mode II (Sliding Shear/Twisting) in [100] Direction
The cleavage plane considered is (001). Here, theẑ-axis is placed along the straight crack front, [010], while the coordinatesx [100],ŷ [001] are used to define the directions along the length of the crack (propagation direction) and the direction transverse to it, respectively, in the middle plane of the plate.û,v andŵ represent the components of the displacements inx [100],ŷ [001], andẑ [010] directions, respectively. The stress-strain relations for an orthorhombic single crystal are given by The three equilibrium equations for a linear elastic orthotropic/orthorhombic solid can now be expressed in terms of the displacement functions,û,v andŵ, as follows: The characteristic equations for the coupled partial differential Equations (A63)-(A65) can be written as follows: in which the normalized elastic parameter,κ = 1/χ, is given bŷ in which E 3 is y-direction Young's modulus, and G 13 is the shear modulus in the x-z plane, while ν 21 is the minor Poisson's ratio in the x-y plane. ν 31 denotes the minor Poisson's ratio in the x-z plane, while ν 23 and ν 32 represent the major and minor Poisson's ratios, respectively, in the y-z plane.κ can also be expressed in terms of the planar anisotropic ratio (in thex [100]-ŷ [001] plane),Â, as follows: κ =Â c 33 √ c 11 c 33 + c 13 whereÂ is defined asÂ = 2c 55 √ c 11 c 33 − c 13 .
For the extension-bending (mode I) and in-plane shear-twisting (mode II) loadings, it can easily be seen that for orthorhombic single crystals withÂ< 1 or, equivalently, κ < √ c 33 /c 11 , the (001) plane is the easy cleavage plane (and [100] -direction is the easy propagation direction). Conversely,Â> 1 or, equivalently,κ > √ c 33 /c 11 yields complex roots, implying that neither (001) is the easy cleavage plane nor is [100] the easy propagation direction, and the crack will likely deviate from this plane and this direction under mode I/II loadings.K