- Email: [email protected]
Investigation of finite element mesh independence in rate dependent materials
Computational Materials Science 37 (2006) 442–453 www.elsevier.com/locate/commatsci
Investigation of finite element mesh independence in rate dependent materials F.J. Harewood a, P.E. McHugh
a,b,*
a b
National Centre for Biomedical Engineering Science, National University of Ireland, Galway, Ireland Department of Mechanical and Biomedical Engineering, National University of Ireland, Galway, Ireland Received 10 August 2005; received in revised form 14 November 2005; accepted 14 November 2005
Abstract Crystal plasticity theory is commonly used in finite element analyses to predict large strain ductility in single crystal and polycrystal deformation. In the rate-dependent formulation of the theory it is possible, for cases of simple deformation, to achieve an analytical solution that is independent of any effects due to the finite element mesh spacing. In this study single crystal and polycrystal models were subjected to alternative loading conditions. The effect of the mesh density on the generation of strain localisations and shear bands was investigated with regard to consistency of results. It was found that, prior to the initiation of a narrow shear band, it was possible to achieve a numerical result independent of mesh spacing. In the larger polycrystal analyses, an element size was identified that enabled the generation of a mesh independent solution. This allowed the accurate prediction of the mechanical behaviour of the model up to, and including, the failure point. The implications of this for small-scale metallic device design are discussed. 2005 Elsevier B.V. All rights reserved. Keywords: Micromechanical modelling; Polycrystalline material; Crystal plasticity theory; Plastic strain localisation; Finite element method
1. Introduction A fundamental issue of importance in using the finite element method is the effect that the number of elements has on the solution of a problem. A compromise must be made between the high level of resolution and accuracy delivered by a large number of elements and the increased computational expense necessary to perform a large simulation. The user must determine the level of resolution needed and balance this with the available computational resources. In regions of high solution gradient, unless there is a sufficient number of elements, the accuracy of the results can be compromised. Experimental testing of metallic single crystals and polycrystals has shown that when the sample is subjected to * Corresponding author. Address: National Centre for Biomedical Engineering Science, National University of Ireland, Galway, Ireland. Tel.: +353 (0) 91 524411x3152. E-mail address: [email protected] (P.E. McHugh).
0927-0256/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2005.11.004
large strains, the initial relatively homogeneous mechanical response breaks down and strain localisations are observed [6,15,20–22,24]. These localised deformations commonly coalesce to form shear bands that traverse the metallic grains. Shear band localisations determine the location of, and are a precursor to, ductile failure in a metal. Finite element studies have been carried out that quantify the mechanical performance of single metallic crystals e.g. [4,5,7,15,18]. Crystal plasticity theory is often used for such analyses and can be used to yield accurate predictions of large strain ductility and failure in boundary value problems [6,8,19]. The various formulations of the theory can be divided into two broad categories; rate-independent and rate-dependent. In large strain crystalline finite element analyses, the presence of shear bands can lead to a high level of mesh dependence in the rate-independent formulation, whereas a solution independent of mesh spacing can be achieved using the rate-dependent formulations [14]. Several authors have examined how the presence of crystallographic orientation gradients, i.e. the variation in
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
crystallographic orientation across the grain due to inhomogeneous deformation, in plastically-strained viscoplastic materials, affects the local deformation behaviour [6,18,22,24]. This has also been investigated using the finite element method [4,17]. Buchheit et al. [4] carried out a study on the effect of mesh density on the orientation gradients in a polycrystal. Grains were classified into two groups based on their initial orientations; intrinsically stable and intrinsically unstable. The former group did not develop orientation gradients during straining while the latter group experienced texture development. In these cases a high level of mesh refinement was required to capture such gradients. These studies show that the accuracy of the single crystal formulation is reliant on the ability of a mesh to resolve a gradient in the material. In this study a single crystal compression analysis is carried out in order to generate strain localisations and a shear band in the material. The influence that the mesh density has on the macroscopic mechanical response and the ingrain contour profile is investigated. The effect that the level of rate sensitivity of the material has on the formation of shear bands is also examined. Following this, the study is extended to include the formation of deformation localisations in a larger two-dimensional model of a 75 lm thick wire. The model is tensed beyond the ultimate tensile strength and the layout of the shear bands is considered and linked to the smaller single grain analyses. Larger 2D models such as these have been shown to yield accurate indications of failure in the design of small scale metallic devices [19]. To maximise the efficiency of these analyses, it is important to define the limits of accuracy in the finite element models used. The ultimate goal of the study is to generate a finite element model that, when strained, is unaffected and independent of the mesh density and can be used to model the mechanical performance of small-scale metallic objects where crystal plasticity is an appropriate constitutive formulation. The following two Sub-sections (1.1 and 1.2) detail the key differences in the way shear bands are dealt with and the numerical complications involved when the rate-dependent and the rate-independent formulations are used. In Section 2 the crystal plasticity theory, as implemented in the finite element solver ABAQUS/Standard, is outlined. Sections 3 and 4 describe the two types of finite element simulations and their results. Finally the implications of the results are discussed.
443
a length scale, the element size, in the set of finite element equations. The solution typically seen is the narrowest band that the size of finite element can resolve. This is inherently non-unique and has little or no physical meaning. To deal with this problem improvements on the formulation have been introduced [1,2,7]. Garikipati and Hughes [7] developed a continuum theory where the calculations are performed on a multiple scale framework. Areas of high strain gradient are solved on a fine scale while other areas of the model are solved on a coarser scale. De Borst et al. [5] examined several different types of continuum models including the Cosserat continuum. In this formulation additional couple stresses and microcurvatures of the material were added. They have the effect of adding a length scale to the formulation that determines the thickness of shear bands. De Borst et al. found that the ratedependent formulation was the most appealing as it did not suffer from the severe mesh dependence that affects the rate-independent case. 1.2. Rate-dependent material In a rate-dependent formulation the formation of strain gradients and bands is due to a material imperfection or a geometrical inhomogeneity. Assuming that the material in question contains no impurities it can be said that when a shear band occurs it is the geometry that introduces a length scale to the analysis [5,14–16]. This is in contrast to the rate-independent case where it is the mesh spacing that provides the governing length scale. Needleman [14] presents a discussion on the occurrence of strain localisations in a simple shearing problem where only one stress component is active. A block with an impurity in the form of a band of material with different properties to the surrounding material is subjected to simple shear. In the case of rate-dependent materials, an imperfection, be it geometrical or material, acts as a length scale in influencing the form of a shear band. This becomes a governing length scale in the problem. Therefore it was possible to generate results from a finite element model where the width of the shear band was independent of the element size used. In the present paper the mesh sensitivity for alternative loading modes and in the context of a representative microscale model of a metallic structure is explored. 2. Theory
1.1. Rate-independent material In the case of classical strain rate-independent continua the existence of high strain gradients can only be satisfied when the stiffness tensor loses ellipticity [3,14,16]. This gives rise to an arbitrary choice of solution. The formulation contains no intrinsic length scale to determine the shear band width. The finite element solution of strain localisations in a rate-independent material proves to have chronic mesh sensitivity. This is due to the introduction of
The finite element analyses performed in this study incorporate elastic and plastic constitutive laws in the context of finite deformation kinematics. The elasticity is considered to be isotropic and linear in terms of finite deformation quantities and can be described using the Young’s modulus, E, and Poisson’s ratio, m. This is a reasonable approach as the elastic strains in the model are very small in comparison to the plastic strains. Plasticity is described using a rate-dependent crystal plasticity theory.
444
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
Crystal plasticity theory is a physically based plasticity theory that represents the deformation of a metal at the microscale. The flow of dislocations in a metallic crystal along slip systems is represented in a continuum framework. Plastic strain is assumed to be solely due to crystallographic dislocation slip. The rate-dependent, viscoplastic single crystal theory, as presented in Huang [10] and Peirce et al. [16] and employed in several works such as McHugh and Connolly [12], Savage et al. [19] and McGarry et al. [11], is used. The rate-dependence is implemented in the formulation through a power-law that relates the resolved shear stress, s(a), to the slipping rate, c_ ðaÞ , on each slip system a ðaÞ n s ðaÞ ðaÞ ð2:1Þ c_ ¼ a_ sgnðs Þ ðaÞ g where a_ and n are a reference strain rate and rate sensitivity exponent, respectively, and g(a) is the slip system strain hardness. As n tends to infinity the material reaches the rate-independent limit. Slip system strain hardening is specified by the following hardness function as defined by Peirce et al. [16]: h0 c a ð2:2Þ gðca Þ ¼ g0 þ ðg1 g0 Þh0 tanh g 1 g0 where the material hardening parameters g0, g1 and h0 are the initial hardness of a slip system, the maximum slip system hardness when og/oca = 0, and the value of og/oca when ca = 0, respectively. These material constants can be derived from the strain hardening part of an experimental tensile stress–strain (r–e) curve for a material. The values used in this study were taken from the tensile data of a 316L stainless steel wire [23]: g0 ¼ 150 MPa a_ ¼ 0:0106 s
g1 ¼ 360 MPa
3. Single grain compression 3.1. Model description The objective of the first set of simulations is to investigate the effect of mesh refinement on shear banding and whether it is possible to reproduce a scenario where the solution is independent of element size, as has been seen in simple shearing cases [14]. A shear band is initiated in a 2D geometry similar to that presented by Harren et al. [8]. As stated in the introduction, an inhomogeneity is only possible when a material or geometrical imperfection is present in the model. The 2D crystal model has a geometrical imperfection along the right hand edge. It is described by x ¼ 1 þ ð5 103 Þ cosðpyÞ
ð3:1Þ
with x and y being the spatial coordinates (Fig. 1). The bottom left hand corner of the block is constrained in all directions while the nodes along the bottom edge are free to move in the x-direction only. The model is subjected to a compressive engineering strain (eeng) of 0.5 in the ydirection. The engineering strain is defined as the ratio of the increase in height to the original height: (h1 h0)/h0, where h1 is the compressed height and h0 is the original height. A rate sensitivity exponent of n = 50 is used. This gives a moderate rate-dependent response and is a physically representative value for the material being modelled [23].
h0 ¼ 380 MPa
1
Several values of n are used. A quantity of importance in the above equation is the accumulated slip, ca. This is a measure of the total crystallographic plastic strain at a point and is defined as follows: X Z t c_ ðaÞ dt ca ¼ ð2:3Þ a
0
where t is the time and the summation is over all the slip systems at a point. The material chosen, 316L stainless steel, has a face-centred-cubic crystalline structure. Slip is assumed to occur along specific slip directions as defined by its lattice and the slip system {1 1 1}h1 1 0i. The above constitutive theory is implemented in the finite element code ABAQUS/Standard as a user-defined subroutine (UMAT), as coded in [10]. This subroutine, written in Fortran, implements the theory in the form of a stress update algorithm that is called at each integration point in the mesh and for every iteration during a finite element simulation [9].
Fig. 1. Single grain model—40 · 80 mesh. Units are in mm.
Fig. 2. Cube depicting slip direction and crystal orientation for the single crystal compression test, a 54.7. All slip directions are parallel to the three shown.
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
445
crystal axes, with respect to the global axes, are specified as material properties in the user interface. A crystal orientation is used whereby the [1 0 1] and the ½1 0 1 crystallographic directions are parallel to the global x and z-axes, respectively. In this arrangement, in a two-dimensional analysis, several slip directions become coincident (Fig. 2). Thus, twelve possible slip directions are reduced to three; two oblique and the horizontal direction. The angle subtended by the oblique and the horizontal slip directions in 2D is a 54.7. In this arrangement the directions of slip are assured to be oblique in the model and the formation of slip bands is promoted. Fig. 3. Force–displacement curves for 10 · 20, 20 · 40, 40 · 80 and 100 · 200 meshes for the rate sensitivity parameter, n = 50. Point A indicates a strain of eeng = 0.38 which is the upper limit for generation of mesh independent results.
The model is meshed with eight-noded generalised plane strain quadrilateral elements (CPEG8R), each having four (2 · 2) integration points. Four increasingly fine mesh densities are used; models with 10 · 20, 20 · 40, 40 · 80, and 100 · 200 elements are generated. The orientations of the
3.2. Results The force applied to the top edge of the model against the vertical displacement of the edge for all four meshes is graphed in Fig. 3. The four curves compare well at lower levels of compression. One would be tempted to say that both convergence and a result independent of element size had been achieved. However the two denser meshes show a softening that is due to the formation of a sharp shear band across the thickness of the material. This softening can be considered to represent failure of the material in this case.
Fig. 4. Contour plots of the accumulated shear strain in the 40 · 80 mesh at increasing levels of compression: (a) eeng = 0.18, (b) 0.32, (c) 0.5.
446
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
It is accompanied by considerable mesh dependence where the shear band thickness is determined by the element size. The contour plots in this paper show the accumulated slip, ca (Eq. (2.3)). This value is a measure of the crystallographic plastic strain and is analogous, but not equal, to the equivalent plastic strain. Fig. 4(a)–(c) shows the 40 · 80 element mesh of the single grain at increasing levels of strain. At eeng = 0.3 (Fig. 4(b)) the localised areas of shearing form bands at regular intervals along the length
of the grain, running at oblique angles to the direction of straining. The two shear bands have similar widths and profiles. As the strain level increases the shear bands become narrower and the strain gradients become steeper, with the band on the right hand side becoming particularly localised (Fig. 4(c)). A visual comparison of the four meshes at eeng = 0.5 shows that, as the shear band becomes narrow, there is a high level of mesh sensitivity (Fig. 5(a)–(d)). The thickness of the shear band, in terms
Fig. 5. Contour plots (with mesh removed for clarity) of the accumulated shear strain at eeng = 0.5: (a) 10 · 20 mesh, (b) 20 · 40 mesh, (c) 40 · 80 mesh, (d) 100 · 200 mesh, (e) 120 · 240 mesh.
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
Fig. 6. Close-up of the shear band in the 40 · 80 mesh at eeng = 0.5. The highlighted elements indicate that the width of the shear band is set by the span of one element.
of element size, is difficult to assess as the mesh is highly distorted within the band. It can be said that there is only one element completely within the shear band and it spans the width of the band, as illustrated in Fig. 6. This is the case with the two denser meshes that exhibit material failure (Fig. 3). It is interesting to note the formation of shear bands at the right hand edge of the model. The 100 · 200 mesh shows two strong shear bands and the initiation of a third one at the top corner (Fig. 5(d)). All three bands are caused by the compression of the material along the concave part of the geometrical imperfection. The gradients of strain at the same location in the 40 · 80 mesh are not as sharp (Fig. 5(c)). It is noticeable that the increase in mesh density merely improves the resolution of the solution, it does not increase the number of shear bands present. Whereas, the 20 · 40 mesh has only two shear bands along the right hand edge (Fig. 5(b)). It would appear that the multiple shear bands are part of the numerical solution to the problem and that the number of shear bands is independent of
447
mesh density above a certain level, i.e. the initiation of shear bands on the right hand edge does not appear to be a function of the number of elements used. It is estimated that a mesh density higher than the 100 · 200 case would not yield an increase in the number of shear bands but would cause narrower shear bands. In order to illustrate this point, a more dense mesh with 120 · 240 elements was constructed. The resulting contour plot (Fig. 5(e)) shows that a higher mesh density does not yield a higher number of shear bands than the 100 · 200 case. It is reasonable to assume that a limit, in terms of shear banding, has been reached. 3.3. Rate sensitivity To assess the influence of the strain rate sensitivity parameter, n, on shear band formation, simulations are performed using the above models where n is varied. More rate-dependent values of n = 10 and 20 and the quite rateindependent value, n = 100, are used. The three cases are modelled using the 20 · 40 and 40 · 80 meshes. It is known that material rate sensitivity delays the development of shear bands [14,16]. This is shown where the increasing strain rate sensitivity blunts the shear strain gradient. The resulting deformed mesh for n = 20 shows a much more diffuse shear band running northwest to southeast (Fig. 7). Unlike in the original case, the denser mesh does not produce a narrow shear band whose width is dependent on the element size. There is a better agreement between the resulting contours in the two models. Although, the profiles of the localised strains are not converged, the macroscopic results are not sensitive to these differences and compare favourably (Fig. 8). When the rate sensitivity is increased further, with n = 10, the strain localisations become even more diffuse (Fig. 9). As a result the denser mesh does not identify
Fig. 7. Contour plots of the accumulated shear strain for rate sensitivity value n = 20 at eeng = 0.5: (a) 20 · 40 mesh, (b) 40 · 80 mesh.
448
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
Fig. 8. Force–displacement curves for 20 · 40 and 40 · 80 meshes for the rate sensitivity parameter, n = 20.
any gradients that are not present in the sparser one. The contour plots and the force–displacement curves (Fig. 10) compare well for the two mesh densities and can be said to have converged. When a relatively rate-independent value is used, n = 100, the formation of narrow shear bands is quite evident. As would be expected, neither the applied force– displacement responses (Fig. 11) nor the contour plots (Fig. 12) compare well. This is an indication of the mesh sensitivity of results that is inherent in the rate-independent finite element analyses. The above simulations show that as well as being deformation dependent, the mesh sensitivity of a model is also heavily reliant on the degree of rate sensitivity of the material.
Fig. 10. Force–displacement curves for 20 · 40 and 40 · 80 meshes for the rate sensitivity parameter, n = 10.
Fig. 11. Force–displacement curves for 20 · 40 and 40 · 80 meshes for the rate sensitivity parameter, n = 100.
3.4. Discussion The finite element analysis of shear band formation in a rate-dependent material is shown to be heavily reliant on two factors. The first of these factors is the level of rate sen-
sitivity of the material and the second is the amount of deformation of the model. Increasing the rate sensitivity has the effect of delaying the onset and sharpening of the
Fig. 9. Contour plots of the accumulated shear strain for rate sensitivity value n = 10 at eeng = 0.5: (a) 20 · 40 mesh, (b) 40 · 80 mesh.
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
449
Fig. 12. Contour plots of the accumulated shear strain for rate sensitivity value n = 100 at eeng = 0.5: (a) 20 · 40 mesh, (b) 40 · 80 mesh.
strain gradient. Thus the model does not require the high level of resolution that a dense mesh provides. The sharp strain gradient present in the denser meshes and the more rate insensitive materials is a highly mesh sensitive state. The width of the shear band is such that one element spans it, regardless of mesh density. This is caused by excessive element distortion and is accompanied by the loss of load-bearing capability in the model. This can be viewed as failure of the material. Although several values of the rate sensitivity parameter were used for the same type of analysis, this is not a possibility for practical simulations when an accurate material description must be used. Therefore, the mesh density must be adjusted in order to yield mesh independent results. For the n = 20 compression case presented in Section 3.3 there is a good level of agreement in terms of the contour plots (Fig. 7) and the macroscopic mechanical responses (Fig. 8). Although the profiles of ca for the two meshes are not identical this is not necessary to get a converged, correct solution for the macroscopic behaviour. A similar situation occurs in the n = 50 case (Section 3.2) where there is a good level of agreement up to a strain level of eeng = 0.38, as indicated by point A in Fig. 3. We can infer from this that a solution can be generated that is effectively independent of the element size used so long as the
overall deformation of the model falls below a certain strain limit (eeng < 0.38, in this case). 4. Polycrystal tension 4.1. Model description The stress–strain behaviour of small scale 316L stainless steel devices is dependent on their size. This size dependence affects, in particular, the point of failure [13]. Polycrystalline finite element models incorporating a crystal plasticity constitutive subroutine have been shown to capture this size effect successfully [19]. In order to investigate the effect of mesh refinement in larger models a 2D polycrystalline model of a 316L stainless steel wire of 75 lm thickness is constructed. The model is based on that presented in Savage et al. [19]. The granular microstructure is explicitly represented with each grain idealised as a hexagon of area 92 lm2 (Fig. 13). As outlined in Section 2, the crystal plasticity constitutive model used is implemented via a user material subroutine (UMAT) in ABAQUS. Although the model is based on a 2D geometry each grain has a full 3D crystallographic description that is randomly oriented for each grain. In addition, the use of 8-noded generalised plane strain elements (CPEG8R) allows for
Fig. 13. Model of 316L stainless steel wire with a thickness of 75 lm. The granular microstructure is idealised as a series of hexagons (highlighted).
450
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
Fig. 14. Macroscopic stress–strain curves for all the meshes analysed in orientation set 2. The legend refers to the number of elements per grain in each mesh. Inset: stress–strain curves showing an agreement in the failure strains for the three denser meshes.
uniform out-of-plane deformation. Therefore the model is considered ‘‘quasi-3D’’. A rate sensitivity parameter of n = 50 is used throughout. To achieve a reasonable spread of data, three sets of randomly assigned orientations are used. Several different mesh resolutions are analysed. Mesh densities of 8, 50, 98, 162 and 242 elements per full grain are constructed. The left hand bottom corner of the model is constrained in all degrees of freedom and the nodes along the left hand edge are free to move in the y-direction only. The wire is tensed by applying a uniform horizontal (x-direction) displacement to the right hand edge of the model. The engineering strain is defined as the ratio of the increase in length to the original length: (l1 l0)/l0, where l1 is the stretched length and l0 is the original length. The engineering stress is the ratio of force applied to the wire to the original cross-sectional area of the wire: F/A0. This produces engineering stress–engineering strain (reng–eeng) curves and is a convenient method for presenting tensile data. In these simulations the failure strain is taken as the strain at the point of ultimate tensile stress (UTS). 4.2. Results Each grain displays a different hardness in the direction of tension due to the random crystal orientation. The softer grains act as initiation sites for strain localisations. This results in non-uniform deformation and coalescence of areas of plastic strain, eventually forming a shear band that
traverses the width of the wire. Although an explicit material separation model, representing fracture, is not included in the current crystal plasticity formulation, the shear band indicates the site of necking and eventual failure of the wire. It is important to note that pockets of localised shearing occur at many places throughout the wire but eventual necking occurs at the weakest part, i.e. where the grain lattice orientations are most favourable for plastic slip. Following the failure strain point (UTS), localisations continue to grow and coalesce leading to macroscopic softening. The results presented are for the second set of random orientations assigned to the grains and are typical of all three sets. The macroscopic reng–eeng curves, comparing the response of the different mesh densities, are shown in Fig. 14. The three more densely meshed geometries (98, 162 and 242 elements per grain) show good agreement and have converged in terms of failure strain (Fig. 14). The hardness of the wire reduces as the number of elements in the model increases. This is due to the increase in number of degrees of freedom in the model. It is also important to investigate the contours of the localisations present in the granular microstructure of the material. Fig. 15 shows the profile of ca at eeng = 0.6 from one of the lattice orientation sets. The shear localisations lie in bands at oblique angles to the tensile axis. To illustrate a convergence in the results a single grain that lies within the shear band is monitored throughout the analysis for the different meshes. The contours of ca in a single grain
Fig. 15. Accumulated shear strain during necking of a 75 lm wire. Shear localisations form at angles oblique to the tensile axis.
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
451
Fig. 16. Accumulated shear strain, ca, in a grain at the point of failure, eeng = 0.35. The grain lies within the shear band that traverses the wire. The contour limits of each plot are the same.
at the point of failure are shown in Fig. 16 for all five mesh densities. The profiles of the shear band in the three denser meshes are quite similar. As the straining continues the gradients become steeper and narrower. The same grain is examined at a strain of eeng = 0.44; Fig. 17(a)–(e) shows ca for all five meshes. The shear band profiles do not appear to agree as well as in Fig. 16. In this case an increase in number of elements continues to yield an improvement in the results and further mesh refinement would be required to achieve convergence. As each grain undergoes higher deformation the strain gradients become more severe and a smaller element size is needed. This is similar to the single grain compression simulation. In contrast to the single grain compression case the point of failure, as identified in the macroscopic response, does not correspond to development of a very narrow, element size dependent shear band. 4.3. Discussion Whereas in Section 3 analyses are carried out that introduce shear bands through a geometrical imperfection, in
this section the cause of the strain localisations is the material inhomogeneity. Due to the microstructural arrangement in the wire, each grain experiences non-uniform deformation under large straining. Each grain can be viewed as an isolated block of material that is subjected to non-uniform boundary deformation. Unlike the single crystal compression case the point of macroscopic failure does not correspond to the initiation of a narrow shear band. Here, the softening of the material is due to the strains in the softer grains becoming very large in comparison to the harder grains and leading to a reduction in the cross-sectional area of the wire. The favourable comparison of the reng–eeng responses in the three more densely meshed grains up to the point of UTS (Fig. 14) is a similar case to that seen in Fig. 3. Although the contour plots of ca are not the same (Fig. 16), the macroscopic responses have converged in terms of the failure point (Fig. 14). In the polycrystal tension simulation the wire is not tensed beyond a strain of eeng = 0.6 (Fig. 15). Were the wire to be stretched further one would see the formation of a shear band with a thickness of one element as is seen in the single grain analyses. Prior to this situation it is possible
452
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
Fig. 17. Accumulated shear strain, ca, in a grain after the point of failure, eeng = 0.42. The grain lies within the shear band that traverses the wire. The contour limits of each plot are the same.
to get a solution that is accurate and unaffected by the element size used. In the current analyses a mesh density of 98 elements per grain has been shown to provide a sufficient level of accuracy in order to provide the failure point. This translates to a mesh to grain area ratio of 0.01. Although the in-grain contour profiles of ca are not completely converged there is still a good agreement between this case and the subsequent more finely meshed models (Fig. 16(c)–(e)). The more persuasive argument for choosing this mesh spacing is by examining the macroscopic performances (Fig. 14). The model with the greater number of elements continues to yield an overall reduction in the stiffness of the wire. This is due to the increased number of degrees of freedom in the mesh. Whereas the strain at which the UTS occurs does not change for the three more finely meshed models. 5. Conclusion This paper shows that computational models that are capable of producing results independent of the element size can be developed. In single crystal compression analyses, excessive straining causes a situation where shear
bands show considerable mesh dependence. Prior to this situation a converged solution was found in both macroscopic and in-grain mechanical performance. The study was then extended to include the tensile behaviour of a 75 lm thickness wire. The microstructure of the model was explicitly represented and the material properties were chosen to approximate 316L stainless steel behaviour. It was found that the point of failure, taken as the strain at the UTS, did not correspond to a situation of critical mesh dependence. A mesh density per grain could be identified that gave a relatively accurate profile of the intragranular plastic strains and a converged value for the macroscopic failure strain. The accurate prediction of failure is vital in the design of small scale metallic components such as cardiovascular stents in the medical device industry or components in the microelectronics industry. The results of the work show that, for the application considered here, by employing an element area to grain area ratio of 0.01 there is an acceptable ‘window’ of macroscopic strain that can be applied to the model while maintaining accuracy of solution. This window of strain extends to the point of failure. This mesh density can be used for further analysis to investigate the
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
mechanical performance of small-scale metallic devices under various loading conditions, with confidence in the level of accuracy provided by the results. Acknowledgements The authors wish to acknowledge funding from Embark, Irish Research Council for Science, Engineering and Technology: Funded by the National Development Plan. The simulations in this work were performed on the SGI 3800 high performance computer at NUI, Galway. References [1] L. Anand, M. Kothari, Journal of the Mechanics and Physics of Solids 44 (4) (1996) 525–558. [2] D.J. Bammann, Materials Science and Engineering A 309–310 (2001) 406–410. [3] T. Belytschko, Z.P. Bazant, Yul-Woong Hyun, Ta-Peng Chang, Computers and Structures 23 (2) (1986) 163–180. [4] T.E. Buchheit, G.W. Wellman, C.C. Battaile, International Journal of Plasticity 21 (2) (2005) 221–249. [5] R. De Borst, L.J. Sluys, H.-B. Muhlhaus, J. Pamin, Engineering Computations 10 (1993) 99–121. [6] F. Delaire, J.L. Raphanel, C. Rey, Acta Materialia 48 (5) (2000) 1075–1087. [7] K. Garikipati, T.J.R. Hughes, Computer Methods in Applied Mechanics and Engineering 159 (3–4) (1998) 193–222. [8] S.V. Harren, H.E. Deve, R.J. Asaro, Acta Metallurgica 36 (9) (1988) 2435–2480.
453
[9] Hibbitt, Karlsson, Sorenson, ABAQUS Theory Manual, Pawtucket, RI, USA, 1997. [10] Y. Huang, A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program, Harvard University Report, MECH 178, 1991. [11] J.P. McGarry, B.P. O’Donnell, P.E. McHugh, J.G. McGarry, Computational Materials Science 31 (3–4) (2004) 421–438. [12] P.E. McHugh, P.J. Connolly, Computational Materials Science 27 (4) (2003) 423–436. [13] B.P. Murphy, P. Savage, P.E. McHugh, D.F. Quinn, Annals of Biomedical Engineering 31 (6) (2003) 686–691. [14] A. Needleman, Computer Methods in Applied Mechanics and Engineering 67 (1) (1988) 69–85. [15] A. Needleman, V. Tvergaard, Element Analysis of Localization in Plasticity, in: Finite Elements—Special problems in Solid Mechanics, Prentice-Hall, New York, 1983, pp. 94–157. [16] D. Peirce, R. Asaro, A. Needleman, Acta Metallurgica et Materialia 31 (12) (1983) 1951–1976. [17] D. Raabe, Z. Zhao, S.-J. Park, F. Roters, Acta Materialia 50 (2) (2002) 421–440. [18] M. Sachtleber, Z. Zhao, D. Raabe, Materials Science and Engineering A 336 (1–2) (2002) 81–87. [19] P. Savage, B.P. O’Donnell, P.E. McHugh, B.P. Murphy, D.F. Quinn, Annals of Biomedical Engineering 32 (2) (2004) 202–211. [20] B.M. Schroeter, D.L. McDowell, International Journal of Plasticity 19 (9) (2003) 1355–1376. [21] A. Staroselsky, L. Anand, Journal of the Mechanics and Physics of Solids 46 (4) (1998) 671–696. [22] A. Tatschl, O. Kolednik, Materials Science and Engineering A 356 (1– 2) (2003) 447–463. [23] X. You, in preparation. [24] N. Zhang, W. Tong, International Journal of Plasticity 20 (3) (2004) 523–542.
Investigation of finite element mesh independence in rate dependent materials F.J. Harewood a, P.E. McHugh
a,b,*
a b
National Centre for Biomedical Engineering Science, National University of Ireland, Galway, Ireland Department of Mechanical and Biomedical Engineering, National University of Ireland, Galway, Ireland Received 10 August 2005; received in revised form 14 November 2005; accepted 14 November 2005
Abstract Crystal plasticity theory is commonly used in finite element analyses to predict large strain ductility in single crystal and polycrystal deformation. In the rate-dependent formulation of the theory it is possible, for cases of simple deformation, to achieve an analytical solution that is independent of any effects due to the finite element mesh spacing. In this study single crystal and polycrystal models were subjected to alternative loading conditions. The effect of the mesh density on the generation of strain localisations and shear bands was investigated with regard to consistency of results. It was found that, prior to the initiation of a narrow shear band, it was possible to achieve a numerical result independent of mesh spacing. In the larger polycrystal analyses, an element size was identified that enabled the generation of a mesh independent solution. This allowed the accurate prediction of the mechanical behaviour of the model up to, and including, the failure point. The implications of this for small-scale metallic device design are discussed. 2005 Elsevier B.V. All rights reserved. Keywords: Micromechanical modelling; Polycrystalline material; Crystal plasticity theory; Plastic strain localisation; Finite element method
1. Introduction A fundamental issue of importance in using the finite element method is the effect that the number of elements has on the solution of a problem. A compromise must be made between the high level of resolution and accuracy delivered by a large number of elements and the increased computational expense necessary to perform a large simulation. The user must determine the level of resolution needed and balance this with the available computational resources. In regions of high solution gradient, unless there is a sufficient number of elements, the accuracy of the results can be compromised. Experimental testing of metallic single crystals and polycrystals has shown that when the sample is subjected to * Corresponding author. Address: National Centre for Biomedical Engineering Science, National University of Ireland, Galway, Ireland. Tel.: +353 (0) 91 524411x3152. E-mail address: [email protected] (P.E. McHugh).
0927-0256/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2005.11.004
large strains, the initial relatively homogeneous mechanical response breaks down and strain localisations are observed [6,15,20–22,24]. These localised deformations commonly coalesce to form shear bands that traverse the metallic grains. Shear band localisations determine the location of, and are a precursor to, ductile failure in a metal. Finite element studies have been carried out that quantify the mechanical performance of single metallic crystals e.g. [4,5,7,15,18]. Crystal plasticity theory is often used for such analyses and can be used to yield accurate predictions of large strain ductility and failure in boundary value problems [6,8,19]. The various formulations of the theory can be divided into two broad categories; rate-independent and rate-dependent. In large strain crystalline finite element analyses, the presence of shear bands can lead to a high level of mesh dependence in the rate-independent formulation, whereas a solution independent of mesh spacing can be achieved using the rate-dependent formulations [14]. Several authors have examined how the presence of crystallographic orientation gradients, i.e. the variation in
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
crystallographic orientation across the grain due to inhomogeneous deformation, in plastically-strained viscoplastic materials, affects the local deformation behaviour [6,18,22,24]. This has also been investigated using the finite element method [4,17]. Buchheit et al. [4] carried out a study on the effect of mesh density on the orientation gradients in a polycrystal. Grains were classified into two groups based on their initial orientations; intrinsically stable and intrinsically unstable. The former group did not develop orientation gradients during straining while the latter group experienced texture development. In these cases a high level of mesh refinement was required to capture such gradients. These studies show that the accuracy of the single crystal formulation is reliant on the ability of a mesh to resolve a gradient in the material. In this study a single crystal compression analysis is carried out in order to generate strain localisations and a shear band in the material. The influence that the mesh density has on the macroscopic mechanical response and the ingrain contour profile is investigated. The effect that the level of rate sensitivity of the material has on the formation of shear bands is also examined. Following this, the study is extended to include the formation of deformation localisations in a larger two-dimensional model of a 75 lm thick wire. The model is tensed beyond the ultimate tensile strength and the layout of the shear bands is considered and linked to the smaller single grain analyses. Larger 2D models such as these have been shown to yield accurate indications of failure in the design of small scale metallic devices [19]. To maximise the efficiency of these analyses, it is important to define the limits of accuracy in the finite element models used. The ultimate goal of the study is to generate a finite element model that, when strained, is unaffected and independent of the mesh density and can be used to model the mechanical performance of small-scale metallic objects where crystal plasticity is an appropriate constitutive formulation. The following two Sub-sections (1.1 and 1.2) detail the key differences in the way shear bands are dealt with and the numerical complications involved when the rate-dependent and the rate-independent formulations are used. In Section 2 the crystal plasticity theory, as implemented in the finite element solver ABAQUS/Standard, is outlined. Sections 3 and 4 describe the two types of finite element simulations and their results. Finally the implications of the results are discussed.
443
a length scale, the element size, in the set of finite element equations. The solution typically seen is the narrowest band that the size of finite element can resolve. This is inherently non-unique and has little or no physical meaning. To deal with this problem improvements on the formulation have been introduced [1,2,7]. Garikipati and Hughes [7] developed a continuum theory where the calculations are performed on a multiple scale framework. Areas of high strain gradient are solved on a fine scale while other areas of the model are solved on a coarser scale. De Borst et al. [5] examined several different types of continuum models including the Cosserat continuum. In this formulation additional couple stresses and microcurvatures of the material were added. They have the effect of adding a length scale to the formulation that determines the thickness of shear bands. De Borst et al. found that the ratedependent formulation was the most appealing as it did not suffer from the severe mesh dependence that affects the rate-independent case. 1.2. Rate-dependent material In a rate-dependent formulation the formation of strain gradients and bands is due to a material imperfection or a geometrical inhomogeneity. Assuming that the material in question contains no impurities it can be said that when a shear band occurs it is the geometry that introduces a length scale to the analysis [5,14–16]. This is in contrast to the rate-independent case where it is the mesh spacing that provides the governing length scale. Needleman [14] presents a discussion on the occurrence of strain localisations in a simple shearing problem where only one stress component is active. A block with an impurity in the form of a band of material with different properties to the surrounding material is subjected to simple shear. In the case of rate-dependent materials, an imperfection, be it geometrical or material, acts as a length scale in influencing the form of a shear band. This becomes a governing length scale in the problem. Therefore it was possible to generate results from a finite element model where the width of the shear band was independent of the element size used. In the present paper the mesh sensitivity for alternative loading modes and in the context of a representative microscale model of a metallic structure is explored. 2. Theory
1.1. Rate-independent material In the case of classical strain rate-independent continua the existence of high strain gradients can only be satisfied when the stiffness tensor loses ellipticity [3,14,16]. This gives rise to an arbitrary choice of solution. The formulation contains no intrinsic length scale to determine the shear band width. The finite element solution of strain localisations in a rate-independent material proves to have chronic mesh sensitivity. This is due to the introduction of
The finite element analyses performed in this study incorporate elastic and plastic constitutive laws in the context of finite deformation kinematics. The elasticity is considered to be isotropic and linear in terms of finite deformation quantities and can be described using the Young’s modulus, E, and Poisson’s ratio, m. This is a reasonable approach as the elastic strains in the model are very small in comparison to the plastic strains. Plasticity is described using a rate-dependent crystal plasticity theory.
444
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
Crystal plasticity theory is a physically based plasticity theory that represents the deformation of a metal at the microscale. The flow of dislocations in a metallic crystal along slip systems is represented in a continuum framework. Plastic strain is assumed to be solely due to crystallographic dislocation slip. The rate-dependent, viscoplastic single crystal theory, as presented in Huang [10] and Peirce et al. [16] and employed in several works such as McHugh and Connolly [12], Savage et al. [19] and McGarry et al. [11], is used. The rate-dependence is implemented in the formulation through a power-law that relates the resolved shear stress, s(a), to the slipping rate, c_ ðaÞ , on each slip system a ðaÞ n s ðaÞ ðaÞ ð2:1Þ c_ ¼ a_ sgnðs Þ ðaÞ g where a_ and n are a reference strain rate and rate sensitivity exponent, respectively, and g(a) is the slip system strain hardness. As n tends to infinity the material reaches the rate-independent limit. Slip system strain hardening is specified by the following hardness function as defined by Peirce et al. [16]: h0 c a ð2:2Þ gðca Þ ¼ g0 þ ðg1 g0 Þh0 tanh g 1 g0 where the material hardening parameters g0, g1 and h0 are the initial hardness of a slip system, the maximum slip system hardness when og/oca = 0, and the value of og/oca when ca = 0, respectively. These material constants can be derived from the strain hardening part of an experimental tensile stress–strain (r–e) curve for a material. The values used in this study were taken from the tensile data of a 316L stainless steel wire [23]: g0 ¼ 150 MPa a_ ¼ 0:0106 s
g1 ¼ 360 MPa
3. Single grain compression 3.1. Model description The objective of the first set of simulations is to investigate the effect of mesh refinement on shear banding and whether it is possible to reproduce a scenario where the solution is independent of element size, as has been seen in simple shearing cases [14]. A shear band is initiated in a 2D geometry similar to that presented by Harren et al. [8]. As stated in the introduction, an inhomogeneity is only possible when a material or geometrical imperfection is present in the model. The 2D crystal model has a geometrical imperfection along the right hand edge. It is described by x ¼ 1 þ ð5 103 Þ cosðpyÞ
ð3:1Þ
with x and y being the spatial coordinates (Fig. 1). The bottom left hand corner of the block is constrained in all directions while the nodes along the bottom edge are free to move in the x-direction only. The model is subjected to a compressive engineering strain (eeng) of 0.5 in the ydirection. The engineering strain is defined as the ratio of the increase in height to the original height: (h1 h0)/h0, where h1 is the compressed height and h0 is the original height. A rate sensitivity exponent of n = 50 is used. This gives a moderate rate-dependent response and is a physically representative value for the material being modelled [23].
h0 ¼ 380 MPa
1
Several values of n are used. A quantity of importance in the above equation is the accumulated slip, ca. This is a measure of the total crystallographic plastic strain at a point and is defined as follows: X Z t c_ ðaÞ dt ca ¼ ð2:3Þ a
0
where t is the time and the summation is over all the slip systems at a point. The material chosen, 316L stainless steel, has a face-centred-cubic crystalline structure. Slip is assumed to occur along specific slip directions as defined by its lattice and the slip system {1 1 1}h1 1 0i. The above constitutive theory is implemented in the finite element code ABAQUS/Standard as a user-defined subroutine (UMAT), as coded in [10]. This subroutine, written in Fortran, implements the theory in the form of a stress update algorithm that is called at each integration point in the mesh and for every iteration during a finite element simulation [9].
Fig. 1. Single grain model—40 · 80 mesh. Units are in mm.
Fig. 2. Cube depicting slip direction and crystal orientation for the single crystal compression test, a 54.7. All slip directions are parallel to the three shown.
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
445
crystal axes, with respect to the global axes, are specified as material properties in the user interface. A crystal orientation is used whereby the [1 0 1] and the ½1 0 1 crystallographic directions are parallel to the global x and z-axes, respectively. In this arrangement, in a two-dimensional analysis, several slip directions become coincident (Fig. 2). Thus, twelve possible slip directions are reduced to three; two oblique and the horizontal direction. The angle subtended by the oblique and the horizontal slip directions in 2D is a 54.7. In this arrangement the directions of slip are assured to be oblique in the model and the formation of slip bands is promoted. Fig. 3. Force–displacement curves for 10 · 20, 20 · 40, 40 · 80 and 100 · 200 meshes for the rate sensitivity parameter, n = 50. Point A indicates a strain of eeng = 0.38 which is the upper limit for generation of mesh independent results.
The model is meshed with eight-noded generalised plane strain quadrilateral elements (CPEG8R), each having four (2 · 2) integration points. Four increasingly fine mesh densities are used; models with 10 · 20, 20 · 40, 40 · 80, and 100 · 200 elements are generated. The orientations of the
3.2. Results The force applied to the top edge of the model against the vertical displacement of the edge for all four meshes is graphed in Fig. 3. The four curves compare well at lower levels of compression. One would be tempted to say that both convergence and a result independent of element size had been achieved. However the two denser meshes show a softening that is due to the formation of a sharp shear band across the thickness of the material. This softening can be considered to represent failure of the material in this case.
Fig. 4. Contour plots of the accumulated shear strain in the 40 · 80 mesh at increasing levels of compression: (a) eeng = 0.18, (b) 0.32, (c) 0.5.
446
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
It is accompanied by considerable mesh dependence where the shear band thickness is determined by the element size. The contour plots in this paper show the accumulated slip, ca (Eq. (2.3)). This value is a measure of the crystallographic plastic strain and is analogous, but not equal, to the equivalent plastic strain. Fig. 4(a)–(c) shows the 40 · 80 element mesh of the single grain at increasing levels of strain. At eeng = 0.3 (Fig. 4(b)) the localised areas of shearing form bands at regular intervals along the length
of the grain, running at oblique angles to the direction of straining. The two shear bands have similar widths and profiles. As the strain level increases the shear bands become narrower and the strain gradients become steeper, with the band on the right hand side becoming particularly localised (Fig. 4(c)). A visual comparison of the four meshes at eeng = 0.5 shows that, as the shear band becomes narrow, there is a high level of mesh sensitivity (Fig. 5(a)–(d)). The thickness of the shear band, in terms
Fig. 5. Contour plots (with mesh removed for clarity) of the accumulated shear strain at eeng = 0.5: (a) 10 · 20 mesh, (b) 20 · 40 mesh, (c) 40 · 80 mesh, (d) 100 · 200 mesh, (e) 120 · 240 mesh.
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
Fig. 6. Close-up of the shear band in the 40 · 80 mesh at eeng = 0.5. The highlighted elements indicate that the width of the shear band is set by the span of one element.
of element size, is difficult to assess as the mesh is highly distorted within the band. It can be said that there is only one element completely within the shear band and it spans the width of the band, as illustrated in Fig. 6. This is the case with the two denser meshes that exhibit material failure (Fig. 3). It is interesting to note the formation of shear bands at the right hand edge of the model. The 100 · 200 mesh shows two strong shear bands and the initiation of a third one at the top corner (Fig. 5(d)). All three bands are caused by the compression of the material along the concave part of the geometrical imperfection. The gradients of strain at the same location in the 40 · 80 mesh are not as sharp (Fig. 5(c)). It is noticeable that the increase in mesh density merely improves the resolution of the solution, it does not increase the number of shear bands present. Whereas, the 20 · 40 mesh has only two shear bands along the right hand edge (Fig. 5(b)). It would appear that the multiple shear bands are part of the numerical solution to the problem and that the number of shear bands is independent of
447
mesh density above a certain level, i.e. the initiation of shear bands on the right hand edge does not appear to be a function of the number of elements used. It is estimated that a mesh density higher than the 100 · 200 case would not yield an increase in the number of shear bands but would cause narrower shear bands. In order to illustrate this point, a more dense mesh with 120 · 240 elements was constructed. The resulting contour plot (Fig. 5(e)) shows that a higher mesh density does not yield a higher number of shear bands than the 100 · 200 case. It is reasonable to assume that a limit, in terms of shear banding, has been reached. 3.3. Rate sensitivity To assess the influence of the strain rate sensitivity parameter, n, on shear band formation, simulations are performed using the above models where n is varied. More rate-dependent values of n = 10 and 20 and the quite rateindependent value, n = 100, are used. The three cases are modelled using the 20 · 40 and 40 · 80 meshes. It is known that material rate sensitivity delays the development of shear bands [14,16]. This is shown where the increasing strain rate sensitivity blunts the shear strain gradient. The resulting deformed mesh for n = 20 shows a much more diffuse shear band running northwest to southeast (Fig. 7). Unlike in the original case, the denser mesh does not produce a narrow shear band whose width is dependent on the element size. There is a better agreement between the resulting contours in the two models. Although, the profiles of the localised strains are not converged, the macroscopic results are not sensitive to these differences and compare favourably (Fig. 8). When the rate sensitivity is increased further, with n = 10, the strain localisations become even more diffuse (Fig. 9). As a result the denser mesh does not identify
Fig. 7. Contour plots of the accumulated shear strain for rate sensitivity value n = 20 at eeng = 0.5: (a) 20 · 40 mesh, (b) 40 · 80 mesh.
448
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
Fig. 8. Force–displacement curves for 20 · 40 and 40 · 80 meshes for the rate sensitivity parameter, n = 20.
any gradients that are not present in the sparser one. The contour plots and the force–displacement curves (Fig. 10) compare well for the two mesh densities and can be said to have converged. When a relatively rate-independent value is used, n = 100, the formation of narrow shear bands is quite evident. As would be expected, neither the applied force– displacement responses (Fig. 11) nor the contour plots (Fig. 12) compare well. This is an indication of the mesh sensitivity of results that is inherent in the rate-independent finite element analyses. The above simulations show that as well as being deformation dependent, the mesh sensitivity of a model is also heavily reliant on the degree of rate sensitivity of the material.
Fig. 10. Force–displacement curves for 20 · 40 and 40 · 80 meshes for the rate sensitivity parameter, n = 10.
Fig. 11. Force–displacement curves for 20 · 40 and 40 · 80 meshes for the rate sensitivity parameter, n = 100.
3.4. Discussion The finite element analysis of shear band formation in a rate-dependent material is shown to be heavily reliant on two factors. The first of these factors is the level of rate sen-
sitivity of the material and the second is the amount of deformation of the model. Increasing the rate sensitivity has the effect of delaying the onset and sharpening of the
Fig. 9. Contour plots of the accumulated shear strain for rate sensitivity value n = 10 at eeng = 0.5: (a) 20 · 40 mesh, (b) 40 · 80 mesh.
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
449
Fig. 12. Contour plots of the accumulated shear strain for rate sensitivity value n = 100 at eeng = 0.5: (a) 20 · 40 mesh, (b) 40 · 80 mesh.
strain gradient. Thus the model does not require the high level of resolution that a dense mesh provides. The sharp strain gradient present in the denser meshes and the more rate insensitive materials is a highly mesh sensitive state. The width of the shear band is such that one element spans it, regardless of mesh density. This is caused by excessive element distortion and is accompanied by the loss of load-bearing capability in the model. This can be viewed as failure of the material. Although several values of the rate sensitivity parameter were used for the same type of analysis, this is not a possibility for practical simulations when an accurate material description must be used. Therefore, the mesh density must be adjusted in order to yield mesh independent results. For the n = 20 compression case presented in Section 3.3 there is a good level of agreement in terms of the contour plots (Fig. 7) and the macroscopic mechanical responses (Fig. 8). Although the profiles of ca for the two meshes are not identical this is not necessary to get a converged, correct solution for the macroscopic behaviour. A similar situation occurs in the n = 50 case (Section 3.2) where there is a good level of agreement up to a strain level of eeng = 0.38, as indicated by point A in Fig. 3. We can infer from this that a solution can be generated that is effectively independent of the element size used so long as the
overall deformation of the model falls below a certain strain limit (eeng < 0.38, in this case). 4. Polycrystal tension 4.1. Model description The stress–strain behaviour of small scale 316L stainless steel devices is dependent on their size. This size dependence affects, in particular, the point of failure [13]. Polycrystalline finite element models incorporating a crystal plasticity constitutive subroutine have been shown to capture this size effect successfully [19]. In order to investigate the effect of mesh refinement in larger models a 2D polycrystalline model of a 316L stainless steel wire of 75 lm thickness is constructed. The model is based on that presented in Savage et al. [19]. The granular microstructure is explicitly represented with each grain idealised as a hexagon of area 92 lm2 (Fig. 13). As outlined in Section 2, the crystal plasticity constitutive model used is implemented via a user material subroutine (UMAT) in ABAQUS. Although the model is based on a 2D geometry each grain has a full 3D crystallographic description that is randomly oriented for each grain. In addition, the use of 8-noded generalised plane strain elements (CPEG8R) allows for
Fig. 13. Model of 316L stainless steel wire with a thickness of 75 lm. The granular microstructure is idealised as a series of hexagons (highlighted).
450
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
Fig. 14. Macroscopic stress–strain curves for all the meshes analysed in orientation set 2. The legend refers to the number of elements per grain in each mesh. Inset: stress–strain curves showing an agreement in the failure strains for the three denser meshes.
uniform out-of-plane deformation. Therefore the model is considered ‘‘quasi-3D’’. A rate sensitivity parameter of n = 50 is used throughout. To achieve a reasonable spread of data, three sets of randomly assigned orientations are used. Several different mesh resolutions are analysed. Mesh densities of 8, 50, 98, 162 and 242 elements per full grain are constructed. The left hand bottom corner of the model is constrained in all degrees of freedom and the nodes along the left hand edge are free to move in the y-direction only. The wire is tensed by applying a uniform horizontal (x-direction) displacement to the right hand edge of the model. The engineering strain is defined as the ratio of the increase in length to the original length: (l1 l0)/l0, where l1 is the stretched length and l0 is the original length. The engineering stress is the ratio of force applied to the wire to the original cross-sectional area of the wire: F/A0. This produces engineering stress–engineering strain (reng–eeng) curves and is a convenient method for presenting tensile data. In these simulations the failure strain is taken as the strain at the point of ultimate tensile stress (UTS). 4.2. Results Each grain displays a different hardness in the direction of tension due to the random crystal orientation. The softer grains act as initiation sites for strain localisations. This results in non-uniform deformation and coalescence of areas of plastic strain, eventually forming a shear band that
traverses the width of the wire. Although an explicit material separation model, representing fracture, is not included in the current crystal plasticity formulation, the shear band indicates the site of necking and eventual failure of the wire. It is important to note that pockets of localised shearing occur at many places throughout the wire but eventual necking occurs at the weakest part, i.e. where the grain lattice orientations are most favourable for plastic slip. Following the failure strain point (UTS), localisations continue to grow and coalesce leading to macroscopic softening. The results presented are for the second set of random orientations assigned to the grains and are typical of all three sets. The macroscopic reng–eeng curves, comparing the response of the different mesh densities, are shown in Fig. 14. The three more densely meshed geometries (98, 162 and 242 elements per grain) show good agreement and have converged in terms of failure strain (Fig. 14). The hardness of the wire reduces as the number of elements in the model increases. This is due to the increase in number of degrees of freedom in the model. It is also important to investigate the contours of the localisations present in the granular microstructure of the material. Fig. 15 shows the profile of ca at eeng = 0.6 from one of the lattice orientation sets. The shear localisations lie in bands at oblique angles to the tensile axis. To illustrate a convergence in the results a single grain that lies within the shear band is monitored throughout the analysis for the different meshes. The contours of ca in a single grain
Fig. 15. Accumulated shear strain during necking of a 75 lm wire. Shear localisations form at angles oblique to the tensile axis.
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
451
Fig. 16. Accumulated shear strain, ca, in a grain at the point of failure, eeng = 0.35. The grain lies within the shear band that traverses the wire. The contour limits of each plot are the same.
at the point of failure are shown in Fig. 16 for all five mesh densities. The profiles of the shear band in the three denser meshes are quite similar. As the straining continues the gradients become steeper and narrower. The same grain is examined at a strain of eeng = 0.44; Fig. 17(a)–(e) shows ca for all five meshes. The shear band profiles do not appear to agree as well as in Fig. 16. In this case an increase in number of elements continues to yield an improvement in the results and further mesh refinement would be required to achieve convergence. As each grain undergoes higher deformation the strain gradients become more severe and a smaller element size is needed. This is similar to the single grain compression simulation. In contrast to the single grain compression case the point of failure, as identified in the macroscopic response, does not correspond to development of a very narrow, element size dependent shear band. 4.3. Discussion Whereas in Section 3 analyses are carried out that introduce shear bands through a geometrical imperfection, in
this section the cause of the strain localisations is the material inhomogeneity. Due to the microstructural arrangement in the wire, each grain experiences non-uniform deformation under large straining. Each grain can be viewed as an isolated block of material that is subjected to non-uniform boundary deformation. Unlike the single crystal compression case the point of macroscopic failure does not correspond to the initiation of a narrow shear band. Here, the softening of the material is due to the strains in the softer grains becoming very large in comparison to the harder grains and leading to a reduction in the cross-sectional area of the wire. The favourable comparison of the reng–eeng responses in the three more densely meshed grains up to the point of UTS (Fig. 14) is a similar case to that seen in Fig. 3. Although the contour plots of ca are not the same (Fig. 16), the macroscopic responses have converged in terms of the failure point (Fig. 14). In the polycrystal tension simulation the wire is not tensed beyond a strain of eeng = 0.6 (Fig. 15). Were the wire to be stretched further one would see the formation of a shear band with a thickness of one element as is seen in the single grain analyses. Prior to this situation it is possible
452
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
Fig. 17. Accumulated shear strain, ca, in a grain after the point of failure, eeng = 0.42. The grain lies within the shear band that traverses the wire. The contour limits of each plot are the same.
to get a solution that is accurate and unaffected by the element size used. In the current analyses a mesh density of 98 elements per grain has been shown to provide a sufficient level of accuracy in order to provide the failure point. This translates to a mesh to grain area ratio of 0.01. Although the in-grain contour profiles of ca are not completely converged there is still a good agreement between this case and the subsequent more finely meshed models (Fig. 16(c)–(e)). The more persuasive argument for choosing this mesh spacing is by examining the macroscopic performances (Fig. 14). The model with the greater number of elements continues to yield an overall reduction in the stiffness of the wire. This is due to the increased number of degrees of freedom in the mesh. Whereas the strain at which the UTS occurs does not change for the three more finely meshed models. 5. Conclusion This paper shows that computational models that are capable of producing results independent of the element size can be developed. In single crystal compression analyses, excessive straining causes a situation where shear
bands show considerable mesh dependence. Prior to this situation a converged solution was found in both macroscopic and in-grain mechanical performance. The study was then extended to include the tensile behaviour of a 75 lm thickness wire. The microstructure of the model was explicitly represented and the material properties were chosen to approximate 316L stainless steel behaviour. It was found that the point of failure, taken as the strain at the UTS, did not correspond to a situation of critical mesh dependence. A mesh density per grain could be identified that gave a relatively accurate profile of the intragranular plastic strains and a converged value for the macroscopic failure strain. The accurate prediction of failure is vital in the design of small scale metallic components such as cardiovascular stents in the medical device industry or components in the microelectronics industry. The results of the work show that, for the application considered here, by employing an element area to grain area ratio of 0.01 there is an acceptable ‘window’ of macroscopic strain that can be applied to the model while maintaining accuracy of solution. This window of strain extends to the point of failure. This mesh density can be used for further analysis to investigate the
F.J. Harewood, P.E. McHugh / Computational Materials Science 37 (2006) 442–453
mechanical performance of small-scale metallic devices under various loading conditions, with confidence in the level of accuracy provided by the results. Acknowledgements The authors wish to acknowledge funding from Embark, Irish Research Council for Science, Engineering and Technology: Funded by the National Development Plan. The simulations in this work were performed on the SGI 3800 high performance computer at NUI, Galway. References [1] L. Anand, M. Kothari, Journal of the Mechanics and Physics of Solids 44 (4) (1996) 525–558. [2] D.J. Bammann, Materials Science and Engineering A 309–310 (2001) 406–410. [3] T. Belytschko, Z.P. Bazant, Yul-Woong Hyun, Ta-Peng Chang, Computers and Structures 23 (2) (1986) 163–180. [4] T.E. Buchheit, G.W. Wellman, C.C. Battaile, International Journal of Plasticity 21 (2) (2005) 221–249. [5] R. De Borst, L.J. Sluys, H.-B. Muhlhaus, J. Pamin, Engineering Computations 10 (1993) 99–121. [6] F. Delaire, J.L. Raphanel, C. Rey, Acta Materialia 48 (5) (2000) 1075–1087. [7] K. Garikipati, T.J.R. Hughes, Computer Methods in Applied Mechanics and Engineering 159 (3–4) (1998) 193–222. [8] S.V. Harren, H.E. Deve, R.J. Asaro, Acta Metallurgica 36 (9) (1988) 2435–2480.
453
[9] Hibbitt, Karlsson, Sorenson, ABAQUS Theory Manual, Pawtucket, RI, USA, 1997. [10] Y. Huang, A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program, Harvard University Report, MECH 178, 1991. [11] J.P. McGarry, B.P. O’Donnell, P.E. McHugh, J.G. McGarry, Computational Materials Science 31 (3–4) (2004) 421–438. [12] P.E. McHugh, P.J. Connolly, Computational Materials Science 27 (4) (2003) 423–436. [13] B.P. Murphy, P. Savage, P.E. McHugh, D.F. Quinn, Annals of Biomedical Engineering 31 (6) (2003) 686–691. [14] A. Needleman, Computer Methods in Applied Mechanics and Engineering 67 (1) (1988) 69–85. [15] A. Needleman, V. Tvergaard, Element Analysis of Localization in Plasticity, in: Finite Elements—Special problems in Solid Mechanics, Prentice-Hall, New York, 1983, pp. 94–157. [16] D. Peirce, R. Asaro, A. Needleman, Acta Metallurgica et Materialia 31 (12) (1983) 1951–1976. [17] D. Raabe, Z. Zhao, S.-J. Park, F. Roters, Acta Materialia 50 (2) (2002) 421–440. [18] M. Sachtleber, Z. Zhao, D. Raabe, Materials Science and Engineering A 336 (1–2) (2002) 81–87. [19] P. Savage, B.P. O’Donnell, P.E. McHugh, B.P. Murphy, D.F. Quinn, Annals of Biomedical Engineering 32 (2) (2004) 202–211. [20] B.M. Schroeter, D.L. McDowell, International Journal of Plasticity 19 (9) (2003) 1355–1376. [21] A. Staroselsky, L. Anand, Journal of the Mechanics and Physics of Solids 46 (4) (1998) 671–696. [22] A. Tatschl, O. Kolednik, Materials Science and Engineering A 356 (1– 2) (2003) 447–463. [23] X. You, in preparation. [24] N. Zhang, W. Tong, International Journal of Plasticity 20 (3) (2004) 523–542.