Determination of Stresses in Composite Plates with Holes and Cracks Based on Singular Integral Equations

Abstract

The problems of determination of stresses at crack in bounded plates with holes of different shapes under the action of concentrated forces or distributed forces at its boundary are considered. The study is performed by the singular integral equation method. They were determined based on the established interdependences between the Lekhnitskii potentials and the stress and strain. The numerical method for solving integral equations is developed based on the quadrature method for the systems of holes and cracks. The eigen solutions of the problem were taken into account in this method. The research of stresses at cracks in samples which are used in experimental studies of crack fracture resistance was performed.

1. Introduction

The boundary integral equation method (BIEM) is widely used to study the stress-strain state (SSS) of anisotropic plates with holes [1, 2, 3] and cracks [4, 5, 6, 7, 8, 9]. The integral equations for anisotropic plates are usually determined based on the Somigliana identity. Such equations for plates with given stresses at the boundaries of the plate are hypersingular. At the same time, the same problem for isotropic plates is reduced to singular integral equations [10, 11], for which simple numerical algorithms for solving with given precision are obtained.

In [12, 13], the simple dependencies between the Lekhnitskii complex potentials and stress and strain are obtained. In a simple form based on them and the Cauchy theorem, the integral equations are written for anisotropic plates with holes [12, 13] and cracks [14, 15, 16]. We will use the established dependencies for the construction and regularization of integral equations for anisotropic plates with holes and cracks.

For conducting experimental studies of crack fracture resistance on experimental samples in relation to isotropic materials, theoretical estimates for stresses at cracks are performed.

For such materials, the stresses in samples of different shapes with cracks under the action of stretching or compressing concentrated forces are studied in detail [10]. The experimental samples for the experimental determination of the characteristics of crack fracture resistance of various types of materials are made based on performed studies. We perform similar studies for composite samples.

2. The integral representations for anisotropic plates with holes and cracks

We consider a plate which is weakened with a system hole with boundaries $L_1,\dots ,L_J$ (j = 1, …, J), and cracks are placed along curves ${\mathrm{\Gamma }}_k\left(k=1,\dots ,K\right)$ . The $L_0$ is the outer boundary of plates. Assume ( Figure 1 ) that a plate is loaded with concentrated forces (X_j , Y_j ), j = 1, ..., M acting at the points (a_j , b_j ); tractions $\left(X_T,Y_T\right)$ are applied to the crack edges, which are accepted the same on its opposite edges; and tractions $\left(X_L,Y_L\right)$ are applied to the boundaries of the holes and plate.

Figure 1.

Scheme of the problem.

2.1 Governing equations

Let us start from the Lekhnitskii complex potentials $\mathrm{\Phi }\left(z_1\right),\mathrm{\Psi }\left(z_2\right)$ , where $z_j=x+s_{j}y$ and $s_j,j=1,2$ are roots with positive imaginary part of the characteristic equation $\mathrm{\Delta }\left(s\right)=0$ [10]:

where

$$\mathrm{\Delta }\left(s\right)={\alpha }_{11}{s}^4-2{\alpha }_{16}{s}^3+\left(2{\alpha }_{12}+{\alpha }_{66}\right)s^2-2{\alpha }_{26}s+{\alpha }_{22}$$ andE1

${\alpha }_{\mathit{ij}}$ are elastic compliances which are included in the Hooke’s law [10]:

$${\epsilon }_x=a_{11}{\sigma }_x+a_{12}{\sigma }_y+a_{16}{\tau }_{\mathit{xy}},{\epsilon }_y=a_{12}{\sigma }_x+a_{22}{\sigma }_y+a_{26}{\tau }_{\mathit{xy}},{\gamma }_{\mathit{xy}}=a_{16}{\sigma }_x+a_{26}{\sigma }_y+a_{66}{\tau }_{\mathit{xy}},$$

where ${\epsilon }_x,{\epsilon }_y,{\gamma }_{\mathit{xy}}$ are strains and ${\sigma }_x,{\sigma }_y,{\tau }_{\mathit{xy}}$ are stresses.

Consider an arbitrary path $\mathrm{\Gamma }$ , which belongs to the domain D occupied by the plate, and select a positive direction of traversal ( Figure 2 ).

Figure 2.

$\overrightarrow{q_{\mathrm{\Gamma }}}$ is the stress vector at plane $\mathit{AB}$ .

Then introduce in consideration the stress vectors $\overrightarrow{q_{\mathrm{\Gamma }}}$ at the plane tangent to the curve. The normal to it is located right relative to the selected direction of traversal. The projections $\left(X_{\mathrm{\Gamma }},Y_{\mathrm{\Gamma }}\right)$ of stress vectors $\overrightarrow{q_{\mathrm{\Gamma }}}$ and derivatives of displacements $\left(u,v\right)$ with respect to an arc coordinate at the curve through Lekhnitskii complex potentials are determined by the formula [17]:

$$Y_{\mathrm{\Gamma }}=-2Re\left[\mathrm{\Phi }\left(z_1\right)z_1^{\prime }+\mathrm{\Psi }\left(z_2\right)z_2^{\prime }\right],X_{\mathrm{\Gamma }}=2Re\left[s_1\mathrm{\Phi }\left(z_1\right)z_1^{\prime }+s_2\mathrm{\Psi }\left(z_2\right)z_2^{\prime }\right],$$ E2

$$u\prime =2Re\left[p_1\mathrm{\Phi }\left(z_1\right)z_1^{\prime }+p_2\mathrm{\Psi }\left(z_2\right)z_2^{\prime }\right],v\prime =2Re\left[q_1\mathrm{\Phi }\left(z_1\right)z_1^{\prime }+q_2\mathrm{\Psi }\left(z_2\right)z_2^{\prime }\right],$$ E3

where $u\prime =\mathit{du}/\mathit{ds},v\prime =\mathit{dv}/\mathit{ds}$ and $z_j^{\prime }=\mathit{dx}/\mathit{ds}+s_j\mathit{dy}/\mathit{ds}$ , where ds is a differential of arc at $\mathrm{\Gamma }$ .

The stress vectors $q_{\mathrm{\Gamma }}\left(z\right)=X_{\mathrm{\Gamma }}+{\mathit{iY}}_{\mathrm{\Gamma }}$ at path $\mathrm{\Gamma }$ are determined using the formulas (2) by the formula:

$$q_{\mathrm{\Gamma }}=\left(s_1-i\right)z_1^{\prime }\mathrm{\Phi }\left(z_1\right)+\left(\overline{s_1}-i\right)\overline{z_1^{\prime }}\overline{\mathrm{\Phi }\left(z_1\right)}+\left(s_2-i\right)z_2^{\prime }\mathrm{\Psi }\left(z_2\right)+\left(\overline{s_2}-i\right)\overline{z_2^{\prime }\mathrm{\Psi }\left(z_2\right)}.$$ E4

Assume that the vectors $\left(X,Y\right)$ and $\left(u,v\right)$ are known at path $\mathrm{\Gamma }$ . Then based on Eqs. (2) and (3) at $\mathrm{\Gamma }$ one has [12, 15]

$$\mathrm{\Phi }\left(z_1\right)=\frac{-v\prime +s_{1}u\prime +p_{1}X+q_{1}Y}{{\mathrm{\Delta }}_1{z}_1^{\prime }},\mathrm{\Psi }\left(z_2\right)=\frac{-v\prime +s_{2}u\prime +p_{2}X+q_{2}Y}{{\mathrm{\Delta }}_2{z}_2^{\prime }}$$ E5

where ${\mathrm{\Delta }}_j=\mathrm{\Delta }\prime \left(s_j\right),j=1,2$ .

2.2 Integral equations for anisotropic bounded plate with holes and cracks

Let us write a general solution of the problem based on [12, 15] through the Lekhnitskii potentials in the form

$$\mathrm{\Phi }\left(z_1\right)=\underset{L}{\int }\left[u\prime {\mathrm{\Phi }}_1\left(z_1,t_1\right)+v\prime {\mathrm{\Phi }}_2(z_1,t_1)\right]\mathit{ds}$$

$$+\underset{\mathrm{\Gamma }}{\int }\left[g_1^{\prime }{\mathrm{\Phi }}_1\left(z_1,t_1\right)+g_2^{\prime }{\mathrm{\Phi }}_2(z_1,t_1)\right]\mathit{ds}+{\mathrm{\Phi }}_S\left(z_1\right)+{\mathrm{\Phi }}_{\mathrm{\Delta }}\left(z_1\right),$$ E6

$$\mathrm{\Psi }\left(z_2\right)=\underset{L}{\int }\left[u\prime {\mathrm{\Psi }}_1\left(z_2,t_2\right)+v\prime {\mathrm{\Psi }}_2(z_2,t_2)\right]\mathit{ds}$$

$$+\underset{\mathrm{\Gamma }}{\int }\left[g_1^{\prime }{\mathrm{\Psi }}_1\left(z_2,t_2\right)+g_2^{\prime }{\mathrm{\Psi }}_2(z_2,t_2)\right]\mathit{ds}+{\mathrm{\Psi }}_S\left(z_2\right)+{\mathrm{\Psi }}_{\mathrm{\Delta }}\left(z_2\right),$$

where $L=L_0+L_1+\dots +L_J$ , $\mathrm{\Gamma }={\mathrm{\Gamma }}_1+{\mathrm{\Gamma }}_2+\dots +{\mathrm{\Gamma }}_K$ , s is an arc coordinate, and ${\mathrm{\Phi }}_{\mathrm{\Delta }}\left(z_1\right)$ and ${\mathrm{\Psi }}_{\mathrm{\Delta }}\left(z_2\right)$ are the known functions, which are determined by the following formulas:

$${\mathrm{\Phi }}_{\mathrm{\Delta }}\left(z_1\right)=\underset{L}{\int }\left[X_L{\mathrm{\Phi }}_3\left(z_1,t_1\right)+Y_L{\mathrm{\Phi }}_4(z_1,t_1)\right]\mathit{ds},{\mathrm{\Psi }}_{\mathrm{\Delta }}\left(z_1\right)=\underset{L}{\int }\left[X_L{\mathrm{\Psi }}_3\left(z_1,t_1\right)+Y_L{\mathrm{\Psi }}_4(z_1,t_1)\right]\mathit{ds},$$ E7

$${\mathrm{\Phi }}_j=\frac{A_j}{t_1-z_1},{\mathrm{\Psi }}_j=\frac{B_j}{t_2-z_2},$$ E8

$$A_1=-\frac{{\mathit{\text{is}}}_1}{2\pi {\mathrm{\Delta }}_1},A_2=\frac{i}{2\pi {\mathrm{\Delta }}_1},A_3=-\frac{{\mathit{ip}}_1}{2\pi {\mathrm{\Delta }}_1},A_4=-\frac{{\mathit{iq}}_1}{2\pi {\mathrm{\Delta }}_1},$$

$$B_1=-\frac{{\mathit{\text{is}}}_2}{2\pi {\mathrm{\Delta }}_2},B_2=\frac{i}{2\pi {\mathrm{\Delta }}_2},B_3=-\frac{{\mathit{ip}}_2}{2\pi {\mathrm{\Delta }}_2},B_4=-\frac{{\mathit{iq}}_2}{2\pi {\mathrm{\Delta }}_2}.$$

Here, $u\prime ,v\prime$ are the values of the derivatives of the displacements with respect to the arc coordinate at the boundary of the plate and holes, $g_1=u^{+}-u^{-},g_2=v^{+}-v^{-}$ are the displacements discontinuity at the cracks, $u^{\pm },v^{\pm }$ are limit values of displacements in the approach to the section at the left and the right relative to the selected direction, and the potentials ${\mathrm{\Phi }}_S,{\mathrm{\Psi }}_S$ correspond to the concentrated forces and have the form [12]:

$${\mathrm{\Phi }}_S\left(z_1\right)=\frac{i}{2\pi {\mathrm{\Delta }}_1}\sum _{j=1}^M\left(p_1{X}_j+q_1{Y}_j\right)\frac{1}{z_1-z_{1j}},{\mathrm{\Psi }}_S\left(z_2\right)=\frac{i}{2\pi {\mathrm{\Delta }}_2}\sum _{j=1}^M\left(p_2{X}_j+q_2{Y}_j\right)\frac{1}{z_2-z_{2j}},$$ E9

where in $z_{\mathit{kj}}=a_j+s_k{b}_j$ , j = 1, 2, and k = 1, 2.

Note that when the boundary is traction free, then ${\mathrm{\Phi }}_{\mathrm{\Delta }}={\mathrm{\Psi }}_{\mathrm{\Delta }}=0$ .

Let us substitute the potentials (6) into the formulas (4) for projections of stress vectors determined at the boundaries path $L$ and $\mathrm{\Gamma }$ . Using Plemelj-Sokhotski formula, we obtain a system of integral equations [12, 15]:

$$\begin{array}{l}\underset{L}{\int }\left[u\prime \left(s\right)Q_1(Z,T)+v\prime \left(s\right)Q_2(Z,T)\right]\mathit{ds}\\ +\underset{\mathrm{\Gamma }}{\int }\left[g_1^{\prime }\left(s\right)Q_1(Z,T)+g_2^{\prime }\left(s\right)Q_2(Z,T)\right]\mathit{ds}=Q\left(Z\right),Z\in L\cup \mathrm{\Gamma },\end{array}$$ E10

where $Q_j\left(Z,T\right)$ are stress vectors $q_L$ at point $Z$ with coordinates $\left(x,y\right)$ . $L\cup \mathrm{\Gamma }$ , the stress vector is determined by the formula (4) accordingly through complex potentials ${\mathrm{\Phi }}_j\left(z_1,t_1\right),{\mathrm{\Psi }}_j\left(z_2,t_2\right),j=1,2$ ; $T$ is a point with coordinates $\left(\xi ,\eta \right)$ , which belongs to the contour $L\cup \mathrm{\Gamma }$ ; $Q\left(Z\right)=Q_L\left(Z\right)-Q_S\left(Z\right)-Q_{\mathrm{\Delta }}\left(Z\right)$ with $Z\in L$ and $Q\left(Z\right)=Q_T\left(Z\right)-Q_S\left(Z\right)-Q_{\mathrm{\Delta }}\left(Z\right)$ with $Z\in \mathrm{\Gamma }$ ; $Q_L=X_L+{\mathit{iY}}_L$ ; and $Q_m=X_m+{\mathit{iY}}_m$ , where $X_m=2Re\left[s_1{\mathrm{\Phi }}_m\left(z_1\right)z_1^{\prime }+s_2{\mathrm{\Psi }}_m\left(z_2\right)z_2^{\prime }\right],Y_m=-2Re\left[{\mathrm{\Phi }}_m\left(z_1\right)z_1^{\prime }+{\mathrm{\Psi }}_m\left(z_2\right)z_2^{\prime }\right]$ and

$$m=S,\mathrm{\Delta }.$$

Using the results [12], we obtained that the unknown functions $u\prime ,v\prime$ at the boundary of each of the holes $L_j,j=0,1,\dots ,J$ in representation (6) are defined up to a summand $\tilde{u}\prime =-{\omega }_j\mathit{dy}/\mathit{ds},\tilde{v}\prime ={\omega }_j\mathit{dx}/\mathit{ds}$ , where ${\omega }_j$ are arbitrary constants. At numerical solution of the problem, the constants ${\omega }_j,j=0,\dots ,J$ are to be necessarily fixed. In addition, to ensure the displacement continuity condition, it is necessary to impose the following conditions on unknown functions:

$$\begin{array}{c}\underset{L_j}{\int }u\prime \mathit{ds}=0,\underset{L_j}{\int }v\prime \mathit{ds}=0,j=0,\dots ,J;.\\ \underset{{\mathrm{\Gamma }}_j}{\int }{g}_1^{\prime }\mathit{ds}=0,\underset{{\mathrm{\Gamma }}_j}{\int }{g}_2^{\prime }\mathit{ds}=0,j=1,\dots ,K\end{array}$$ E11

Let us consider a problem-solving equation (10) for the case of one hole and a crack. Let us assume that the contour on which the crack is placed is described parametrically in the form $x={\alpha }_{\mathrm{\Gamma }}\left(\tau \right)$ , $y={\beta }_{\mathrm{\Gamma }}\left(\tau \right),-1\le \tau \le 1$ , and the equation of the boundary hole is described in the form $x={\alpha }_L\left(\theta \right)$ , $y={\beta }_L\left(\theta \right),0\le \theta <2\pi$ .

Let us assume the representation for the displacement discontinuity at the cracks:

$$g_1^{\prime }s\prime =\frac{{\mathit{dg}}_1}{\mathit{d\tau }}=\frac{U_{\mathrm{\Gamma }}\left(\tau \right)}{\sqrt{1-{\tau }^2}},g_2^{\prime }s\prime =\frac{{\mathit{dg}}_2}{\mathit{d\tau }}=\frac{V_{\mathrm{\Gamma }}\left(\tau \right)}{\sqrt{1-{\tau }^2}}.$$

Let us replace the integrals with Lobatto-type quadrature formulas [15], and the integrals at the boundaries of the holes replaced by the quadrature of a rectangle, which, for periodic functions, are Gauss quadrature-type formulas [12]. Then we obtain the system of equations:

$$H\sum _{k=1}^{N_O}\left(q_{\mathit{\nu k}}^{\left(1\right)}{U}_k^L+q_{\mathit{\nu k}}^{\left(2\right)}{V}_k^L\right)+\sum _{m=1}^{N_{\mathrm{\Gamma }}}{C}_m\left(q_{\mathit{\nu m}}^{\left(1\right)}{U}_m^{\mathrm{\Gamma }}+q_{\mathit{\nu m}}^{\left(2\right)}{V}_m^{\mathrm{\Gamma }}\right)=q_{\nu },\nu =1,\dots ,N_O+N_{\mathrm{\Gamma }}-1,$$ E12

where $q_{\mathit{\nu k}}^{\left(j\right)}=Q_j\left(Z_{\nu },T_k^L\right),p_{\mathit{\nu m}}^{\left(j\right)}=Q_j\left(Z_{\nu },T_m^{\mathrm{\Gamma }}\right),q_{\nu }=Q_L\left(Z_{\nu }\right)-Q_S\left(Z_{\nu }\right)-Q_{\mathrm{\Delta }}\left(Z_{\nu }\right)$ ,

$U_k^L=u\prime \left(x_k^L,y_k^L\right)s_k^{\prime },V_k^L=v\prime \left(x_k^L,y_k^L\right)s_k^{\prime }$ , $U_m^{\mathrm{\Gamma }}=g_1^{\prime }\left(x_m^{\mathrm{\Gamma }},y_m^{\mathrm{\Gamma }}\right)s_m^{\prime },V_m^{\mathrm{\Gamma }}=g_2^{\prime }\left(x_m^{\mathrm{\Gamma }},y_m^{\mathrm{\Gamma }}\right)s_m^{\prime }$ ,

$x_k^L={\alpha }_L\left({\theta }_k\right),y_k^L={\beta }_L\left({\theta }_k\right),{\theta }_k=\mathit{Hk},{\tilde{\theta }}_n={\theta }_n-H/2,H=2\pi /N_O$ ,

$x_m^T={\alpha }_{\mathrm{\Gamma }}\left({\tau }_m\right),y_m^T={\beta }_{\mathrm{\Gamma }}\left({\tau }_m\right),{\tau }_m=-cos\left({\pi }_N\left(m-1\right)\right),m=1,\dots ,N_{\mathrm{\Gamma }};C_m={\pi }_N\mathrm{at}m\ne 1$ and $m\ne N_{\mathrm{\Gamma }};C_1=C_{N_{\mathrm{\Gamma }}}=0,5{\pi }_N$ ; ${\pi }_N=\frac{\pi }{N_{\mathrm{\Gamma }}-1};T_k^L$ is a point with coordinates $\left(x_k^L,y_k^L\right)$ , $T_m^{\mathrm{\Gamma }}$ is a point with coordinates $\left(x_m^{\mathrm{\Gamma }},y_m^{\mathrm{\Gamma }}\right)$ , $Z_{\nu }$ is a point with coordinates $\left({\tilde{x}}_{\nu }^L,{\tilde{y}}_{\nu }^L\right)$ with $1\le \nu \le N_L$ , and is a point with coordinates $\left({\tilde{x}}_{\nu -N_L}^{\mathrm{\Gamma }},{\tilde{y}}_{\nu -N_L}^{\mathrm{\Gamma }}\right)$ with $N_L<\nu \le N_L+N_{\mathrm{\Gamma }}-1$ , where ${\tilde{x}}_k^L={\alpha }_L\left({\tilde{\theta }}_k\right),{\tilde{y}}_k^L={\beta }_L\left({\tilde{\theta }}_k\right),{\tilde{\theta }}_n={\theta }_n-H/2,{\tilde{x}}_m^{\mathrm{\Gamma }}={\alpha }_{\mathrm{\Gamma }}\left({\tilde{\tau }}_m\right),{\tilde{y}}_m^{\mathrm{\Gamma }}={\beta }_{\mathrm{\Gamma }}\left({\tilde{\tau }}_m\right),{\tilde{\tau }}_m=-cos\left[{\pi }_N\left(m-0,5\right)\right]$ .

We obtain the additional equation of system (12) from condition (11)

$$\sum _{k=1}^{N_{\mathrm{\Gamma }}}{C}_k\left(U_k^{\mathrm{\Gamma }}+{\mathit{iV}}_k^{\mathrm{\Gamma }}\right)=0.$$ E13

Analogously to [12], we should remove three equations with $1\le \nu \le N_L$ and add the following three equations to the received incomplete system:

$$\sum _{k=1}^N{U}_k^L=0,\sum _{k=1}^N{V}_k^L=0,U_m^L=0,1\le m\le N$$ E14

The first two equations follow from the displacements continuity conditions (9). The last equation is obtained when fixing an arbitrary constant $\omega$ .

The system of Eqs. (12)–(14) is generalized in the case of hole and crack system in the same way as it was done in [12].

3. Stresses in circular samples with cracks under the action of concentrated forces

Let us consider the circular composite plate with radius a, which is weakened by a central crack with the half-length L. The plate is stretched by the concentrated forces $\pm P$ applied at points $\left(0,\pm y_0\right)$ . Destruction of the plate is happening when the stress intensity factors (SIFs) reach a certain limit value. Therefore, when we calculated the limiting loads, we considered the SIF which explicitly takes into account the length of the cracks. Due to this, we performed calculations of the relative SIF $K_a=\frac{K_I\sqrt{a}}{P\sqrt{\pi }}$ with different relative distances $\alpha =y_0/a$ , depending on the half-length of the crack, which is divided into a.

Calculations are made for the composite plates with elastic constants shown in Table 1 .

Material	$E_1$ (GPa)	$E_2$ (GPa)	$G_{12}$ (GPa)	${\mathit{\nu }}_{21}$	${\mathit{\nu }}_{12}$
LU	10.8	96.0	2.61	0.21	0.024
EF	21	32.8	5.7	0.21	0.134

Table 1.

Elastic constants of LU and EF materials.

The results of the calculations for the plate made of an EF material (with a small degree of anisotropy) with the maximum stiffness in the direction of the OX axis are shown in Figure 3 .

Figure 3.

Relative SIF for a circular plate made of an EF material: a direction with maximum stiffness is parallel to the OX axis.

In general, the character of the distribution for an EF material is not significantly different from that of an isotropic material. It is necessary to increase monotonically the load for a stable growth of the crack when the distances of forces to a crack are smaller than $0,2a$ . With $y_0/a=0,3$ , the stable growth of cracks (without jumping) will occur at $L/a>0,2$ ; with $y_0/a=0,4$ if $L/a>0,25$ ; and with $y_0/a=0,5$ if $L/a>0,35$ . At greater distances to forces, after reaching the corresponding level of values of traction, the circle fractures. For a case where the crack is perpendicular to the direction with the maximum stiffness of the material, the SIF is slightly increasing, especially at greater distances to forces.

For a weakly anisotropic material, the incline of the crack to the main axis of orthotropy had little effect on the SIF $K_I$ , and the SIF $K_{\mathit{II}}$ is practically absent.

The calculations have shown that for the case of placing the crack in parallel to the direction with the maximum stiffness of the material, the above set of specifics of the SIF remain unchanged for substantially anisotropic LU-1 material. When the crack is placed perpendicular to the direction with maximum stiffness, stable crack growth occurs of small cracks ( $L/a\le 0,1$ ) when the relative distance to forces is $y_0/a<0,2$ . Moreover, for these cases, the fracture is spasmodic. In all other cases, the SIF increases monotonically with increasing crack length. Hence, the plate fractures completely after reaching forces of critical value.

Testing the developed algorithm is conducted for the case of isotropic plate with $\nu =0,4488$ and $y_0/a=0,15;0,16;0,18;0,2$ . The calculation results of the relative SIF $K_a=\frac{K_I\sqrt{a}}{P}$ are shown in Figure 4 .

Figure 4.

Relative SIF for isotropic material.

On the right are shown figures from the book [10]. Such calculations were also performed for the same material, and here the corresponding relative SIF is represented by dashed lines. It is seen that the results obtained by different methods coincide.

4. Determination of working intervals of crack lengths at circular samples

Two types of samples are used in experimental studies of crack fracture resistance [10]. The first is a sample for which the SIF grows monotonically with the growth of the crack. In the second type, the range of cracks’ lengths is selected in such a way that the SIF $K_I$ is practically constant. Hence, it is varied in this range from the mean value to the small value (∼2–4%). This range of crack lengths is called working. The samples of the second type are particularly suitable for conducting experimental studies including a wide range of problems in the area of destruction. In particular, with the constant force factor (the SIF is constant) in such samples, the possibility of an effective study of the rate of growth of fatigue cracks with cyclic loads, the study of crack fracture resistance depending on the influence of working environment, etc. arises.

Based on the studies in the literature for isotropic material, it has been established [10] that the range of working lengths is most favorable with $\alpha =0,18$ , although allowed, and $\alpha =0,16;0,2$ , where $\alpha =y_0/a$ . With $0<\alpha <0,2237$ the SIF $K_I$ increases from zero to a certain maximum (depending on $\alpha$ ), then falls to a minimum, and then increases monotonously. With $\alpha >0,2237$ the SIF increases monotonically with the increasing length of the crack. We note that such conclusions are valid for material with a Poisson ratio $\nu =0,4488$ .

Since at big lengths of cracks $L/a>0,5$ the SIF depends little on $\alpha$ , then for the first type of samples $\alpha =0,65$ is taken.

Based on the obtained results, let us perform a similar study of samples of two types of composite materials.

For samples of the second type, we perform calculations only for small ratios (at $\alpha \sim 0,2$ ). The results of calculations of the SIF in a circular isotropic sample with $\nu =1/3$ are shown in Figure 5 . Here and further is assumed that the relative value of the SIF is equal to $K_a=\frac{K_I\sqrt{a}}{P\sqrt{\pi }}$ and the parameter value $\alpha$ is indicated near the curves.

Figure 5.

Relative SIF for isotropic material.

Figure 5 shows that the range with few changed SIF is necessary to determine in the vicinity of the lengths of the cracks with $L/a\sim 0,25$ .

Similar results for the plate made of an EF material are shown in Figures 6 and 7 . Here two cases are considered: the crack is parallel or perpendicular to the direction in which the stiffness of the material is maximal.

Figure 6.

Relative SIF for an EF material: the crack is parallel to the direction of greater stiffness of the material.

Figure 7.

Relative SIF for an EF material: the crack is perpendicular to the direction of greater stiffness of the material.

From the given data, it is seen that the SIF for an EF material is bigger in the case when the crack is perpendicular to the direction of maximum stiffness of the material.

The following conclusions are made based on the data of the calculations: for the crack that is parallel to the direction of maximum stiffness of the material, the minimum deviations from the constant SIF (with an error of not more than 2%) are achieved on the ranges of lengths of cracks with relative dimensions with $\mathrm{\Delta }<0,21$ with $\alpha =0,24-0,26$ , where $\mathrm{\Delta }=\left(L_2-L_1\right)/a$ . The biggest range ( $\mathrm{\Delta }=0,27$ ) with SIF values close to constants with an error $2,1\%$ is achieved with $\alpha =0,24$ . For the crack that is perpendicular to the direction of maximum stiffness of the material, the range with SIF values close to constants is reduced. Moreover, the distance of force application needs to be increased. The biggest range ( $\mathrm{\Delta }=0,27$ ) with SIF values close to constants with an error $3,2\%$ is achieved with $\alpha =0,24$ .

Let us consider the case of samples of the first type 1. For them, the forces are selected that are distant from the crack. For isotropic materials, as a rule, $\alpha =0,65$ is taken. The above-obtained results of the calculations show that the same distance can be chosen for the composite materials with a crack parallel to the principal axes of the orthotropy.

To compare the effect of a sample shape, similar calculations are made for a square plate with a crack. The results of calculations for such an isotropic sample that are similar to the results of calculations for a circle are shown in Figure 8 .

Figure 8.

Relative SIF for a square isotropic sample.

The conclusion is made based on the comparison of Figures 5 and 8 that with small distances of forces from cracks, the shape of the sample has little effect on the SIF.

Similar results of calculations for LU material are shown in Figure 9 for a horizontal crack and in Figure 10 for a diagonal crack.

Figure 9.

LU material, square sample, horizontal crack.

Figure 10.

LU material, square sample, diagonal crack.

5. Determination of the SSS of samples under the action of the tractions applied to the hole’s boundary

Let us apply the developed algorithm to study a square plate with a half-side a; weakened by a central crack with a half-length L ( Figure 11 ), the edges are not loaded.

Figure 11.

Scheme of the sample.

Two identical circular holes of radius R, the centers of which are located at points $\left(0,\pm c\right)$ , are created for stretching in a plate. It was assumed that the load was applied to the boundary of the circular holes. Using [10], we accept that the forces act normally on the domain $\left|\theta -{\theta }_c\right|<\gamma$ and are given on it in the form

$$p=P\frac{{\left(\theta -{\theta }_c\right)}^2}{4R\left(sin\gamma -\gamma cos\gamma \right)},$$ E15

where $\theta$ is angle coordinate on each of the holes, ${\theta }_c$ is angle coordinate of the middle of the domain, and Р is the principal vector applied to the domain of forces, which is directed from the center of the hole at an angle ${\theta }_c$ .

At first, for the purpose of testing the algorithm, the calculations are performed for the case of a localized load at $\gamma =\pi /32$ and $\alpha =c/a=0,2;0,4;0,6$ , with ${\theta }_c=\pi /2$ for the upper hole and ${\theta }_c=-\pi /2$ for the lower hole (thus the stretching of the plate in the direction of the OY axis is considered). The relative SIF $Y=\frac{K_I}{P}\sqrt{\mathit{\pi L}}$ for an isotropic material with a Poisson ratio $\nu =1/3$ at different crack lengths at $R/a=0,1$ is calculated and given in Table 2 .

$L/a$	$Y$	$Y_S$	$Y$	$Y_S$	$Y$	$Y_S$
	$\mathit{\alpha }$ = 0.2		$\mathit{\alpha }$ = 0.4		$\mathit{\alpha }$ = 0.6
0.1	0.175	0.177	0.248	0.247	0.271	0.270
0.2	0.721	1.720	0.583	0.581	0.568	0.565
0.3	1.164	1.159	0.958	0.953	0.894	0.888
0.4	1.503	1.495	1.323	1.314	1.237	1.227
0.5	1.819	1.805	1.675	1.663	1.588	1.575
0.6	2.156	2.131	2.034	2.015	1.954	1.935
0.7	2.552	2.505	2.438	2.406	2.367	2.334
0.8	3.112	2.995	2.999	2.917	2.937	2.854
0.9	4.279	—	4.178	—	4.131	—

Table 2.

Relative SIFs for a square sample, isotropy.

The values of relative SIF $Y_S$ obtained by another method in [10] for the case of stretching by concentrated forces (i.e., with $\gamma \to 0$ ) are given in the same table. As you can see, the results for these cases were close, with the exception of the values $L/a=0,8$ . Some differences in them are due to different localizations of applied tractions.

The results of the calculations of SIF are given in Table 3 for the case $R\to 0,\gamma =0$ , that is, the case is considered when forces are applied at the center of the holes and more smoothly applied efforts (at $\gamma =\pi /8$ ).

	$\mathit{\gamma }\mathbf{=}\mathit{\pi }/8$			$R\mathbf{\to }0,\mathit{\gamma }=0$			$\mathit{\gamma }=\mathit{\pi }/8$ (circle)
$L/a$	$\mathit{\alpha }$ = 0.2	0.4	0.6	$\mathit{\alpha }$ = 0.2	0.4	0.6	$\mathit{\alpha }$ = 0.2	0.4	0.6
0.1	0.186	0.254	0.275	0.727	0.462	0.366	0.2069	0.2637	0.2767
0.2	0.732	0.593	0.576	1.073	0.851	0.717	0.7563	0.6096	0.5786
0.3	1.172	0.969	0.905	1.278	1.167	1.049	1.1914	0.9910	0.9101
0.4	1.510	1.335	1.250	1.489	1.452	1.370	1.5415	1.3678	1.2655
0.5	1.827	1.689	1.605	1.742	1.744	1.694	1.8964	1.7516	1.6483
0.6	2.167	2.050	1.974	2.048	2.067	2.039	2.3113	2.1770	2.0811
0.7	2.569	2.459	2.391	2.436	2.459	2.440	2.8332	2.6963	2.6106
0.8	3.139	3.029	2.971	3.008	3.027	3.011	3.5616	3.4216	3.3492
0.9	4.322	4.225	4.181	4.214	4.225	4.209	4.8896	4.7662	4.7132

Table 3.

Relative SIFs for a square and circular sample, isotropy.

The results for a circular sample with $\gamma =\pi /8$ are shown ( Figure 12 ) in the same table.

Figure 12.

Circular sample with holes and a crack.

The following conclusions are made based on Tables 2 and 3 : SIF does not differ significantly in the case of distributed loads with different degrees of localization, SIF increases somewhat with the growth of the domain of action of tractions, and SIFs are bigger at small crack lengths at point action of tractions (at $R\to 0$ ) and at all lengths with $\alpha >0,6$ . SIFs in a circular sample are bigger than in a square one under the same load conditions.

Similar results for a square sample made from a LU material are given in Table 4 . Here, the relative SIFs $K_a$ in which the crack length is explicitly taken into account and the case where the crack is parallel to the direction of bigger stiffness (data on the left) and is perpendicular to it (data on the right) are given.

	LU (Ох)			LU (Оу)
$L/a$	$\mathit{\alpha }$ = 0.2	0.4	0.6	$\alpha$ = 0.2	0.4	0.6
0.1	0.2338	0.2633	0.2764	0.4182	0.478	0.538
0.2	0.5187	0.4276	0.4059	1.0504	1.0073	0.9911
0.3	0.5952	0.5259	0.4978	1.1416	1.1306	1.1327
0.4	0.6202	0.5791	0.5595	1.2268	1.2138	1.2163
0.5	0.6409	0.6158	0.6046	1.33	1.3094	1.31
0.6	0.6711	0.6539	0.6479	1.4615	1.4313	1.4304
0.7	0.7252	0.7117	0.7085	1.6257	1.5829	1.5809
0.8	0.8392	0.8284	0.8267	1.7876	1.7411	1.7386
0.9	1.1482	1.1419	1.1409	1.9405	1.8846	1.8831

Table 4.

Relative SIFs $K_a$ for square sample, LU material.

Similar results for a circular sample for the same material are given in Table 5 .

	LU (Ох)			LU (Оу)
$L/a$	$\mathit{\alpha }$ = 0.2	0.4	0.6	$\mathit{\alpha }$ = 0.2	0.4	0.6
0.1	0.2363	0.2679	0.2806	0.4192	0.4799	0.5407
0.2	0.5269	0.4369	0.4138	1.0594	1.0163	0.9995
0.3	0.6108	0.5418	0.5112	1.1602	1.1477	1.1466
0.4	0.6467	0.605	0.5814	1.2588	1.2422	1.2389
0.5	0.683	0.6558	0.6396	1.3795	1.3531	1.3459
0.6	0.7341	0.7131	0.702	1.5354	1.4974	1.4888
0.7	0.8131	0.7948	0.7879	1.7482	1.6965	1.6905
0.8	0.9524	0.9376	0.9342	2.0507	1.9973	1.9977
0.9	1.2791	1.2726	1.2704	2.592	2.528	2.5346

Table 5.

Relative SIFs for a circular sample, LU material.

6. Determination of the SSS of samples loaded with concentrated forces at the boundary (compression test)

The methods of studied crack fracture resistance based on sample compression in experimental practice are widely used.

The direct application of the abovementioned variant of the method of integral equations for the case when the concentrated forces act on the boundary of the domain is associated with significant errors, because unknown functions in the vicinity of the points of application of forces have a singularity. Due to this, it is necessary to separate a singular part in the solution for a more precise solution of this type of task.

6.1 Determining a singular part of the solution of this problem

Let us consider a point of the plate boundary $z_0$ , in which the concentrated force (X, Y) is applied. A singular part of the solution (of Lekhnitskii potentials) will be the same as in the half-plane, whose boundary is tangent to the plate at the point of action of the concentrated force. Let us mark this angle through $\phi$ and potentials for the half-plane through ${\mathrm{\Phi }}_0\left(z_1\right),{\mathrm{\Psi }}_0\left(z_2\right)$ .

Let us consider at first a half-plane $y<0$ , which is loaded with force (X, Y) at an arbitrary point $x_0$ on the boundary. The Lekhnitskii potentials for this half-plane will be [17]

$${\mathrm{\Phi }}_0=\frac{A\prime }{z_1-x_0},{\mathrm{\Psi }}_0=\frac{B\prime }{z_2-x_0},$$

where

$$A\prime =-\frac{X+s_{2}Y}{2\mathit{\pi i}\left(s_1-s_2\right)},B\prime =\frac{X+s_{1}Y}{2\mathit{\pi i}\left(s_1-s_2\right)}.$$

It can be shown [17] that the half-plane also corresponds to the ${\mathrm{\Phi }}_0,{\mathrm{\Psi }}_0$ potentials, whose boundary passes through the point $\left(x_0,0\right)$ and is inclined at an arbitrary angle under the action of the same force. Let us consider a semicircle with the center at point $\left(x_0,0\right)$ of radius ${\rho }_0$ , which belongs to the half-plane. It is easy to show that the principal vector of all forces applied to the arc of the semicircle is equal to $\left(-X,-Y\right)$ . This proves that the case of loading by the concentrated force of the half-plane corresponds to the potentials ${\mathrm{\Phi }}_0,{\mathrm{\Psi }}_0$ .

Let us now consider a bounded plate occupying the domain D. The self-balanced concentrated forces $\left(X_j,Y_j\right),\left(j=1,\dots ,J\right).$ are applied to the boundary of this domain at points $z_j=x_j+{\mathit{iy}}_j$ . Let us represent the complex potentials in the form

$$\mathrm{\Phi }\left(z_1\right)={\mathrm{\Phi }}_0\left(z_1\right)+{\mathrm{\Phi }}_{\mathrm{\Delta }}\left(z_1\right),\mathrm{\Psi }\left(z_1\right)={\mathrm{\Psi }}_0\left(z_1\right)+{\mathrm{\Psi }}_{\mathrm{\Delta }}\left(z_1\right),$$ E16

$${\mathrm{\Phi }}_0\left(z_1\right)=\sum _{j=1}^J\frac{A_j^{\prime }}{z_1-z_{1j}},{\mathrm{\Psi }}_0\left(z_2\right)=\sum _{j=1}^J\frac{B_j^{\prime }}{z_2-z_{2j}},$$

where $z_{1j}=x_j+s_1{y}_j,z_{2j}=x_j+s_2{y}_j$ . Here the coefficients $A_j^{\prime },B_j^{\prime }$ are determined based on expressions for $A\prime ,B\prime$ by the substitution of X and Y on $X_j$ and $Y_j$ , respectively. By substituting formulas (16) into boundary conditions, we obtain the boundary problem for obtaining the introduced complex potentials at $\left(x,y\right)\in L$ :

$$\begin{array}{l}\left(1+{\mathit{\text{is}}}_1\right)z_1^{\prime }{\mathrm{\Phi }}_{\mathrm{\Delta }}\left(z_1\right)+\left(1+i\overline{s_1}\right)\overline{z_1^{\prime }}\overline{{\mathrm{\Phi }}_{\mathrm{\Delta }}\left(z_1\right)}+\\ +\left(1+{\mathit{\text{is}}}_2\right)z_2^{\prime }{\mathrm{\Psi }}_{\mathrm{\Delta }}\left(z_2\right)+\left(1+i\overline{s_2}\right)\overline{z_2^{\prime }{\mathrm{\Psi }}_{\mathrm{\Delta }}\left(z_2\right)}=-i{\left(X+\mathit{iY}\right)}_0,\end{array}$$ E17

where L is the boundary of domain D:

$$i{\left(X+\mathit{iY}\right)}_0=\left(1+{\mathit{\text{is}}}_1\right)z_1^{\prime }{\mathrm{\Phi }}_0\left(z_1\right)+\left(1+i\overline{s_1}\right)\overline{z_1^{\prime }}\overline{{\mathrm{\Phi }}_0\left(z_1\right)}+\left(1+{\mathit{\text{is}}}_2\right)z_2^{\prime }{\mathrm{\Psi }}_0\left(z_2\right)+\left(1+i\overline{s_2}\right)\overline{z_2^{\prime }{\mathrm{\Psi }}_0\left(z_2\right)}.$$

It is easy to show that the right-hand side of formula (17) is a continuous and limited function, and therefore the introduced complex potentials with an index $\mathrm{\Delta }$ are continuous and limited in the vicinity of the points of application of forces. In this regard, the above-developed numerical algorithm based on BIEM can be used to determine these potentials.

6.2 Calculation of the SIF for a rectangular sample with compression

Let us consider a square plate with half-side a, which contains a diagonal central vertical crack with half-length $L$ . The crack fracture resistance of such a sample is determined based on compression by force R applied in vertical direction. The tips of the angles, in the vicinity of which forces are applied, can be cut off. Due to this, a sample is considered whose tops have coordinates

$$z_1=c,z_2=\mathit{ic}+z_p,z_3=\mathit{ic}+z_m,z_4,z_5=\overline{z_3}{z}_6=\overline{z_2},$$

where in $c=\sqrt{2}a,z_p=h-\mathit{ih},z_m=-h-\mathit{ih},$ h is the height of the cut triangle and а is the half-side of the square.

The calculations are performed at $h=c/8$ ; moreover, all tops are rounded by the arcs of the circle of the radius $a/10$ (the shape of the sample— Figure 13 ).

Figure 13.

Sample view.

Calculation of the relative SIF $K_a=\frac{K_I\sqrt{a}}{P\sqrt{\pi }}$ for isotopic material and composites of EF and LU with different directions of the orthotropic axis are given in Table 6 . The angle between the crack and the direction with the maximum stiffness of the material is indicated in brackets.

$L/a$	Isotr.	EF (90°)	EF (0°)	LU (90°)	LU (0°)
0.1	0.0824	0.0579	0.0727	0.0226	0.0743
0.2	0.1203	0.085	0.1073	0.0332	0.1182
0.3	0.1557	0.1106	0.1407	0.0434	0.1663
0.4	0.1927	0.1384	0.1776	0.0547	0.2217
0.5	0.2354	0.1709	0.2208	0.068	0.2865
0.6	0.2889	0.2107	0.2737	0.0846	0.3643
0.7	0.3544	0.2609	0.3406	0.1059	0.4609
0.8	0.4418	0.3273	0.4301	0.1343	0.5872

Table 6.

Relative SIF $K_a$ when compressing a sample with a diagonal crack.

The table shows a significant effect on the SIF of the placement of the crack relative to the axis with the maximum stiffness of the material. In particular, for cracks parallel to the maximum stiffness direction, the SIFs appeared to be significantly larger than those returned for 90°. The difference between the SIFs for these two directions is increasing for a substantially anisotropic LU material.

Table 7 shows the results of calculations for the case of stretching the same sample with a horizontal crack ( Figure 14 ).

$L/a$	Isotr.	EF (0°)	EF (90°)	LU (0°)	LU (90°)
0.1	0.2585	0.2886	0.3499	0.2956	0.7572
0.2	0.3696	0.4065	0.4874	0.4114	0.9247
0.3	0.4617	0.4963	0.5884	0.4936	1.0104
0.4	0.5494	0.5752	0.6766	0.5594	1.1008
0.5	0.641	0.6534	0.7668	0.6185	1.2159
0.6	0.7434	0.7397	0.8705	0.6789	1.364
0.7	0.8641	0.842	0.9973	0.7472	1.5542
0.8	1.0129	0.9704	1.1585	0.8319	1.8021

Table 7.

Relative SIFs when stretching a sample with a diagonal crack.

Figure 14.

Sample view.

Based on the comparison of data from Tables 6 and 7 , it follows that, unlike the case of compression, the SIF with stretching is larger for cracks that are perpendicular to the direction with maximum stiffness of the material.

7. Conclusions

An algorithm for calculating stresses at cracks in bounded plate with holes of various shapes due to concentrated forces or distributed forces at its boundary has developed. The solution of integral equations is performed by quadrature Gauss-type formulas for regular and singular integrals.

The research of stresses at cracks in the samples which are used in experimental studies of crack fracture resistance was performed.

The calculation of the stresses at cracks in samples of various forms is performed, in which ones’ experimental research are performed. To study the crack fracture resistance of composite samples, the optimal distances from the central crack to the forces at which the SIF increases monotonously with increasing crack length are determined. In particular, for square samples with a half-side a, forces should be placed at a distance of 0.6a–0.7a from the crack. For the experimental study of the growth rate of fatigue cracks, there are definite ranges of lengths of cracks for which SIFs are practically constant values. At the same time, the distances are determined at which it is expedient to apply forces. The problem for studying samples with cracks with compression is considered.

List of symbols with explanations

s is an arc coordinate.

$\nu$ is a Poisson ratio.

${\epsilon }_x,{\epsilon }_y,{\gamma }_{\mathit{xy}}$ are strains.

${\sigma }_x,{\sigma }_y,{\tau }_{\mathit{xy}}$ are stresses.

${\alpha }_{\mathit{ij}}$ are elastic compliances which are included in the Hooke’s law.

$\mathrm{\Phi }\left(z_1\right),\mathrm{\Psi }\left(z_2\right)$ are Lekhnitskii complex potentials.

${\mathrm{\Phi }}_S,{\mathrm{\Psi }}_S$ are the potentials which correspond to the concentrated forces.

$u\prime ,v\prime$ are the values of the derivatives of the displacements with respect to the arc coordinate at the boundary of the plate and holes.

$u^{\pm },v^{\pm }$ are limit values of displacements in the approach to the section at the left and the right relative to the selected direction.

$\overrightarrow{q_{\mathrm{\Gamma }}}$ is the stress vector $q_{\mathrm{\Gamma }}\left(z\right)=X_{\mathrm{\Gamma }}+{\mathit{iY}}_{\mathrm{\Gamma }}$ at path $\mathrm{\Gamma }$ .

$K_I,K_{\mathit{II}}$ are the stress intensity factors (SIFs).

$K_a$ , $Y$ , $Y_S$ are the relative SIFs.

Р is the principal vector.