PSI - Issue 80

Luke Wyatt et al. / Procedia Structural Integrity 80 (2026) 31–42

35

L. Wyatt et. al. / Structural Integrity Procedia 00 (2023) 000–000

5

where ˆ σ 0 = 0 yy is the initial in plane load on the plate [11]. The global buckling problem can be solved as the eigenproblem ( K − λ i K G ) X i = 0 , σ 0 xx τ yx σ 0 xy τ 0

(19)

th buckling multiplier, and X

th buckling mode, and K

where λ i is the i

i is the i

G is the sti ff ness matrix of U G .

The element geometric sti ff ness matrix, K G | E can be found as

T T

T

t 3 12

T dE ˜ Π

˜ Π

T dE ˜

Π + ˜ Π

t ( ∇· N P

T ˆ σ 0 ( ∇· N P w )

T ˆ σ 0 ( ∇· N P

( ∇· N P

K G | E =

w )

θ x )

θ x )

(20)

E

E

t 3 12

T dE ˜ Π

+ ˜ Π

( ∇· N P

T ˆ σ 0 ( ∇· N P θ y )

θ y )

(21)

E

where N P and N P θ y are the polynomial trial functions of each displacement. The projection matrix ˜ Π is known from the construction of the sti ff ness and mass matrices. K G | E is not stabilised, following the formulation of the VEM for buckling for Kirchho ff plates [12]. The global geometric sti ff ness matrix can be assembled, and Equation 19 can be solved the find the buckling loads and mode shapes. w , N P θ x

4.3. Bending Stress Intensity Factor

For a Mindlin plate containing a crack, the bending stress intensity factor in mode 1, K 1 B can be found by use of the J-integral technique [13]. The mode 1 J-integral, J 1 B can be written as an oriented path integral around a crack tip, where the boundary is denoted Γ and the enclosed area Ω C , (Figure 1a). For a purely mode 1 crack, the local crack tip directions coincide with the global coordinate system. This means that J 1 B can be written J 1 B = Γ ( W − qw ) n x d Γ − Γ M R n · ∂ θ ∂ x + Q R n · ∂ w ∂ x d Γ+ Ω C ∂ q ∂ x wd Ω C , (22) where W is the strain energy, n x is the x component of the outward normal vector on Γ , n [13]. This is valid for cracks with traction free crack faces. As the displacements vary across the element linearly, the strains are constant in each element and can be found as usual [2]. The derivatives of θ with x can be approximated in each element as ∂ x . The implementation of the J-integral was chosen to keep the procedure a simple as possible, a circular path of radius r , was chosen as Γ and the number of integration points were defined. The value of the functions in Equation 22 were calculated for each point, using the interpolation provided by the element which the point was in. The path integral was then found using the trapezium rule (Figure 1b). The area integral can be found by approximating the load as standard in VEM [2], which leaves only the integral of polynomial functions over each element to be calculated. The integral over Ω C can be found by summing the contributions for each E ∈ Ω C . Once J 1 B is found, K 1 B can be found as K 1 B = Et 3 12 J 1 B . (24) ∂ θ ∂ x ≈ ∂ N θ P ∂ x ˜ Π ˜ V E , (23) where ˜ V E are the results of the generalised displacements for the element. The same procedure can be used to find ∂ w