Numerical modeling of fracture processes of bodies with stress concentrators under conditions of proportional loading, taking into consideration the statistical distribution of ultimate strength and partial loss of load bearing capacity

Introduction

Methodology

1. Consideration of the failure processes, taking into account the heterogeneity of the distribution of strength properties of structural elements, during which variating structure's macrolevel behavior takes place (from elastic-brittle to quasi-plastic with realization of the postcritical stage) depending on the value of the ultimate strength dispersion;
2. Three characteristic types of damage accumulation were identified: localized (with the development of a macrodefect), dispersed (with the accumulation of damage in the entire volume), mixed (with the growth of a macrodefect in the most weakened areas). To assess the type of damage accumulation, based on the results of solving boundary value problems of the theory of elasticity, a new approach has been proposed, which consists in the analysis of overload coefficients;
3. The geometry of the stress concentrator has a significant impact on the behavior of the structure at the macro level and the strength reserve, however, there is a threshold value of the variation coefficient of ultimate strength, after which the concentrator ceases to affect the behavior of the structure;
4. Reducing the size of the finite elements leads to a decrease in the load bearing capacity of the body, while there is no convergence of the numerical solution. In this regard, it is suggested that the characteristic size of finite elements should be physically justified and selected for different materials on the basis of experimental studies.
   Nevertheless, several disadvantages and limitations were noted:
5. The works consider only a uniform distribution, poorly suitable for describing a random distribution of material's mechanical properties, and a two-parameter Weibull distribution, using which close to zero values of ultimate strength were generated;
6. In the work, at the realization of the fracture criterion, all stiffness properties of the final element were reduced, which does not quite correspond to the real behavior of elastic-brittle materials, for which the bearing capacity disappears, when stretched across the crack plane, and remains, when deformed along the crack plane (i.e. significant anisotropy of properties appears). This drawback can be significant when modeling fracture processes under biaxial loading conditions;
7. Only uniaxial loading was considered.

1. Material. Elastic-brittle material without plastic deformation, viscoelastic behavior or friction between damaged parts of the structure is considered. The body is divided into subregions, the material inside each subregion is isotropic and homogeneous, the elastic moduli at all points are the same, the ultimate strength is distributed statistically according to the three-parameter Weibull law. Partial loss of the bearing capacity of the material within the subregion occurs when the first principal stress reaches ultimate strength. In this case, in the subregion, the intact isotropic material is replaced by an orthotropic one, the anisotropy axes of which are oriented in accordance with the directions of the principal stresses at the moment of damage (the propagation of the crack through the subregion is simulated). Elastic properties of the orthotropic material ensure absence of resistance during the deformation across the simulated crack, but show resistance during the deformation along it. In case, when in the material, that has partially lost its bearing capacity, the normal stress (oriented along the simulated crack) again exceeds the ultimate strength, the material is replaced with an isotropic one that has completely lost its bearing capacity. Since the material that has completely lost its bearing capacity is extremely malleable, geometric nonlinearity is considered. The described scheme for the material properties' changing, when the failure criterion is met, is illustrated in Fig. 1. The implementation of this scheme is described in the section "Boundary value problem and its solution algorithm".
2. Body with stress concentrator at biaxial loading. A body with an internal stress concentrator in a plane stress state (plate) is considered. Kinematic boundary conditions are applied to the outer boundaries of the body, while changes of all displacements are made proportionally (i.e. complex loading is not considered). At the boundaries of the stress concentrator, boundary conditions of the second type with zero forces are applied;
3. Features of the numerical realization. When dividing a body into finite elements, it is assumed that one subregion (within which the material is homogeneous) corresponds to one finite element. Also, there is no additional consideration of the effect of the size of the finite element mesh on the modeling results, it is initially built small. These assumptions help to avoid the necessity of building a random structure of the material inside the body and reduce the number of tasks to be solved, while ensuring the possibility of checking the adequacy and feasibility of the developed model's applicability. In accordance with the results of work [5], at each stage of solving the boundary value problem, one of the most overloaded elements is deactivated, the loading step is selected automatically, based on the analysis of stress fields, in case of damage or deactivation of one finite element, the boundary conditions remain unchanged until a stable state is obtained. For each of the variation coefficients, five generations of finite elements' ultimate strength will be considered, for each of the generations, different modes of multiaxial loading will be implemented;

Figure 1: The scheme of the material properties reduction after the fracture criterion fulfillment (plane stress state): $x,y$$x,y$x,yx, y$x,y$ are global coordinate system axes; $x^{\prime },y^{\prime }$$x^{\prime },y^{\prime }$x^('),y^(')x^{\prime}, y^{\prime}$x^{\prime },y^{\prime }$ are the anisotropy axes of the damaged material, corresponding to the orientation of the principal stress before damage.

4. Registered data and their analysis. After each of the calculations, data on the displacements of the external boundaries of the body and the emerging reaction forces are registered, which allows building diagrams of the body's load at the macro level. The number of each damaged/deactivated element is also registered, which allows analyzing the kinetics of the fracture process and assessing the integral damage to the body structure;
5. Assessment of the type of damage accumulation. In accordance with the approach proposed in [7], for the obtained generations of ultimate strength of structural elements, an analysis of overload factors during the elastic deformation of the body will be carried out with an estimation of the number of FEs, in which the threshold value is exceeded, and estimation of the average distance from the concentrator to centers of these FEs. The plotted dependences of these values on the variation coefficient of the ultimate strength will be associated with the realization of the previously identified types of damage accumulation.
   Based on the application of the developed methodology, conclusions will be drawn about the adequacy and feasibility of the developed model's applicability.

Boundary value problem and its solution algorithm

Mathematical formulation

Solution algorithm

1. Initial data input: elastic properties 1 - undamaged isotropic material; elastic properties 2 - damaged orthotropic material, for which it is necessary to take into account the rotation of the anisotropy axes in the direction of the principal stresses at the moment of damage; elastic properties 3 - fractured isotropic material. Generating of the ultimate strength values of FE according to the chosen probability distribution law.
2. Creation output files ${\mathrm{\Phi }}_1$${\mathrm{\Phi }}_1$Phi_(1)\Phi_{1}${\mathrm{\Phi }}_1$ and ${\mathrm{\Phi }}_2$${\mathrm{\Phi }}_2$Phi_(2)\Phi_{2}${\mathrm{\Phi }}_2$ with the solution results.
3. Designing and meshing of the body.
4. Creating of the initial boundary conditions (Eq. (4)), $u_i{}^0$$u_i{}^0$u_(i)^(0)u_{i}{ }^{0}$u_i{}^0$ and $S_i{}^0$$S_i{}^0$S_(i)^(0)S_{i}{ }^{0}$S_i{}^0$ values should be small to prevent the failure criterion fulfillment at the first step.
5. Calculating of the stress-strain state.
6. Calculating of the external load value $P_{(i)}$$P_{(i)}$P_((i))P_{(i)}$P_{(i)}$ (if the loading is multiaxial, several values are calculated). Here ( $i$$i$ii$i$ ) is the step number. Output to ${\mathrm{\Phi }}_1:(i),u_{(i)},P_{(i)}$${\mathrm{\Phi }}_1:(i),u_{(i)},P_{(i)}$Phi_(1):(i),u_((i)),P_((i))\Phi_{1}:(i), u_{(i)}, P_{(i)}${\mathrm{\Phi }}_1:(i),u_{(i)},P_{(i)}$.
7. If $P_{(\mathrm{i})}$$P_{(\mathrm{i})}$P_((i))P_{(\mathrm{i})}$P_{(\mathrm{i})}$ value is less than the critical value $P_{\text{crit}}$$P_{\text{crit }}$P_("crit ")P_{\text {crit }}$P_{\text{crit }}$, the fracture modeling process ends (except for the first step) and the resulting data files ${\mathrm{\Phi }}_1,{\mathrm{\Phi }}_2$${\mathrm{\Phi }}_1,{\mathrm{\Phi }}_2$Phi_(1),Phi_(2)\Phi_{1}, \Phi_{2}${\mathrm{\Phi }}_1,{\mathrm{\Phi }}_2$ are saved.
8. Calculating of the field of the overload coefficient $K$$K$KK$K$ (Eq. (2)), defining of its maximum value $K_{\text{max}}$$K_{\text{max }}$K_("max ")K_{\text {max }}$K_{\text{max }}$ and the number ( $m$$m$mm$m$ ) of the element with $K_{\text{max}}$$K_{\text{max }}$K_("max ")K_{\text {max }}$K_{\text{max }}$.

Figure 2: The flow chart of the boundary value problem solution algorithm.

9. If $K_{\text{max}}\ge 1$$K_{\text{max }}\ge 1$K_("max ") >= 1K_{\text {max }} \geq 1$K_{\text{max }}\ge 1$ and the finite element with number ( $m$$m$mm$m$ ) is assigned the properties of the material with elastic characteristics 1 , then the angle of rotation of the local coordinate system relative to the global one is determined in accordance with the directions of the principal stresses and the elastic properties of an orthotropic material 2 are assigned, $N_{\text{dam}}=1,N_{\text{fract}}=0$$N_{\text{dam }}=1,N_{\text{fract }}=0$N_("dam ")=1,N_("fract ")=0N_{\text {dam }}=1, N_{\text {fract }}=0$N_{\text{dam }}=1,N_{\text{fract }}=0$, $Nu{m}_{\text{dam}}=m,Nu{m}_{\text{fract}}=0$$Nu{m}_{\text{dam }}=m,Nu{m}_{\text{fract }}=0$Num_("dam ")=m,Num_("fract ")=0N u m_{\text {dam }}=m, N u m_{\text {fract }}=0$Nu{m}_{\text{dam }}=m,Nu{m}_{\text{fract }}=0$. If $K_{\text{max}}\ge 1$$K_{\text{max }}\ge 1$K_("max ") >= 1K_{\text {max }} \geq 1$K_{\text{max }}\ge 1$ and the finite element with number (m) is assigned the properties of a material with elastic characteristics 2 , then it is assigned the elastic properties of material $3,N_{\text{dam}}=0,N_{\text{fract}}=1,Nu{m}_{\text{dam}}=0,Nu{m}_{\text{fract}}=m$$3,N_{\text{dam }}=0,N_{\text{fract }}=1,Nu{m}_{\text{dam }}=0,Nu{m}_{\text{fract }}=m$3,N_("dam ")=0,N_("fract ")=1,Num_("dam ")=0,Num_("fract ")=m3, N_{\text {dam }}=0, N_{\text {fract }}=1, N u m_{\text {dam }}=0, N u m_{\text {fract }}=m$3,N_{\text{dam }}=0,N_{\text{fract }}=1,Nu{m}_{\text{dam }}=0,Nu{m}_{\text{fract }}=m$. If $K_{\text{max}}<1$$K_{\text{max }}<1$K_("max ") < 1K_{\text {max }}<1$K_{\text{max }}<1$, then $N_{\text{dam}}=0,N_{\text{fract}}=0,Nu{m}_{\text{dam}}=0,Nu{m}_{\text{fract}}=0$$N_{\text{dam }}=0,N_{\text{fract }}=0,Nu{m}_{\text{dam }}=0,Nu{m}_{\text{fract }}=0$N_("dam ")=0,N_("fract ")=0,Num_("dam ")=0,Num_("fract ")=0N_{\text {dam }}=0, N_{\text {fract }}=0, N u m_{\text {dam }}=0, N u m_{\text {fract }}=0$N_{\text{dam }}=0,N_{\text{fract }}=0,Nu{m}_{\text{dam }}=0,Nu{m}_{\text{fract }}=0$, magnifying of the boundary conditions $1/K_{\text{max}}$$1/K_{\text{max }}$1//K_("max ")1 / K_{\text {max }}$1/K_{\text{max }}$ times.
10. Output to ${\mathrm{\Phi }}_2:(i),N_{\text{dam}},N_{\text{fract}},Nu{m}_{\text{dam}},Nu{m}_{\text{fract.}}$${\mathrm{\Phi }}_2:(i),N_{\text{dam }},N_{\text{fract }},Nu{m}_{\text{dam }},Nu{m}_{\text{fract. }}$Phi_(2):(i),N_("dam "),N_("fract "),Num_("dam "),Num_("fract. ")\Phi_{2}:(i), N_{\text {dam }}, N_{\text {fract }}, N u m_{\text {dam }}, N u m_{\text {fract. }}${\mathrm{\Phi }}_2:(i),N_{\text{dam }},N_{\text{fract }},Nu{m}_{\text{dam }},Nu{m}_{\text{fract. }}$ Going to the step 5.

Model setup

Figure 3: The geometry of the plate and the boundary conditions (a); the proportional biaxial loading modes (b).

RESULTS AND DISCUSSION

Analysis of feasibility of consideration of bearing capacity partial loss during damage

Figure 4: Image of the body with damaged (grey color) and fractured (blue color) FEs, red and black segments represent the orientation of the anisotropy axes of the damaged material (a); image of the imitated cracks (anisotropy axes, along which material properties are unchanged), went through the FEs (b).

Figure 5: The calculated loading diagrams (a), solid plots - $P_x\left(U_x\right)$$P_x\left(U_x\right)$P_(x)(U_(x))P_{x}\left(U_{x}\right)$P_x\left(U_x\right)$, dashed plots - $P_y\left(U_y\right)$$P_y\left(U_y\right)$P_(y)(U_(y))P_{y}\left(U_{y}\right)$P_y\left(U_y\right)$; the corresponding diagrams of growth of the relative number of elements with damaged (solid line) and fractured (dashed lines) material (b): black plots - with partial bearing capacity loss taken into account, green plots - without it; the images of the plate with damaged (red color) and completely fractured (black color) elements (c).

Influence of the coefficient of variation of the ultimate strength distribution on the fracture processes under biaxial loading

Influence of loading mode and variation coefficient of strength properties on the loading diagrams has been considered. Typical loading diagrams are given in Tab. 1. The results demonstrate that an increase in the variation coefficient in all cases leads to a gradual decrease in the bearing capacity of the body and decrease in displacement corresponding to the point of the maximum load, while at the macro level a non-linear deformation begins to appear with the realization of the postcritical stage of body's deformation. However, due to the absence of elements with low strength (as, for example, in the twoparameter Weibull distribution), the loading diagrams are practically linear until the maximum load is reached, pseudoplastic behavior is not observed, unlike in the work [7]. The occurrence of situations, in which a gradual drop in the load on one of the axes leads to a smooth increase in the load on the other of the axes, due to the preservation of resistance in the direction transverse to the direction of macrodefect's growth, was noted. Changing the loading mode with equal variation coefficients of strength properties distribution did not qualitatively affect the type of loading diagrams, however, a nonmonotonic change in the maximum loads along the $x$$x$xx$x$ and $y$$y$yy$y$ axes sustained by the body was observed. In all loading modes, maximum forces $P_x$$P_x$P_(x)P_{x}$P_x$ and $P_y$$P_y$P_(y)P_{y}$P_y$ were achieved simultaneously.

Table 1: The calculated loading diagrams for the various loading modes and coefficients of variation of strength properties.

Figure 6: The strength surface of the plate with the stress concentrator.

Figure 7: The calculated loading diagrams, states 1-3 are shown on the dependence $P_y\left(U_y\right)$$P_y\left(U_y\right)$P_(y)(U_(y))P_{y}\left(U_{y}\right)$P_y\left(U_y\right)$ (a), the dependence of $P_y$$P_y$P_(y)P_{y}$P_y$ on $P_x$$P_x$P_(x)P_{x}$P_x$ (b), the diagrams of the growth of the relative number of damaged (gray line) and completely fractured (black line) elements (c), the kinetics of the damaging process (d) for CV=0.2 under proportional loading with $U_x=U_y$$U_x=U_y$U_(x)=U_(y)U_{x}=U_{y}$U_x=U_y$.

Figure 8: The calculated loading diagrams, states 1-8 are shown on the dependence $P_y\left(U_y\right)$$P_y\left(U_y\right)$P_(y)(U_(y))P_{y}\left(U_{y}\right)$P_y\left(U_y\right)$ (a), the dependence of $P_y$$P_y$P_(y)P_{y}$P_y$ on $P_x$$P_x$P_(x)P_{x}$P_x$ (b), the diagrams of the growth of the relative number of damaged (gray line) and completely fractured (black line) elements (c), the kinetics of the damaging process (d) for CV $=0.5$$=0.5$=0.5=0.5$=0.5$ under proportional loading with $U_x=U_y$$U_x=U_y$U_(x)=U_(y)U_{x}=U_{y}$U_x=U_y$.