NuMeRICAL SIMuLATION OF TWO-DIMeNSIONAL PROBLeMS OF CReeP CRACK GROWTH WITH MATeRIAL DAMAGe CONSIDeRATION

D . B r e s l a v s k y Doctor of Technical Sciences, Professor, Head of Department* E-mail: brdm@kpi.kharkov.ua A . K o z l y u k Postgraduate student* E-mail: alenakozlyk@gmail.com O . T a t a r i n o v a PhD, Associate Professor* E-mail: ok.tatarinova@gmail.com *Department of Computer Modeling of Processes and Systems National Technical University «Kharkiv Polytechnic Institute» Kyrpychova str., 2, Kharkiv, Ukraine, 61002 Надано опис розробленого підходу до чи­ сельного моделювання процесу розповсю­ дження тріщини при повзучості з урахуван­ ням прихованої пошкоджуваності матеріалу. Підхід базується на скінченноелементному моделюванні повзучості й пошкоджуванос­ ті та застосуванні алгоритму перебудови сітки з виключенням зруйнованих елементів. Визначено час та форму області руйнуван­ ня, отримано коефіцієнти до рівняння роз­ повсюдження макротріщини повзучості Ключові слова: повзучість, пошкоджува­ ність, розповсюдження тріщини повзучості, скінченноелементна розрахункова модель


introduction
Ensuring the reliability and durability of elements of modern aerospace equipment, steam and gas turbines, etc. requires the creation of effective methods of computational analysis.Thanks to the application of such methods at the design stage, it becomes possible to significantly reduce the costs of experimental proof of products.The operating conditions of such structural elements are often characterized by elevated temperatures and loads, leading to creep deformation and associated damage of the material.
Long-term operation in creep conditions, the accumulation of damages in the material lead to the appearance of cracks in the structural elements.At this stage, it is necessary to take a decision or to urgently replace such an element, or to extend its working time.In connection with this, an actual and practically important task is the development of methods and algorithms for analyzing the total fracture time and the shape of the destroyed area.The latter is necessary to prevent the violation of the product functionality at the time, is less time for its complete fracture.Such methods and algorithms should be used both at the design stage of new designs and when estimating the residual resource of existing ones.

literature review and problem statement
The issue of the crack growth under the creep conditions of the material began to attract the attention of researchers from the middle of the last century, when the approaches to Fracture Mechanics began to be applied [1].After the developed method for solving the problem of a crack growth under elastic-plastic deformation, based on the use of the J-integral [2], similar models were created in creep theory in [3][4][5].In many cases, the Finite Element Method (FEM) [6] was used to determine the components of the stress-strain state in the structural element or its model under creep conditions [7,8].Later, for the first time in the work of J. Hayhurst's school [7], approaches of Continuum Damage Mechanics were used to estimate fracture under creep [9].A review of the work carried out before the end of the 1980s in the direction of the one under consideration is contained in the survey [10].
In connection with the importance of the question of a complete analysis of the long-term strength of the elements of modern industrial structures, research in the field of fracture under creep conditions continues in recent years.As at the initial stage, work is under way to evaluate the two stages of fracture -hidden and open, in the crack growth process.
In work [11], the effect of preliminary damage on a crack growth in a weld joint was analyzed.The calculations involve the so-called «physical based» state equation with several damage parameters.The place of crack initiation was estimated only according to the calculation of the damage distribution without estimating its further spread, which for the joint was considered and determined by complete fracture.
The application of the considered approaches to the analysis of fracture of standard specimens with a cut modeling the initial crack was carried out in [12].Further crack growth is calculated using two methods.The first leaves elements in the finite element model, in which the critical value of the damage parameter is reached.The second removes the nodes («node release method») from the model, and a special function is created for the implementation of the program in the ABAQUS software complex.It is assumed in advance that the crack growth follows a given line («postulated crack path»), which is justified by the experimental data.This approach is not universal for estimating the crack growth direction in structural elements of real geometry and the shape of the destroyed region.
The materials of the report [13] contain a method, which is also based on the damage analysis in a rod with a given direction of failure.The means of transition from the continuum damage model to the crack opening are discussed.In the paper, only one-dimensional problem is considered and there is no method of transferring the applied method to the case of a complex stress state.
The monograph [14] considers classical methods and models of both Fracture Mechanics and Continuum Damage Mechanics.A separate subsection is devoted to methods for estimating the behavior of a creep crack.The authors consider the theoretical achievements of Fracture Mechanics in the form of fracture motion equations, kinetic equations for scalar and tensor damage parameters, and the like.These approaches and methods can be applied to the analysis of the behavior of real structures, but the monograph does not specify an example of the corresponding solved problems.
The authors of [15] carried out experimental studies on the analysis of the crack occurrence in a standard specimen with a notch simulating an initial crack.With the help of the FEM simulation of the damage growth process, the fracture of specimens was carried out, which for the considered conditions and determines the fracture time.An analysis of the shape of the destroyed material in time is not carried out in this study.
Thus, at present, the problem of the fracture of structural elements operating under creep conditions in a complex stress state requires the decomposition of existing and the development of new approaches.Such approaches have links between the methods of Continuum Damage Mechanics and Fracture Mechanics.
The end of the process of hidden damage and the appearance of macrodefects in the form of a crack for many materials is not a moment of complete fracture of the component.The crack growth time that has arisen can be commensurate with the time of the previous development of hidden damage.
To more accurately estimate the lifetime of the structural element, which is used for the creep of its material, a complex analysis is necessary.Simulation should consist of two stages: the first is the growth of hidden damage, the second is the growth of the macrocrack, taking into ac-count both acquired damage and those that arise and grow at this stage.
The creation of such approach to simulation allows to develop a technology for a more accurate assessment of deformation and the long-term strength of structural elements.This makes it possible to design more reliable equipmentsteam and gas turbines, gas turbine engines, liquid jet engines and the like.

The aim and objectives of research
The aim of research is proposition of an effective algorithm for modeling the fracture process of a structural element, the material of which is in creep conditions and the accumulation of associated damage.
To achieve the stated aim, the following tasks are set: -to create auxiliary algorithms for sequential reconstruction of the FE model by removing elements in which hidden damage has ended, reformulating the boundary conditions in displacements and stresses, rearranging the stiffness matrix of the system to solve the problem by the FEM using; -to develop a methodology for the application of the developed numerical algorithms and the FEM CREEP software package to continue the calculation to the point at which the crack motion that occurs in this case is finished by the complete resolution of the structural element to independent parts; -to perform calculations necessary for fracture analysis using the example of a model of a notched specimen made from a high-temperature nickel-based alloy; -to provide a mathematical description of the numerically obtained process of fracture motion using the Fracture Mechanics equation.

mathematical formulation of the creep problem in the plane stress state
Let's give the mathematical formulation of the two-dimensional creep problem, which is accompanied by damage.Let's consider a region V with isotropic properties in the Cartesian coordinate system x i (i = 1, 2).Let's assume that traction p i (x) is applied on a part of the surface of the body S 1 as well as the displacements u u i S i | 1 = are known on the other part of the surface of the body S 2 .Let's suggest that the volume forces are absent.The body is in a homogeneous temperature field.By the Lagrangian description of the kinematics of an arbitrary point of the body, let's consider small strains and displacements.Let's denote by u the displacement vector with the components u i (x, t), σ, ε are stress and strain tensors with components σ σ ij ji x t = ( , ) and ε ε ij ji x t = ( , ), which are functions of the coordinates of the body point x i (i = 1, 2) and time t.Let's assume that at any time the strain tensor is the sum of elastic and irreversible creep strain tensors: where e el -the elastic strain tensor; c -the creep strain tensor with components c c x t ij ji = ( , ), c x ji ( , ) , 0 0 = i, j = 1, 2.. Let's note that at the initial time t = 0, when c ij (x, 0) = 0, an elastic deformation occurs.
The system of equations for determining the stress-strain state of a body under creep is as follows [18]: where, in addition to the previously defined notation, n -the unit vector of the outer normal to the surface S 2 of the body with the components n i , i = 1, 2.
To formulate creep state equations accompanied by damage, let's apply the classical model of the incremental theory, the effective stresses for describing the complex stress state are specified by the Bailey-Norton and Rabotnov-Kachanov models [19]: Here B, D, n, k, m are constants, which are determined experimentally with creep and long-term strength tests at a given test temperature [19]; s ij -components of the stress deviator; σ i -von Mises stress; σ eq -the equivalent stress, which is determined by the accepted fracture criterion.In this study, according to [16], σ eq = σ i accepted.
In system (3), the kinetic equation for the scalar damage parameter w = w(x, t) is applied.It is considered that at the beginning of the solution of the system of differential equations (1), (2), which corresponds to elastic deformation, the parameters w at all points of the region are equal to 0. This corresponds to the model of virgin material.Due to the creep growth and aftereffect in the case of a complex stress state, the material begins to be damaged, and after some time at some point the damage parameter reaches the so-called critical value w * .In the classical model of Yu.Rabotnov [19], it is equal to 1, which corresponds to the process of ending the hidden damage and the appearance of a macrocrack at this point.The solution of system (1), (2) should continue further with the changed geometry of the region.The new initial conditions must correspond to the values of the components of the stress-strain state acquired at that time and the damage parameter at all points in the problem definition area.
The specified approach should be applied until the fracture criterion defined for the given problem is satisfied.Let's take as the fracture criterion the attainment by the fracture region (macrocrack) of a certain critical length.

method for solving the creep-damage problem
Let's solve the initial-boundary value problem with the help of the Finite Element Method in combination with the third-order predictor-corrector method for time integration.In the program complex FEM CREEP [17], which is used for numerical simulation, the following vector-matrix formulation of the problem is used: After the procedure of ensembling the system of finite elements [20], the following system of differential equations is obtained: where [K] -system stiffness matrix; {u} -global vector of nodal displacements; {F } -vector of nodal loads from traction; {F с } -vector of nodal forces caused by creep strains.When time integration, at each step, the components of stress and strain vectors and damage parameters in finite elements as well as nodal displacements are determined numerically.
The FEM CREEP software is designed for a triangular three-node finite element.For the generation of finite element models, special preprocessor programs for automatic grid generation are used [17].

Algorithm for estimating the macrocrack growth from generation point
Let's describe the developed algorithm.Its block diagram is shown in Fig. 1.
The solution of the problem of creep and fracture of the structural element must pass before the introduction of the critical size of the macro-destruction region, which describe by the crack length l.Until the end of the development of hidden damage in the material, l = 0.After the damage para meter in any finite element has received a critical value w = w * , the current crack length l n is determined.It is believed that at this stage (n = 1) the solution of the problem is equal to the characteristic size of the destroyed finite element.Let's note that in numerical calculations, w * = 0.97..0.99 is adopted.
Then, for each finite element, the values of the components of the stress-strain state and the parameters of the damage properties acquired at that time are stored in the respective arrays.The «destroyed» element and the nodes directly belonging to it are removed from the finite element grid.In order not to miss the case in which a stress state is simulated around the resulting fragment of the structural element with insufficient accuracy, a grid with a sufficiently large number of elements is pre-selected in advance.The structural stiffness matrix of the structure and the vector of the right parts are rearranged, if necessary, after an appropriate check of the presence of remote nodes in them, the boundary conditions are replaced.A new system stiffness matrix is formed, and the creep calculation process is started by solving a system of differential equations of smaller dimension, which is determined by the number of remote nodes.The process is up to the moment of fracture of the next element or group, the value of the damage parameter in which reaches the range w * = 0.97..0.99.For this n-th stage of life of a structural element, the size of the destroyed element or elements is determined by the new current value of crack length l n .The calculation algorithm is repeated until the current crack length has reached its critical value l * .

numerical simulation of the creep crack growth process in notched specimens
As an example of the described numerical simulation technique, consider the creep and fracture of tensile plane specimens with triangular, symmetrical with reference to longitudinal axe, notches with an angle of 45° (Fig. 2).The tensile load is 10 MPa.Samples are made of heat-resistant nickelbased alloy EI-867 (XH62MBKЮ).The test temperature is 950 °C.To obtain the values of the constants entering into the governing equations (3), experimental data obtained in [20] are used.By processing the experimental curves given in this paper, the following values are obtained: n = m = 2.4, B = 9.87⋅10 -10 MPa -n /h, D = 5.45⋅10 -8 MPa -m /h, l = 3.4.
It is established that the hidden damage accumulation is completed in the element that is located at the top of the notch, at t = 7.7 h.Further calculations continue in time for a number of rebuilt models with the removal of the «de-stroyed» elements.In Fig. 3 as an example, the form of the finite element grid for the time t = 12.3 h is given.Due to the symmetry of the sample, the calculation model is built for its fourth part.After numerous experiments in the study of the convergence of solutions to calculations, a grid with 6,884 elements was adopted.
Destruction of the specimen by removing the material of one or several elements occurs before the time t = 12.35 h (Fig. 4), after which it passes into the so-called «avalanche» stage.By the time t = 12.36 h in a large number of elements located diagonally from the left bottom to the right upper part of the sample, the value of the damage parameter reaches a critical value (Fig. 5).The length of the crack becomes equal to the length of this diagonal, taking into account the symmetry of the model means complete fracture of the sample by dividing it into two parts.A numerical analysis is performed that makes it possible to fix the average crack length for each n-th stage of the fracture simulation algorithm.In this study, the length of the crack is determined by counting the distance from the top of the notch to the center of the «destroyed» element as far removed from it.As can be seen from Fig. 5, the fracture pattern in this case can be approximately described by a one-dimensional model of crack motion, and it is possible to directly estimate the geometry of the destroyed zone.Let's note that the numerically established direction of the crack propagation corresponds to the experimental data for specimens with notches [7].
The calculations show the comparability of the latent (7.7 h) and open fracture (4.66 h).For similar time relationships, such analysis is necessary to establish the part of the time with which the structural element can still exist with macrocracks.It should be noted that, of course, in practice there are considerably more durations of work for creep.
To describe the motion of macrocracks, differential equations that establish the dependence of length on such quantities as the value of the contour integral, the Norton law exponent, the dimensions of the specimen etc, are used.The approach of this paper proposes to obtain an equation describing the crack growth, according to the numerical simulation of the creep of a damaged element.An example of such calculation for a notched specimen is given above.In this case the constants describing the creep and damage of the given material -B, D, n, k, m, are indirectly included in the created equation, because only due to their definite values the numerical data are obtained, by which processing the coefficients of the equations are determined.
For the considered model of the notched specimen, calculations with other axial tension values are additionally carried out, and similar results of time dependence of geometric changes in the fracture region are obtained.Differences occurred only in the growth speed of the destroyed zones, hence in the growth speed of the macrocrack.
The results of numerical modeling of fracture of specimens with notches for three cases, tensile stresses p (12, 10 and 9 MPa) are shown in Fig. 6 by a plot of the length of the crack versus time (curves 1, 2, 3, respectively).As can be seen from Fig. 6, each graph can be represented by two partsthe first with a constant speed (different slope for different stresses of practically straight segments) and the stage of avalanche fracture.
This result agrees well with the results known in the literature.In the monograph [21] it is shown that the motion of shallow cracks is described precisely by such graphs -segments of lines of constant inclination.For this case, it is pro-posed to construct a differential equation as the dependence of the crack propagation rate on the tensile stress or the stress intensity factor and the constants B, n, the half-width of the sample, and others.The following approach is also applied.

Fig. 6. Dependence of the maximum crack length on time
The numerical results obtained for the first, linear portion of the fracture growth graphs, are processed using an equation of this type: In equation ( 6), А, s, α -constants, which are determined by solving a system of three equations for three pairs of values (l i , t i ) obtained from three curves 1, 2, 3 in Fig. 6.The values of the constants are found by processing: A = 11.562mm/h, s = 2.846, α = 5.17⋅10 -4 MPa -1 .Fig. 7 contains a plot of the length of a crack versus time for the basic time of macro fracture, constructed from the obtained values of the constants.At the same time, the deviations of the data obtained by the finite element calculation, and by the relation (6), do not exceed 0.1 %.For the obtained total fracture times, the difference between the values determined with and without avalanche regions is about 1 %.In this regard, the total fracture time of a notched specimen with a completely satisfactory degree of accuracy is determined by the sum of the hidden damage accumulation time and determined by equation (6).

discussion of the numerical simulation results
Thus, the created approach is in a certain sense an alternative proposed earlier.The widely used fracture motion equations for creep [3][4][5][6][7][8][9][10][11][12][13][14][15]21] include as constants the characteristics of the sample material and its dimensions, others experimentally determined values such as the C* -integral, the rate of displacements or strains, and the like.All these important parameters are indirectly included in the resulting fracture motion equation (6) and are determined from the data of numerical analysis.The material properties for the motion of the crack front during creep are taken into account due to its constants, which enter the state equations (3), and the constantly increasing values of the damage field ahead of the crack front.Geometry is accurately taken into account at the stage of solving the boundary value problem.
The advantage of the proposed approach is the possibility of a sufficiently accurate description of the fracture pattern by analyzing the computational model with the removed finite elements.This advantage is due to the properties of the developed numerical simulation algorithm, built on the possibility of obtaining the initial-boundary creep problem solution for any moment and fracture stages -hidden or open.The exact shape of the destroyed area is determined by the cumulative form of the deleted elements.The prospect of the proposed solution can be considered obtaining the constants of the fracture motion equations from numerical modeling data for the most common types of structural elements.Such equations, after appropriate experimental checks, can be used in the design without using a calculation based on the FEM.
Also an important conclusion after the performed numerical simulation will be given from the ratio of the time of hidden damage and the crack motion.In the case when the latter turns out to be small, it is possible to recommend for performing calculations for the analysis of only hidden damage.
The creep problem of a sharp-cut specimen is considered in this study with the aim of presenting a method that can be applied to any structural elements of complex geometry.It is clear that the proposed approach at the current stage of development has certain limitations.It requires special control of the quality of the finite element grid.The number of elements in the path of the crack motion, which is assumed before the simulation start, should be one that can adequately describe the stress concentration in the fracture region.It is possible that in some cases it will be necessary to modernize the algorithm and apply the adjustment of the grid with a decrease in the size of the element.
Also in similar studies, it is necessary to accurately solve the problems of adequate description of stress concentration.When simulating such kind of concentrators, it is necessary to involve a comparison with the corresponding experiment and additional FE procedures.
In many cases, when the element is destructed, the displacement of the points during creep is small and can be described by geometric linear relationships of the Solid Mechanics, large deformations can occur due to crack opening.Thus, for such cases, the formulation of the problem is considered has to be extended by transferring the obtained values of the components of the stress-strain state and the parameter of the model's damage to the model, constructed on the general formulation of the problem with allowance for finite strains.

conclusions
1.The description of the developed numerical algorithms and the general algorithm for numerical modeling of crack propagation in the case of creep taking into account the material damage is given, the limits of the admissibility of its application are discussed.
2. The suitability of the approach is established, according to which the total fracture time of the notched specimens chosen for its demonstration is defined as the sum of the time until the end of the hidden damage in a certain finite element and the time that is necessary to destroy the material in the direction of motion of macrodefect.The use of the Finite Element Method together with the numerical integration of the system of differential equations allows to track the complete picture of the fracture in the plane and in time by eliminating the «destroyed» elements.The determination of the dimensions of these elements makes it possible to numerically describe the process of macrocrack growth by applying the differential equation of its motion over time.
3. On the example of the creep problem of a notched specimen from a heat-resistant nickel-based alloy, it is demonstrated that the proposed algorithm for modeling the motion of a creep cracks allows its qualitative and quantitative description to be performed.This is ensured by using the algorithm for extracting «destroyed» finite elements.Also, thanks to this, it is possible to describe the shape of the destroyed area.Due to the solution of the initial-boundary creep-damage problem, the values of the extraction time of elements are set that coincides with the propagation time of the crack to the current length.
4. On the basis of numerical simulation, coefficients to the crack growth equation in time are determined.For the section of rectilinear growth, the widespread form of the dependence of the of crack length rate in dependence on the tensile stress as a function of the hyperbolic sine and power law is applied.The possibility of applying the proposed procedure as an alternative to the existing ones, which are usually performed by valuable and long-term experimental tests, is discussed.
where γ and b indicate that the value refers to the node or element of the finite element grid; [B] -strain matrix; [C]matrix containing elastic characteristics of the material (matrix of elasticity); {u γ } -vector of nodal displacements; {σ b }, {ε b }, {с b } -vectors of stresses and elastic and creep strains in the element; [N 2d ] -matrix, which expresses by means of shape functions the value of any displacement of the element through the node; р -traction.The integral in parentheses on the left-hand side of the first equation (4) corresponds to the stiffness matrix of the element [K b ].

Fig. 7 .
Fig. 7. Dependence of crack length on time for the basic time of macrofracture