PSI - Issue 81

Oleh Yasniy et al. / Procedia Structural Integrity 81 (2026) 244–250

247

The sign in Eq. (5) was determined based on the sign of / II I K K . A boundary element method was used to perform calculations for two-dimensional cracks, without accounting for stress variation through the thickness. To include the effect of through-thickness stresses, a parameter max C was introduced to limit the surface length of the crack. It was assumed that cracks do not change their shape (aspect ratio) during growth. Therefore, a small value of max C results in short, shallow cracks, while a large max C produces long and deep cracks. Additionally, the distance between the tips of two neighboring cracks could not exceed max C , regardless of the magnitude of the stress intensity factor. This approach allows the model to approximate three-dimensional effects in a two-dimensional framework by controlling crack aspect ratios and limiting interactions between closely spaced cracks. By adjusting max C , the model can represent realistic variations in crack depth and length, providing more accurate predictions of crack network evolution under thermo-fatigue loading. 2.3. Boundary element modeling of stress state To compute the stress fields at crack initiation nuclei and to determine the stress intensity factors, the dual boundary element method (DBEM) (Pasternak et al., 2013) was applied. In this method, based on Somigliana’s identity for plane elasticity problem s, a system of 2 n equations is constructed: n equations, as in the classical BEM, for displacements, and an additional n equations for stresses, obtained by differentiating Somiglian a’s identity. Here, n is the problem’s dimensionality. Thus, for fracture mechanics problems, the integral equations of the dual BEM take the following form (Pasternak et al., 2013): • when the collocation point y lies on a smooth boundary  of the body:

1 2

  y

      , d   x y x x

      ,   x y x x   , x y x x   ( ) u d   j T





u

U t

T u d



i

ij

j

ij





(6)

.

  , U t

     d





( ) x y x x

;



ij

j

ij

j

C 

C 





• when the collocation point y lies on a smooth surface

C  of a crack:

1 2

    y u

      , d   x y x x

      , u d   x y x x   , x y x x   ( ) u d   j





U t

T

i

ij

j

ij





  , U t

     d





( ) x y x x

,

T



ij

j

ij

k

C 

C 





(7)

.

1 2

    y t

        , d    x y y x x     ,  x y y x x      ( ) j k t d

        , n u d    x y y x x     , n  x y y x x .   ( ) u d   j k





D n t

S

i

ijk

ijk









D n

S



ijk

j

k

ijk

j

k

C 

C 





Here ij U , ij T , ijk D , ijk S are the kernels of the boundary integral equations (Pasternak et al., 2013), i u and i t are displacement and traction components, i i i u u u      , i i i t t t      , i i i u u u      , i i i t t t      ; i n are the components of the normal vector to the surface. The symbols “+” and “−” indicate quantities associated with the surfaces C   and C   created by the crack faces. The indices correspond to the projections of vectors onto the axes of the global coordinate system 1 2 Ox x . In the formulas, Einstein summation convention is applied over repeated indices. For efficient evaluation of the singular and hypersingular integrals arising in the problem, as well as for accounting for the square-root singularity of the solution at crack tips, numerical quadrature schemes, nonlinear polynomial mappings, and special basis functions proposed in (Pasternak et al., 2013) were employed. This approach also made it possible to determine the stress intensity factors (SIFs) of the cracks with high accuracy. 2.4. Numerical modeling The procedure for modelling the initiation and subsequent growth of thermo-fatigue cracks was organized as follows. At the preparatory stage, each potential crack initiation nucleus was randomly assigned a value of the maximum admissible damage level ( ) i R according to a uniform distribution. The node with the smallest value of ( ) i R was then identified, and a crack with an initial length of 0.2 mm was introduced at this location. In parallel, the number of loading cycles required for the initiation of this crack was determined based on the assigned damage parameter and the local stress – strain state. The main cycle of the algorithm consisted of the following operations. The input parameters were the number of loading cycles performed in a single iteration and the maximum allowable distance max C between crack tips. For the system of cracks within the body, the range of the stress intensity factor (SIF), the parameter max C , and the crack increment length were evaluated, with the latter being determined on the basis of the Paris law (4). If the crack length exceeded the value of the parameter max C , further execution of the main cycle no longer led to its growth, and the crack was treated as having reached its limiting size under the given modelling assumptions.