A two-dimensional formulation of the 3-PG Model of implicit gradient elasticity has been developed in Part I. The predicted near-tip fields for Mode-I and Mode-II crack problems have been derived in Part II. It has been found that both the classical Cauchy stress and the nonclassical double stress are singular with the order of singularity r−12. In the present chapter, the two-dimensional model formulation is implemented in a finite element code. For verification of the resulting finite element model, a square section with a circular hole subjected to displacement-controlled tension loading is considered and discussed. The main concerns of the chapter are, on the one hand, to validate the analytical solutions of Part II. On the other hand, the chapter aims to investigate the effect of nonclassical material parameters on the stress intensity factors.
- implicit gradient elasticity
- finite elements
- square section with a hole
- mode-I and mode-II crack problems
- stress intensity factors
- angular functions
Mode-I and Mode-II crack problems have been discussed analytically in Part for the 3-PG Model of implicit gradient elasticity. Solutions for the near-tip fields have been obtained by employing the method of asymptotic expansion of Williams’ type (see Williams ). It has been proved that both the classical Cauchy stress and the nonclassical double stress are singular with the order of singularity . Even more, the first two terms in the asymptotic expansion of the Cauchy stress are identical to those in the context of classical elasticity. The leading terms of the asymptotic expansions of the classical and the nonclassical stresses are represented by the so-called stress intensity factors.
The present chapter deals with numerical simulations that employ the 3-PG Model to solve crack problems of Mode-I and Mode-II types. A finite element model for plane strain is developed in the framework of a weak formulation based on the principle of virtual work. To verify the finite element formulation and implementation, a representative example is considered: a square section with a circular hole subjected to tension loading. The predicted stress concentration factors are compared with the corresponding stress concentration factors predicted by classical elasticity. The main objectives are to confirm the assumptions and the analytical results of Part as well as to assess the effect of the nonclassical material parameters of the 3-PG Model on the stress intensity factors for Mode-I and Mode-II crack problems. It is perhaps of interest to remark that similar investigations for the case of micropolar elastic continua are provided in Diegele et al. .
The scope of the chapter is organized as follows: Section 2 gives some details about the implementation of the 3-PG Model in a finite element code. A square section with a circular hole subjected to tension loading is discussed in Section 3. The numerical simulations verify the finite element model and its ability to predict length scale effects. Further, they provide a first comparison to classical elasticity by calculating the corresponding stress intensity factors. Section 4 is devoted to an analysis of edge-cracked specimens. The analysis comprises, among others, the effect of material parameters on the stress intensity factors. Moreover, it indicates a very good agreement between the numerical and the analytical solutions of the angular functions. This confirms, a posteriori, the assumed symmetry conditions of the micro-deformation. The chapter closes with some concluding remarks in Section 5.
Throughout the chapter, the same notation as introduced in Part applies.
2. Finite element formulation
The following formulations refer to three dimensions and apply especially to two-dimensional cases when conditions for plane strain are imposed. Let us consider once more the equilibrium equations (see Section “The 3-PG Model as particular case of micro-strain elasticity” in Part I) to be solved,
and specify the corresponding boundary conditions (see Section “Boundary conditions” of Part I) as follows:
These equations reflect Dirichlet boundary conditions for the macro-displacement and the micro-deformation and Neumann boundary conditions for the Cauchy stress and the couple stress.
The first step toward a finite element formulation is to elaborate a weak form of the above boundary value problem.
2.1 Weak form of the boundary value problem
Let the fields and as well as the so-called virtual fields and belong to the following function spaces:
As usually, in favor of a short notation, we use the integral over to indicate the summation of single integrals over , which generally do not coincide. The meaning of the integration over is analogous. Eq. (13) is the weak form of the boundary value problem and is the starting point of the finite element formulation.
According to the finite element method (see, e.g., Hughes ), the domain is approximated by finite elements , so that
The elements are connected to each other at selected points called nodal points, or simply nodes, and the following notation holds.
: Set of global node numbers with macro-displacement degree of freedom.
: Set of global node numbers with micro-deformation degree of freedom.
The exact solutions and are approximated by
where and are the unknown values of and at node . The so-called shape functions and belong to finite-dimensional function spaces and , which approximate the function spaces and , respectively. Similarly, and are approximated by
with and being constants. The functions and are elements of finite-dimensional function spaces and , which approximate the function spaces and , respectively.
with the meaning of and being obvious.
We next employ the elasticity laws (see Section “The 3-PG Model as particular case of micro-strain elasticity” of Part I)
in order to replace the stresses and as well as the double stresses in Eq. (19):
The form of this equation suggests to define the element stiffness matrices:
and the element force vectors:
and to recast Eq. (25) as
As and may be chosen arbitrary, we obtain the system of equations:
The proceeding approximation of the weak form, constrained to plane strain states, has been implemented in the finite element code FEAP. Isoparametric elements are employed, that is, the space coordinates are represented by using the shape functions, ,
where are constants.
3. Square section with a circular hole
The main objective of this section is to validate the implemented finite element code. To this end, we consider the plane strain problem shown in Figure 1a, where the square section with a circular hole, located at the center of the section, is stretched in the -direction. The length of the section is , while the radius of the hole is . With respect to the Cartesian coordinate system , the boundaries are assumed to be free of classical and nonclassical tractions. At the boundary , the macro-displacement component , the component of the classical traction, and the nonclassical traction components are assumed to vanish. At the boundary , the macro-displacement in the -direction is given by , while and are imposed. The whole circular hole is assumed to be free of classical and nonclassical tractions. For a small circular hole, a nearly uniform stress component
will be required to realize the given boundary conditions.
The most simple case in classical elasticity, analogous to the boundary value problem above, is to consider the square section in the context of a plane stress problem subjected to the traction boundary condition . Attention is focused on the distribution of along the section -(see Figure 1a), as a function of the local coordinate
with . A so-called stress concentration factor is defined by
and turns out to be (see, e.g., Gould , p. 124).
Now consider the square section in Figure 1a in the context of a plane strain problem with the boundary conditions stated in the first paragraph within classical elasticity. The stress distribution along the section -has been determined by employing the standard elastic element of FEAP using the classical material parameters
As depicted schematically in Figure 1b, a radial mesh of stripes of quadratic elements with reduced integration, that is, elements with nodes, is used. Every stripe consists of elements and the elements are chosen to decrease in size the closer to the circular hole the elements are placed. In the context of classical elasticity, a mesh consisting of such stripes is used, whereas the mesh for the simulations of the 3-PG Model only needs such stripes of elements.
Figure 2 illustrates the distribution of the dimensionless stress along the section -. The value of at represents the stress concentration factor and turns out to be for classical elasticity.
We consider next corresponding distributions predicted by the 3-PG Model. To this end, we employ the finite element model developed in the last section. The classical material parameters are given by Eq. (39). For the purposes of the present chapter, we find it convenient to use the nonclassical parameters and
Note that since (see Section “The 3-PG Model as particular case of micro-strain elasticity” of Part I), the constrain applies. It becomes apparent, from the elasticity laws (20)–(22), that for , the 3-PG Model will approach to classical elasticity. Hence, we expect the distributions of the dimensionless stress along the section -to be close to the classical one whenever is sufficiently small. Indeed, Figure 2 confirms this expectation: It can be seen, that for sufficiently large values of , the two graphs almost coincide. Another way to illustrate this issue is to consider the effect of the nonclassical material parameter on the values of the stress concentration factor . We expect that for and , the values of will approach to the classical value . This is exactly what we can observe in Figure 3. Alternatively, we can consider the effect of the material parameter on the values of the stress concentration factor . For the assumed geometry of the specimen, the imposed boundary conditions and for a fixed value , the effect of on is illustrated in Figure 4 by the graph referred to as specimen . A convenient way to illustrate the effect of is to consider further boundary value problems. Thus, we consider in addition three further problems, denoted as specimen , , and , which arise by multiplying the given geometry of the specimen and the imposed boundary conditions with the factors and , respectively. The corresponding graphs of as a function of are displayed in Figure 4 and are referred to as specimens , , and , respectively. Keeping in mind, that is an internal material length, we rescale the abscissa by considering the graphs of as a function of . We expect, that all distributions should coincide, and in fact this is shown in Figure 5.
Altogether, the calculated responses reflect the expected results, which in turn provides a validation of the developed and implemented finite element approximation of the 3-PG Model.
4. Finite element analysis of crack problems
4.1 Edge-cracked specimen
The remaining analysis is referred to the edge-cracked specimen shown in Figure 6a. The assumed width and length of the specimen are, respectively, and , while the crack length is chosen to be . The origin of the Cartesian and the cylindrical coordinate systems, to which we refer, is located at the crack tip.
The specimen is discretized by two meshes (see Figure 6b), a rectangular mesh for the main part of the specimen, consisting of 320 quadratic eight-node elements (with reduced integration) and a radial mesh around the crack tip. The radial mesh, schematically shown in Figure 6c, consists of stripes of quadratic elements with reduced integration, that is, elements with eight nodes, are used. Every strip consists of 45 elements, which decrease in size the closer to the crack tip they are placed. The elements, which contain the crack tip itself, are singular and are the so-called quarter point elements. That means, that one edge of each quadratic element degenerates into a point and two intermediate nodes of two adjacent edges are repositioned. To be more specific, they are moved to a position only a quarter of the edge length away from the crack tip, as indicated in Figure 6d. Due to this repositioning, the corresponding shape functions differ from those of regular 8–node elements. Specifically, the shape functions of quarter point elements consist of linear shape functions extended by a term proportional to the square root of the corresponding component of the position vector (see Henshell and Shaw  or Barsoum ).
For all calculations in the remainder of the chapter, the values of the classical material parameters and are given by Eq. (39) and the following Dirichlet boundary conditions are imposed:
which excludes the possibility of rigid body motions.
The crack faces are subjected to the following boundary conditions for the classical stress: Mode-I near-tip fields are produced by imposing an internal pressure to act, while Mode-II near-tip fields are enforced by subjecting the crack faces to a shear stress loading of . All nonclassical traction components are supposed to vanish for both Mode-I and Mode-II crack types.
4.2 Stress intensity factors
Finite element simulations of Mode-I and Mode-II crack problems allow to verify, numerically, the order of singularity and the angular functions as well as to determinate the stress intensity factors. First, we verify the order of singularity and show how to determine stress intensity factors numerically. Motivated by the analytical solutions in Section “Discussion of the asymptotic solutions” of Part II, assume that the stress and the double stress in the vicinity of the crack tip are given by
These equations indicate, respectively, a linear respone with slope . This means that, as in the case of classical elastic fracture mechanics, one can fit the exponent on the basis of values of stress components calculated by the finite element method at the nodes ahead of the crack tip (see, e.g., Figure 7). Computed responses with the finite element model have confirmed the value with great accuracy. Therefore, the value will be fixed, to avoid errors owing to inaccurate numerical determination of when discussing predicted responses. All stress intensity factors are determined from Eqs. (49)–(53) for the fixed value by applying the least square method. It can be recognized from Figure 7, that the linear response of the stress applies for a radius . The numerical determination of the angular functions below is referred to a fixed radius . The verification of the angular functions and the effect of material parameters on the stress intensity factor is discussed seperately for Mode-I and Mode-II crack problems.
4.3 Results for mode-I
Figures 8–11 display plots of the angular functions , and for both the analytical and the finite element solutions. The plots of the latter are constructed by dividing the values of , , and at a radius by , and , respectively. The stress intensity factors and have been determined by least square fitting as described in the last section. In the finite element computations, the values and for the nonclassical material parameters are chosen. The values of are determined to be , , and . The general observation is that there is good agreement between the analytical and the finite element predictions. The nonvanishing values of verify the existence of this constant terms. The fact that the analytical and the numerical results of the angular functions agree very well verifies the assumed symmetry conditions for (see Section “Symmetry conditions” in Part I).
It is worth remarking that even though the asymptotic solutions of have the same form as the ones of classical elasticity, the corresponding values the stress intensity factors and of the 3-PG Model will be, in general, different from the stress intensity factors and of classical elasticity. This difference results from the fact that the 3-PG Model and the classical elasticity imply, in general, different distributions of the components for identical classical boundary conditions.
Recall also, that the 3-PG Model includes two material parameters more than the classical elasticity, namely, and . The effect of these nonclassical material parameters is illustrated in Figures 12 and 13. These figures reveal that for and very large values of or for and very small values of , the values of converge to the values of . In other words, the responses of the 3-PG Model approach the ones according to classical elasticity.
The effect of and on the nonclassical stress intensity factor is illustrated in Figures 14 and 15. The principal observation is that for smaller values of , the value of gets smaller as well. On the other hand, if , then the values of increase with increasing values of .
4.4 Results for mode-II
The graphs of the angular functions for both the analytical and the finite element solutions are shown in Figure 16. It can be seen that the graphs fit very well.
Since there are two stress intensity factors of Mode-II, , and , we find it convenient, in the vicinity of the crack tip, to consider the angular distributions of and themselves instead of the corresponding distributions of the angular functions , , , and . For , Figures 17–19 illustrate angular distributions of and for both the analytical and the numerical solutions. Once more, we can recognize, that the analytical and the numerical results agree with great accuracy, which also verifies the assumed symmetry conditions of Mode-II for (see Section “Symmetry conditions” in Part I). Further, the values of have been determined to be , , and and this, in turn, verifies the existence of the constant terms .
The effect of and on the nonclassical stress intensity factors and is illustrated in Figures 22–25. Again, this effect is quite similar to the effect of and on , but the stress intensity factors and both are negative in contrast to the positive stress intensity factor .
With regard to Mode-II crack problems, it is of interest to recall that the general solutions of depend on three constants of integration, that is, , , and . By virtue of the boundary conditions on the crack faces, we obtained the relation , which is independent of the crack geometry, that is, the crack length. We do not know further conditions to relate and with each other. However, keeping in mind that, in general, stress intensity factors depend on the crack geometry and the applied loading conditions, one may ask if there exists a relation between and , or equivalently between and , which is independent of the crack length. To clarify this question, we consider the effect of the crack length and the applied shear stress on the ratio . Evidently, this ratio will depend on the material parameters and . The effect of and is illustrated in Figure 26a–d. In all figures, the width and the length of the specimen are the same and equal to the values given in Section 4.1. However, and are different for each figure. It can be recognized from Figure 26a and b that the ratio does not depend on the applied loading, , which might be expected, for the problem is linear. But the ratio depends on the crack length, , as can be seen by comparing Figure 26a and c or Figure 26b and d. This result may be seen as a justification of considering and as different stress intensity factors.
5. Concluding remarks
A finite element model for the 3-PG Model has been developed and verified with reference to the expected responses of a square section with a hole subjected to displacement-controlled tension loading. Then, the finite element model has been employed to discuss Mode-I and Mode-II crack problems. From these investigations, one can draw the following conclusions.
The finite element and the analytical solutions fit with good accuracy.
In particular, the finite element simulations confirm the order of singularity, , and the assumed symmetry conditions for the micro-deformation very well. Further, they confirm the existence of a constant term in the asymptotic expansion of the micro-deformation.
By using numerical simulations with the developed finite element model, the effect of the nonclassical material parameters on the classical and the nonclassical stress intensity factors has been investigated. Especially, limiting cases of material parameters leading to responses of classical elasticity have been analyzed.
For Mode-II loading conditions, two independent stress intensity factors have been assumed to be present in the near-tip fields of the double stress. The numerical investigation seem to verify this remarkable feature.
The authors would like to thank the TU Darmstadt for support of publishing this work in open access.
Each one of the authors, Broese, Frischmann, and Tsakmakis, contributed the same amount of work for the three chapters “Mode-I and Mode-II Crack Tip Fields in Implicit Gradient Elasticity Based on Laplacians of Stress and Strain,” Parts I–III.
Dr. C. Broese: theory and numerical simulations.
Dr. J. Frischmann: analytical solution and numerical simulations.
Prof. Dr. Tsakmakis: theory and analytical solution.
Therefore, it is a joint work by all three authors.
The first and second authors acknowledge with thanks the Deutsche Forschungsgemeinschaft (DFG) for partial support of this work under Grant TS 29/13-1.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this chapter.