PSI - Issue 80

Yohei Sonobe et al. / Procedia Structural Integrity 80 (2026) 368–377 Y. Sonobe et al. / Structural Integrity Procedia 00 (2023) 000–000

371

4

2.2. Basic density function and calculation of mode I SIF K I The magnitude of the body force doublet at an arbitrary source point Q ( η,ζ ) on the crack face, γ ( η,ζ ) , is directly proportional to the three-dimensional crack opening displacement (COD) at the same point, δ u ( η,ζ ) (Noda and Oda, 1992):

(1 − 2 ν )(1 + ν ) E (1 − ν )

(6)

γ ( η,ζ )

δ u ( η,ζ ) =

where E is the Young’s modulus. In the BFM, to achieve a highly accurate analysis, γ ( η,ζ ) is replaced by the product of a fundamental density function, B ( η,ζ ) , based on a reference solution for the COD and an unknown weight function, W ( η,ζ ) , which accounts for three-dimensional effects: γ ( η,ζ ) = B ( η,ζ ) W ( η,ζ ) (7) B ( η,ζ ) is derived from a 2D plane strain analog of a deep surface crack. This reference solution considers a crack subjected to the same n -th order pressure distribution as given by Eq. (1). The COD for this 2D plane strain case is denoted δ u ∞ ( η, n ) , which in turn defines the fundamental density function.

E (1 − ν ) (1 − 2 ν )(1 + ν )

(8)

δ u ∞ ( η, n )

B ( η, n ) =

The reference COD, δ u ∞ ( η, n ) , is calculated using the complex variable method (England, 2012) as

tanh − 1   

a 2 − t η    dt

( a 2 − t 2 )( a 2 − η 2 )

π E

t a

a − a

4(1 − ν 2 ) p 0

n

(9)

δ u ∞ ( η, n ) =

Once the unknown weight function W ( η,ζ ) in Eq. (7) is determined, the mode I SIF, K I ( ζ ) , along the crack front is calculated from the following equation, using W ( ± a ,ζ ) , and the stress intensity factor of the reference solution, K ∞ I ( n ) : K I ( ζ ) = W ( ± a ,ζ ) K ∞ I ( n ) (10) where K ∞ I ( n ) can be defined as follows:

E √ π a 4(1 − ν 2 )

δ u ∞ ( η, n ) a 2 − η 2

(11)

K ∞ I ( n ) =

lim η →± a

2.3. Discretization of the governing integral equation By substituting Eq. (7) into the governing integral equation Eq. (5), the following equation is obtained. a − a B ( η, n ) ∞ 0 W ( η,ζ ) σ γ xx ( y , z ,η,ζ ) d ζ d η = − p 0 y a n (12) Here, when B ( η, n ) is defined from the 2D reference solution as described in the previous section, it is known that W ( η,ζ ) rapidly approaches p 0 in the deep region of the crack ( ζ a ). Utilizing this property, in this analysis, the domain for discretization is limited to a cutoff depth ζ = ζ thresh , and it is assumed that W ( η,ζ ) = p 0 holds for the deeper region ( ζ > ζ thresh ). For the specific discretization, only the region where y ≥ 0 is considered, taking advantage of the symmetry of the crack face. The analysis domain, which is a half-plane of the crack face ( 0 ≤ y ≤ a , 0 ≤ ζ ≤ ζ thresh ), is divided into N e rectangular elements, Ω j ( j = 1 , . . . , N e ). The weight function W ( η,ζ ) is assumed to be a constant value W j within each element, and the center of each element, ( y i , z i ) , is used as a collocation point where the boundary condition is satisfied.