Effect of random variation in input parameter on cracked orthotropic plate using extended isogeometric analysis (XIGA) under thermomechanical loading

2.1 Fracture analysis

The orthotropic structure with a crack is defined within global cartesian coordinates (X, Y), aligning its elastic symmetry axes with local Cartesian axes (x1, y1). The crack tip is represented using local polar coordinates (r, θ). as illustrated in Figure 1.

Figure 1

Cracked orthotropic body experiencing traction forces. This area is enclosed by the contour r and contains arbitrary boundary conditions within the designated area A.

The body encounters various forces alongside different displacement and traction boundary conditions. Expressing orthotropic stress–strain relationship by:

The relevant compliance coefficients of body are C _ij (i, j = 1, 2, 6) in local Cartesian coordinates. The stress–strain equation is given as [49]

Further, in-plane elastostatic set of equations is given as follows:

The partial differential equation governing the fourth order for an orthotropic body, derived from compatibility and equilibrium considerations, can be stated as follows: [37,38]

Complex or purely imaginary roots of Eq. (7) and are derived as ( μ k = μ k x 1 + i μ k y 1 ) which always occur in conjugate pairs as μ 1 , μ 1 ¯ and μ 2 , μ 2 ¯ .

Complex variables and analytical functions ( z k = x 1 + μ k y 1 k ( k = 1 , 2 ) ) are employed to delineate the two-dimensional displacement and stress fields surrounding the crack tip [50] The pure mode I stress components given by,

and the corresponding displacements are given as

Pure mode II stress and displacement are:

and the displacement is given as

In this scenario, K _I denotes the mode I stress intensity factors (SIF), and K _II denotes the mode II stress intensity factors (SIF). “Re” indicates the real part of the statement, while p _k and q _k are further explained below.

2.2 XIGA formulation

To effectively address problems involving discontinuities such as cracks, integrating concepts from XFEM with IGA significantly enhances their capabilities. XFEM is adept at handling static and dynamic discontinuities within structures, while IGA demonstrates precision and efficiency in analyzing complex geometries. Merging these methodologies, known as XIGA, allows for defining the entire crack independently of the mesh. XIGA essentially represents IGA approximation enriched by the application of partition of unity (PU), which is derived from the XFEM formulation.

To model singular fields and discontinuities, XIGA improves the local isogeometric approximation by introducing extra degrees of freedom at certain points near crack locations within the basic IGA model. Enrichment functions are used to contribute to overall approximation. In the XIGA method, crack-tip enrichment functions are utilized for capturing singularity within the stress fields at the crack tip, whereas Heaviside functions are employed for the representation of the crack face (Figure 2).

Figure 2

Sketch depicting the arrangement of various methods.

The generalized expression for the displacement approximation at the designated control point corresponding to ξ i is:

where u _i signifies parameters related to the particular control point; n _en refers to the number of basis functions for each element; R i p denotes the rational B-spline basis functions; n _cf denotes the number of basis functions entirely intersected by the crack face; n _ct represents the number of basis functions partially intersected by the crack tip; n _h is the count of basis functions partially intersected by the hole; and n _in indicates the number of basis functions partially intersected by the inclusion.

Additional control points influence enriched degrees of freedom as follows:

1. a _j accounts for the Heaviside function H ( ξ ) ;

2. b k a influences degrees of freedom associated with asymptotic crack-tip enrichment functions β α ;

3. c j and d j contribute to control points associated with specific functions, namely χ(ξ) and ψ(ξ) respectively, characterizing material interfaces.

Blending functions like B h and B t are employed for the Heaviside function and crack-tip function, respectively. The Heaviside function H ( ξ ) equals +1 when the parametric coordinate ξ corresponds to a particular Gauss point situated above the crack face. Conversely, when ξ lies below the crack face, H ( ξ ) equals −1.

Strain is comprised of two components in a thermomechanical problem.

Here, ε t is the total strain, ε m is the mechanical strain, ε th is the thermal strain, where ε th can be written as for plane stress state,

and for plane strain state,

Thermal expansion coefficient is denoted as α i j , and Poisson’s ratio is denoted by μ i j . According to the generalized Hooke’s law, the mechanical part of strain, ε m is shown below.

The discretized form can be written as

The displacement vector u h encompasses standard as well as additional degrees of freedom.

The displacement u h represents the standard degrees of freedom, while ‘a’ and ‘bα’ (where α ranges from 1 to 4) symbolize additional degrees of freedom utilized for the crack face and fields surrounding the crack tip. The mechanical force vector f and global stiffness matrix K are both assembled from their respective element’s contributions.

Similarly, the discretize form of thermal is

where T h is the temperature vector

Thermal force vector f th and thermal stiffness matrix Q are formed using a different matrix of basis function derivatives.

2.3 Stochastic response of MMSIF for different input random variables

This article explores the stochastic behavior of MMSIFs using the XIGA approach. It employs two probabilistic methods, namely SOPT, which relies on perturbation techniques, and direct MCS, a sampling-based method. SOPT, employing Taylor series expansion, establishes connections between certain features of random parameters and their corresponding responses. But, this method’s applicability is confined because of its dependency on lower-order polynomials. The subsequent section provides a numerical formulation outlining these methods.

3.1 Edge crack with various orientations of orthotropic axis subjected to tensile loading

We apply the proposed method to analyze an edge crack in a rectangular plate. The plate has a height-to-width ratio of 2 and a crack length-to-width ratio of 0.45, as shown in Figure 3(a). The material used is graphic epoxy with properties that are described below:

Figure 3

(a) The loading and geometry for a single edge crack with various orientations of orthotropic axis in rectangular plate, (b) NURBS mesh with enriched nodes, and (c) integral domain to calculate SIF.

E ₁ = 114.8 GPa, E ₂ = 11.7 GPa, G ₁₂ = 9.66 GPa and µ ₁₂ = 0.21.

The mesh features a uniform distribution of nodes in extended isogeometric analysis. NURBS basis functions employed are (p = q = 3) cubic, as depicted in Figure 3(b), with enriched nodes shown in Figure 3(c). The integral domain for calculating SIF is illustrated. Our results are compared with Aliabadi and Mohammadi. Figure 4 presents the SIF values for various α angles, obtained using the proposed XIGA method, Aliabadi’s FEM, and Mohammadi’s XFEM analyses.

Figure 4

Results of normalized SIFs for K _I and K _II of a single edge crack in a rectangular plate with various orientations of the orthotropic axes, at angles of α = 0° 30° 45° 60° 90° and ( K I ¯ = K I / σ πa and K II ¯ = K II / σ π a ) .

The maximum difference in mode I SIF is observed at α = 0° where a 4.3% error is observed, and for mode II, this occurs when α = 20° where a 3.6% error is observed. The results obtained by this method demonstrate that for both modes, there is an increase in normalized SIF value from α = 0° to α = 45° and then are decreased from α = 45° to α = 90°.

The effect of crack length on K _I and K _II was investigated for α = 0° orientation which is shown in Table 1 below which shows an exponential increment in K _I with an increase in crack length, and minor variation was found in K _II values.

Further study was carried out for different loading on the top edge of the plate (Tensile, Shear, Combined, and Triangular), and values of Normalized SIF, K _I, and K _II were recorded for different orientation angles (α = 0° to α = 90°), which are presented in Table 2.

The crack length taken into consideration was a/w = 0.45 and h/w = 2 as depicted in Figure 3(a), and values of K _I opening mode are more in combined loading as compared with only shear loading which is more than only tensile loading for all orientation angles. K _II shearing mode was very low in tensile loading but in shear and combined loading it was comparatively higher for all the orientation angles. For a triangular tensile load, it was obvious that the values of K _I and K _II would give lower values as the total load was half that of the uniformly distributed load.

Results of normalized K _I and K _II for different orientation angles under tensile, shear, and combined loading are depicted in Figures 5–7.

Figure 5

Normalized K _I & K _II for different orientation angles under tensile loading.

Figure 6

Normalized K _I & K _II for different orientation angles under shear loading.

Figure 7

Normalized K _I & K _II for different orientation angles for combined tensile and shear loading.

3.2 Isotropic plate with center crack under thermomechanical loading

Here, we investigate a square plate with a center crack through numerical simulations under thermomechanical loading, employing the XIGA approach. The dimensions considered are 2,000 mm × 2,000 mm, and the domain is discretized using a 50 × 50 control net. We use a NURBS basis function order of three for the simulations.

The plate contains a preexisting center crack with a length of 600 mm, as shown in Figure 8. The material properties considered are as follows:

• elastic modulus (E) = 200 GPa;

• Poisson ratio ( ν ) = 0.3;

• thermal expansion coefficient ( α ) = 11.7 × 10 − 6 C − 1 .

Figure 8

Center cracked square plate under thermomechanical loading.

The scenario involves a square plate under tensile load of 30 MPa applied at the top edge which is perpendicular to the crack. Additionally, temperature boundary conditions are applied: 0°C at the crack face and 10°C at the top edge of the plate. The bottom edge of the plate is constrained, as depicted in Figure 8.

For each step of crack extension, a crack increment of 100 mm is applied, extending proportionally from both the left and right tips. The stress intensity factors are then calculated using the domain-based interaction integral approach.

Figure 9 illustrates the stresses in x-direction (u _x ) for a/w = 0.4 subjected to thermomechanical loads, while Figure 10 shows stresses in the y-direction. Additionally, Figures 11 and 12 display the displacements along the x direction and y direction. The variation of SIFs (K _I and K _II) with the extension of the right crack tip is shown in Figure 13. The graph indicates that K _I is higher than K _II. Consequently, under combined loading conditions, the center crack in the plate exhibits opening mode behavior.

Figure 9

Stresses along x-axis.

Figure 10

Stresses along y-axis.

Figure 11

Displacement along x-axis.

Figure 12

Displacement along y-axis.

Figure 13

SIFs (K _I and K _II) variation of the plate subjected to thermomechanical loads is presented (validation study).

3.3 Orthotropic plate with center crack under mechanical, thermal, and thermomechanical uniaxial loading

We further analyze a center crack in a square plate under thermomechanical loads, using the XIGA approach for an orthotropic material. The square plate has initial dimensions of 2,000 mm × 2,000 mm, and the domain is discretized using a control net of size 50 × 50.

The square plate contains a preexisting center crack with a length of 600 mm, as depicted in Figure 14. This crack is analyzed when subjected to thermomechanical loading conditions using the XIGA method. The material properties considered for the analysis are as follows:

• elastic modulus (E11) = 114.8 GPa;

• elastic modulus (E22) = 11.7 GPa;

• shear modulus (G12) = 9.66 GPa;

• Poisson ratio ( ν 12 ) = 0.21;

• thermal expansion coefficient ( α 11 ) = 11.7 × 10 − 6 C − 1 ;

• thermal expansion coefficient ( α 22 ) = 2 × 10 − 5 C − 1 .

Figure 14

Center-cracked orthotropic square plate under thermomechanical loading.

Boundary conditions are as follows:

• temperature difference = 100°C;

• tensile stress on top edge = 30 MPa.

Table 3 presents the normalized SIF, K _I and K _II, for different crack lengths. The results show that K _I is consistently higher than K _II. This indicates that the center crack in the square plate primarily experiences opening mode under combined loading conditions, even for orthotropic materials. Furthermore, the analysis suggests that thermomechanical loading is more critical than pure mechanical or pure thermal loading. This is supported by the observation that as the crack length increases, the SIF values increase for all loading types.