Multi-parameter approximation of stress field in a cracked specimen using purpose-built Java applications

ABSTRACT. A quality of multi-parameter approximation of the stress and displacement fields around a crack tip in a non-brittle material test specimen is studied in the paper. The stress field approximation using Williams power series is intended to be utilized for estimation of the nonlinear zone extent which potentially plays a role within methods for determination of true values of fracture parameters of materials exhibiting nonlinear failure. Considering the fact that in the case of elastic-plastic and especially quasi-brittle materials the size of this zone is substantial in comparison to the specimen dimensions, it is necessary to take a large region around the crack tip into account for this task. An automatic utility created to determine the values of coefficients of the higher order terms of Williams power series by usage of over-deterministic method applied on results of finite element analysis of arbitrary mode I test geometry is one of the tested procedures. The second one provides a backward reconstruction of the crack-tip stress field analytically by means of truncated Williams expansion. The developed procedures are based on the support of Java programming language and considerably simplify analyses of the mechanical fields’ description in a farther distance from the crack-tip. The presented research is focused on optimization of selection of FE nodal results for improvement of accuracy of the approximation.


INTRODUCTION
he concept of fracture process zone (FPZ) has been introduced into analysing fracture behaviour of nonlinearly failing (especially softening) materials already several decades ago [1][2][3].Estimation of the size, shape and other relevant properties of the FPZ evolving at the tip of propagating crack in these materials is not commonly utilized within fracture analyses of these materials; however, the situation has changed recently.The FPZ issue is also among the main goals of the presented research.This is motivated by the fact that fracture properties of softening materials exhibit strong size-dependence, know under the term size-effect, linked mainly to the mutual proportions of the internal length of their coarse material structure, the extent of their large FPZ which varies during the crack propagation, and the structural dimensions.Thus, properties of the FPZ are prospected to be utilized, T among others, within methods for evaluation of fracture parameters and subsequently the mechanical characterization of these materials.Due to this effect (connected also with geometry and boundary effect), investigation the role of FPZ in the fracture of quasi-brittle materials, particularly concrete, has attracted lot of interest from engineering community in the field of concrete fracture mechanics.For the FPZ extent estimation, a precise description of the stress state in the cracked body is necessary.For that purpose, multi-parameter fracture mechanics is employed as the proportion of sizes of the FPZ and the whole specimen is much larger in the case of the quasi-brittle materials than e.g.plastic zone in metals [4], and therefore the stress field farther from the crack trip must be described precisely.Main focus of the contribution is on testing of two own developed automatic utilities; the first enables computation of values of coefficients of the higher order terms of Williams power series using the over-deterministic method (ODM) applied to results of FEM analysis, the second provides analytical approximation of the stress field in the cracked body via truncated Williams power expansion taking into account the values of coefficients of the higher order terms calculated using the first programmed application.From this point of view the latter utility provides a kind of backward reconstruction of the stress field analysed by FEM and processed via ODM.The influence of the number of terms considered for the power series reconstruction of stress field as well as the way of selection of the FEM nodes for the calculation of the values of coefficients of those terms are investigated in the paper.Verification of the overall analysis results is conducted on an example of the wedge splitting test specimen as it has become very frequently used in the field of testing of silicate-based composites (let us name e.g.[5][6][7] from recent period).

MULTI-PARAMETER FRACTURE MECHANICS CONCEPT
Near-crack tip fields for mode I fracture problem he stress and displacement fields in a two-dimensional homogeneous elastic isotropic body with a crack loaded by an arbitrary remote loading can be expressed in a form of an infinite power series, the Williams expansion [8].For the stress tensor  and displacement vector u it holds: and , , , , respectively, where r and θ are polar coordinates (with their beginning placed at the crack tip and the crack propagation direction corresponds with the x-axis), f ij, and f i,u are known functions, E and  represent Young's modulus and Poisson's ratio.Values of coefficients A n depend on the cracked specimen geometry (crack length, specimen shape and boundary conditions) and they are usually determined numerically at present.Typically, coefficients An are expressed as functions of the relative crack length  and transformed (normalized with respect to loading and appropriate power of the specimen's dimension) into the form of dimensionless functions [9,10] as follows and  ( ) where  nom is the nominal stress in the central plane of the specimen due to the applied load, W and B is the specimen's width and breadth, respectively.The nominal stress is usually considered as being caused only by the splitting component P sp of the load imposed to the specimen via the wedge (see Fig. 1a, i.e.  nom = P sp /BW ), although it has been shown by the authors [11] that also the vertical component P v plays a significant role in many cases.Note that the approximations of the individual stress tensor components presented further in this paper are expressed using truncated forms of the Williams series.

Over-deterministic method
Computation of the values of the A n coefficients of terms of the Williams expansion is conducted using the technique solving the system of equations resulting from the displacement field expression in Eq. ( 2) referred to as the Over-T Deterministic Method (ODM, [12][13][14]).Basically, the ODM is a linear least-squares formulated regression, it provides the solution of system of 2k equations, where k means the number of selected nodes around the crack tip, for first N selected terms of the power series.Components of the displacement vector u a v in the selected nodes of the finite-element mesh and also their coordinates serve as inputs to this technique and values of the coefficients of the terms of the series A n , where 1, n N  present its outputs.In the case of this study, the number of the selected nodes k was set to 49 (due to keeping consistence with previous studies in [15,16]) and maximal order of the Williams series term was N = 12, so the system of equations was strongly over-determined (the condition of the over-determination is 2k > N).For the FE calculations the ANSYS computational software [17] was employed.

ODeMApp
he ODM procedure has been implemented in Java programming language [18] as a multi-platform and objectoriented computational tool, referred to as ODeMApp (Over Deterministic Method Application).Design of each class of the software application came out from the requirement of transparency and further extensibility.It is a command line operated tool; the user specifies the path to external input files containing necessary FE parameters (the coordinates of selected nodes and their displacements) and additional input information (parameters of material and geometrical description of the test specimen) for performing the solution using the ODM procedure.The program then allows defining the type of the solved 2D problem (plane strain or plane stress).In the next step the user enters the highest index of the Williams series term up to which the coefficients shall be calculated.Then the ODM procedure takes place.The calculated values of A n coefficients are eventually converted to their dimensionless form g n within this procedure.

ReFraPro
ReFraPro is a Java application programmed for the advanced determination of the fracture characteristics of silicate-based materials failing in a quasi-brittle manner [19].The tool reconstructs the progress of quasi-brittle fracture from the measured load-displacement curve and the knowledge of basic mechanical properties of the material (ReFraPro si an acronym for Reconstruction of Fracture Process).The application implements a developed technique for estimation of the size and shape of the FPZ [20][21][22].The technique is based on an amalgamation of several modelling concepts dealing with the failure of structural materials, i.e. multi-parameter linear elastic fracture mechanics, classical non-linear fracture models for concrete (equivalent elastic crack and cohesive crack models), and the plasticity approach.The application utilizes the description of the stress field in a cracked body considering the higher order terms of the Williams power expansion so that the field in a larger distance from the crack tip is well approximated.This description allows creating graphical representations of stress fields (presented below) for particular set of terms of William expansion which are used for this study of the dependency.This paper presents only the part of the application that deals with the analytical reconstruction of the stress field by means of the truncated Williams series (similarly as in [23,24]), other capabilities are presented in e.g.[19,22].

STUDY OF THE STRESS FIELD IN A BODY WITH CRACK -WEDGE SPLITTING TEST SPECIMEN
he wedge-splitting test (WST, [25][26][27]) specimen in one of its variant is used for the study -see Fig. 1a.A cubeshaped specimen with the edge length of W = 100 mm is supported in two points and loaded through steel platens.Specimen dimensions provided in Fig. 1a

Numerical model
The created FE model of the studied test, namely its symmetrical half, with typical applied boundary conditions (the vertical and horizontal components of the applied load; fixed nodes in both directions at the support and perpendicular to the axis of symmetry along the ligament) is illustrated in Fig. 1b.The test specimen is considered to be made of quasi-T T brittle cementitious composite material, e.g.concrete, with modulus of elasticity E = 35 GPa and Poisson's ratio  = 0.2 and the loading platens of elastic isotropic material with E = 210 GPa,  = 0.3.The FE mesh of the task is formed from almost 70,000 FE nodes in order to enable fine enough nodal selection for the ODM application in all directions from the crack tip.Details of the FE mesh around the crack tip and the support is shown in Fig. 1c and d.

Description of parametric study
The parametric study is constructed in a way that shall reflect the two main issues of this research: i) which number of terms of the Williams power series is suitable to consider for reconstruction of stress field in relevant fracture analyses and ii) how to determine the values of coefficients of those terms with sufficient accuracy.In both these issues the distance from the crack tip plays a significant role.The ODeMApp application is created for arbitrary number of nodes picked from the FE model; however, k was set to 49 in this study (explained above, adopted from [15,16]).The space of the parametric study presented in the paper consists of several directions/axes.The basic two are formed by the radial and the tangential (using angular sections) way of FE nodes selection for the ODM technique.The nodal selection concerning the angle range can be optional from 0° to 180° (in the case of the used symmetrical half of the specimen, i.e. 360° for the whole model).In this study, four sub-variants were considered as is shown on Fig. 2a to d (representing models with different values of relative crack length ).Particularly, the Fig. 2a shows the case marked as 90°(10°) which describes sector of angular measure equal to 90° with initial side of the sector being rotated from the crack propagation direction by 10°.Analogically, cases marked as 90°(45°), 90°(80°), and 180°(10°) are depicted in Fig. 2b, c and  d, respectively.Note that especially in the cases 90°(10°) and 90°(80°) the sectors are oriented to directions, where boundary conditions are imposed, i.e. in the case of 90°(10°) towards the groove with loading platens and the 90°(80°) case to support on the bottom surface.For the distance distribution of the nodal selection, three options are considered in this study: the constant (f(x) = const., marked as con), the quadratic (f(x) = x 2 , marked as qua) and the exponential (f(x) = e x , marked as exp) function.For clarity, the distributions are sketched in Fig. 2 as well to complete the graphical description of the conducted study.Particular expression of the functions depends on the chosen number of nodes within the length interval between the crack tip and the boundary for individual directions.In this case four nodes, i.e. nodes in four layers from the crack tip, were chosen.Thus, altogether 12 cases of the FE nodes selection are taken into account in the study: four cases of tangential selection × three cases of radial distribution.Selection of nodes whose state variables' values (displacements) are used within the ODM technique for calculation of the values of coefficients of the higher order terms of the Williams series is performed from a regular FE mesh (see Fig. 1c  and d).The procedure of selection of the nodes according to the two above-described criteria was programmed within the post-processing of results in the used ANSYS computational software.Coordinates and displacements of the selected nodes for each considered nodal selection case (for a model of each value of relative crack length) are written in a text file after the run of FE analysis.These created data files are then loaded into ODeMApp and processed.

Geometry functions
ig. 4 shows the results from the ODeMApp procedure, the dimensionless values of Williams terms coefficients gn (for 1;12 n  ) as functions of the relative crack length  determined for all considered variants of nodal selection.It can be seen that trends of the g-functions for individual nodal selection variants are not uniform, especially in the cases of short and long cracks (the short cracks are represented by  < 0.2 and the long ones by  > 0.8).In the case of short cracks the 90°(10°) variants significantly differ from the rest.On the other hand, the 90°(80°) variants exhibit different trends in the long cracks region, especially for n > 6 (see Fig. 4g to l).The nodal distance distribution function also plays a role; there can be seen a substantial difference between the exponential distribution in comparison with the other two variants.The most suitable distribution function of nodal selection is the constant one -it seems to be best for the majority of angle selections.The most stable angle selection variant is 180° (0°).These statements are in consistence with Fig. 5 where the g-functions with stable progress are kept while the fluctuating variants are deleted.The selection of variants that are regarded as incorrect was based only on visual judgement, no rigorous criterion was applied.However, according to opinion of the authors, such an approach is sufficient for the presented study.
It can be stated that: the way of angular selection of nodes has greater effect than the distance distribution function for the progress of g-functions captured by the following graphs in Fig. 4 and 5; the progress of g-function is strongly influenced by the used type of distance distribution and the place of interest (determined by the initial angle, the selection sector and the relative crack length).In this aspect the results are in accordance with the conclusions in [16]; this is the reason for choosing three different values of  for stress the field analysis using the ReFraPro application.Note that the FE model, before it was used for evaluation of the stress field via ODM, was successfully verified for the value of the first term of William expansion in previous works (e.g.[28]) via its comparison to results from literature [9,10,29,30].

Stress field reconstruction by means of truncated Williams series
For presentation of the stress fields approximation using the ReFraPro application, three values of the relative crack length were selected to cover the whole interval, from the short cracks with  = 0.25, the middle ones with  = 0.5, and the long ones with  = 0.75.Results of this rather thorough numerical study on the stress fields' approximation are displayed in Figs. 6 to 11. Contour plots of the  1 principal stress and the normal stress components (in this case marked as  yy  because of rotated coordinate system) are displayed; they were taken into account as they might serve as the equivalent stress thresholds in the simplest failure condition for estimation of the nonlinear zone (plastic zone) extent.Nomenclature in tables in these figures is as follows: Main columns represent the number of the higher-order terms (variants with N = 1, 2, 4, 7, 11 were selected) considered for the backward reconstruction of stress field by ReFraPro application.Each of these columns has three sub-columns for the considered nodal distance distribution functions (con, qua, exp); rows represent the angle selection variants (90°(10°), 90°(45°), 90°(80°), 180°(0°)).In the left column, a corresponding finite-element solution using ANSYS software (considered as the exact one, or in other words that one which takes into account theoretically infinite number of William series terms) is given for comparison.It is obvious from Fig. 6 that variants of nodal selection have almost no effect on the stress field reconstruction if only first or first two terms were taken into account.Therefore, contour plots of the selected stress tensor components calculated for N = 4, 7, 11 are in shown Figs.7 to 11.The contour plots provided by both tools, the ANSYS FE software and the ReFraPro application, are created for the same load (Psp = 2.9 kN considering the wedge slope angle equal to 15°).The stress scale is uniform for all displayed contour plots: grey colour for values < 0 MPa, blue for interval (0;0.5)MPa, cyan for (0.5;1.25) MPa, green for (1.25;2.5)MPa, yellow for (2.5;5.0)MPa and red for values > 5 MPa.Discussion of the results of the stress field reconstruction is based again only on the visual comparison of the contour plots in this paper.It is intended to apply a suitable method for the inaccuracy quantification in further research.After detailed observations of the results, it can be stated that: -Determination of a particular number of terms that allows a correct enough approximation of the stress field is tricky in a real case since it depends also on the distance from the crack tip at which the field is characterized.-Using a low number of parameters, for example only the first two parameters, the stress intensity factor (corresponding to g 1 ) and the T-stress (g 2 ), the area of reasonably accurate approximation is very small, see the red colour contour in the vicinity of crack tip for N = 1 or 2 (in Fig. 6).-If the knowledge of stress field from a larger distance of the crack tip is necessary, one must take into account several terms of the Williams power series.-However, it's not possible to determine a certain number, because it depends, particularly in the studied WST specimen, on the relative crack length, the investigated component of the stress tensor and the distance from the crack tip.It can be expected that it is dependent also on the shape and boundary conditions of the cracked body.But generally for this study, the number of terms should be greater than 4. Considering the facts from previous paragraph, it is evident that the attention should be paid to the choice of the nodal selection.From the g-function graphs can be seen that the usage of nodal selection from the whole test specimen with constant distance distribution provides the most stable and accurate results.

CONCLUSIONS
n the present work, issues regarding the stress field approximation around the crack tip in a 2D body loaded in mode I are discussed.Two own developed Java applications are utilized for the task; one (ODeMApp) for determination of the values of coefficients of the Williams power series terms from results of FEM analysis using a square-root regression based technique, the second (ReFraPro) for reconstruction of the crack-tip fields via the Williams series considering the values of the terms coefficients determined using the former one.These developed software tools are employed with convenience in a study concerning the way of selection of the nodal results from the FEM solution on the quality of the stress field capturing (using ODeMApp) and reconstruction (using ReFraPro).Several variants of way of for selecting of the nodes were compared and discussed with regard to accuracy of the stress field approximation.Significant advantages of the usage of multi-parameter description in the analysis of the stress field are reported; large errors of the approximation caused by a limited number of terms are demonstrated.Results of this study have significant implications to intentions on estimations of the FPZ extents in quasi-brittle materials.The technique of the FPZ estimation implements the crucial parts of the above-mentioned procedure; functioning of the technique is currently being verified via numerical simulations of the fracture process in the test specimens by means of a lattice-particle (spring network) model and is planned to be validated with the help of acoustic emission technique and/or dynamical X-ray radiography in future research.

ACKNOWLEDGEMENTS
are of following values: f = 30 mm, dn = 20 mm, h = 15 mm, Wef = 85 mm, i = 25 mm.The crack length a is considered as the vertical distance from the crack tip to the loading bearings's axles, i.e. a = c + dn -h.Variable relative crack length 0.025 is used to simulate stages of the progressive fracture through the specimen ligament.

Figure 1 :
Figure 1: a) Used WST geometry, b) FE model (symmetrical half of the specimen with loading platen), c) and d) details of FE mesh.

Figure 2 :Figure 3 :
Figure 2: Variants of selection of nodes within an angle range and distance from the crack tip.

Figure 4 :
Figure 4: Dimensionless values of g n as functions of the relative crack length determined by ODM; plots of all considered variants of the node selection are displayed.

Figure 5 :
Figure 5: Dimensionless values of g n as functions of the relative crack length determined by ODM; only plots of variants of the node selection without a significant error are displayed.

Figure 6 :
Figure 6: Contour plots of the approximation of  1 principal stress in the test specimen with the relative crack length  = 0.25 for all the considered variants of the nodal selection and number of Williams series terms N = 1, 2, 4, 7, 11 compared to the FEM solution.
inancial support from the Czech Science Foundation (projects No. 13-09518S and No. 15-07210S) is gratefully acknowledged.This paper has been worked out under the "National Sustainability Programme I" project "AdMaS UP -Advanced Materials, Structures and Technologies" (No. LO1408) supported by the Ministry of Education, Youth and Sports of the Czech Republic.

Figure 7 :FEMFigure 8 :FEMFigure 9 :FEMFigure 10 :FEMFigure 11 :
Figure 7: Contour plots of the approximation of  1 principal stress in the test specimen with the relative crack length  = 0.5 for all the considered variants of the nodal selection and number of Williams series terms N = 4, 7, 11 compared to the FEM solution.