CC BY
96
УДК 620.172.224.2
Micromechanical modeling approach to derive the yield surface for bcc and fcc steels using statistically informed microstructure models and nonlocal crystal plasticity
N. Vajragupta, S. Ahmed, M. Boeff, A. Ma, and A. Hartmaier
Interdisciplinary Centre for Advanced Materials Simulation, Ruhr-Universität Bochum, Bochum, 44801, Germany
In order to describe irreversible deformation during metal forming processes, the yield surface is one of the most important criteria. Because of their simplicity and efficiency, analytical yield functions along with experimental guidelines for parameterization become increasingly important for engineering applications. However, the relationship between most of these models and microstructural features are still limited. Hence, we propose to use micromechanical modeling, which considers important microstructural features, as a part of the solution to this missing link. This study aims at the development of a micromechanical modeling strategy to calibrate material parameters for the advanced analytical initial yield function Barlat YLD 2004-18p. To accomplish this, the representative volume element is firstly created based on a method making use of the statistical description of microstructure morphology as input parameter. Such method couples particle simulations to radical Voronoi tessellations to generate realistic virtual microstructures as representative volume elements. Afterwards, a nonlocal crystal plasticity model is applied to describe the plastic deformation of the representative volume element by crystal plasticity finite element simulation. Subsequently, an algorithm to construct the yield surface based on the crystal plasticity finite element simulation is developed. The primary objectives of this proposed algorithm are to automatically capture and extract the yield loci under various loading conditions. Finally, a nonlinear least square optimization is applied to determine the material parameters of Barlat YLD 2004-18p initial yield function of representative volume element, mimicking generic properties of BCC and FCC steels from the numerical simulations.
Keywords: representative volume element, virtual microstructures, weighted Voronoi tessellation, nonlocal crystal plasticity finite element method, BarlatYLD 2004-18p yield function
Микромеханическое моделирование поверхности текучести для сталей с ОЦК- и ГЦК-решеткой на основе статистически обоснованных моделей микроструктуры и нелокальной модели пластичности кристаллов
N. Vajragupta, S. Ahmed, M. Boeff, A. Ma, A. Hartmaier
Рурский университет, Бохум, 44801, Германия
Поверхность текучести является одним из наиболее важных критериев при описании необратимой деформации в процессе металлообработки. Благодаря своей простоте и эффективности аналитические функции текучести наряду с экспериментальными способами параметризации приобретают все большее значение для инженерных приложений. Однако большинство указанных моделей не позволяют в полной мере учитывать микроструктурные особенности. Для решения данной проблемы предлагается использовать микромеханическое моделирование, учитывающее важные особенности микроструктуры. В статье предложено развитие подхода к микромеханическому моделированию для калибровки параметров материала для расширенной аналитической функции начала текучести Barlat YLD 2004-18p. На основе метода, использующего статистическое описание морфологии микроструктуры в качестве входного параметра, смоделирован представительный элемент объема. Данный метод сочетает в себе моделирование с помощью частиц и разбиение Вороного для создания реалистичных виртуальных микроструктур, которые принимаются за представительные элементы объема. Использована нелокальная модель пластичности кристаллов для описания пластической деформации представительного элемента объема методом конечно-элементного моделирования в рамках теории пластичности кристаллов. Таким образом, предложен алгоритм построения поверхности текучести на основе конечно-элементной модели в рамках теории пластичности кристаллов. Основным назначением предлагаемого алгоритма является автоматическое определение и выделение площадок текучести в разных условиях нагружения. Проведена оптимизация с помощью нелинейного метода наименьших квадратов и определены параметры материала функции начала текучести Barlat YLD 2004-18p представительного элемента объема, имитирующие общие свойства сталей с ОЦК- и ГЦК-решеткой путем численного моделирования.
Ключевые слова: представительный элемент объема, виртуальные микроструктуры, взвешенная диаграмма Вороного, метод конечных элементов нелокальной теории пластичности кристаллов, функция текучести BarlatYLD 2004-18p
© Vajragupta N., Ahmed S., Boeff M., Ma A., Hartmaier A., 2017
1. Introduction
With respect to the development tools for sheet metal forming processes, irreversible or plastic deformation process must be predicted with a high reliability. For engineering applications, analytical yield functions are usually used because of their high numerical efficiency and straightforward application [1]. Since the first introduction of the von Mises yield function [2], several analytical yield functions that consider anisotropy, resulting for example from textures, have been established to describe plasticity of various materials. For example, the quadratic anisotropic yield function proposed by Hill [3] is suitable for body centered cubic (bcc) metals such as ferritic steels [1], but this function is not adequate for the case of face centered cubic (fcc) metals like aluminum alloys [4]. To overcome such limitation, nonquadratic anisotropic yield functions have been proposed [5, 6]. Furthermore, linear transformations of the stress tensor are used to further enhance the yield functions [7-10].
With respect to engineering applications, material models along with a limited set of simple and low-cost tests, e.g. uniaxial tensile tests, required for parameters calibration are preferable [11]. However, as these advanced analytical yield functions consist of several material parameters, rather complex and time consuming experimental procedures are necessary to fully parameterize the models. Furthermore, the relationship between such advanced yield functions and microstructural features are still unknown and the material parameters do not have a direct physical or microstructural interpretation. From the perspective of materials science, it is mandatory to establish such a relationship between microstructure and mechanical properties to guide the design of new microstructures of materials to meet properties required in their applications.
In consequence, the micromechanical modeling approach, which explicitly takes into account important microstructural features such as orientation distribution functions, grain size and shape distribution functions, is considered here as a solution to provide this missing link. One of the ultimate aims of this modeling approach is to predict macroscopic mechanical behavior from simulations of microstructure models by means of homogenization techniques. Thus, macroscopic stress and strain tensors at the onset of plasticity can be obtained for various loading conditions and used as data for fitting parameters of advanced yield functions. In this way, the finite element simulation of microstructure can be considered as a virtual mechanical lab. With this technique, the high costs of experimental setups used to parameterize macroscopic model can potentially be significantly reduced.
The aim of this study is to develop a microstructure-based simulation strategy to determine parameters for an advanced initial yield function for bcc and fcc steels. This micromechanical modeling strategy involves several steps
including creation of a microstructure model, which quantitatively describes the real microstructure, implementation of adequate constitutive equations capturing the deformation behavior of polycrystalline aggregates, formulation of homogenization techniques to obtain macroscopic stress and strain response from micromechanical simulations, and application of optimization scheme to parameterize the yield function. In this work, we choose the Barlat YLD 2004-18p model [8], because of its flexibility to describe the anisotropic yield behavior.
At the center of micromechanical modeling, the method used to generate the representative volume element (RVE) is crucial as it must be able to capture microstructural features as grain size and grain shape distributions quantitatively. Even though there are several approaches for the construction of RVE described in the literature, the concept of the Voronoi tessellation is still one of most favorable techniques, because it requires only a small number of input parameters. In order to better control grain size using this tessellation technique, several authors make use of the weighted Voronoi tessellation which assigns additional weight property to the seeds to generate RVE [12-14]. The weighted Voronoi tessellation techniques make use of the position along with corresponding weight to subdivide the space. This weight is defined to manipulate the size of segmented volume or grain with respect to the context of microstructure generation. Generally, there are various mathematical expressions describing geometry of the bisector for the weighted Voronoi tessellation [15]. Among these formulations, the radical Voronoi tessellation is one of the straightforward weighted Voronoi tessellations since the bisector is a straight line for 2-dimensional case and a plane for 3-dimensional case. Therefore, the radical Voronoi tessellation has been widely used because of its simplicity for generating geometries [16].
Within the framework of micromechanical modeling, crystal plasticity models that consider the influence of crys-tallographic orientation on deformation for both single crystals and polycrystalline systems have been intensively developed [17-19]. One of the crucial advantages of crystal plasticity models is their capability to solve complex crystal mechanics under complicated internal and/or external boundary conditions, as it allows one to address the issue of inter- and intragrain micromechanical interaction [20]. This leads to an improved understanding on the abrupt mechanical transitions along the interface. Nevertheless, most of the constitutive models in crystal plasticity theory neglect the influence of deformation gradients. However, in some applications, these size effects are important, as observed in experiments, including bending of polycrys-talline nickel [21], microbending experiments of single crystal aluminum and single crystal copper [22, 23], and twisting of polycrystalline copper [24]. To tackle the influence of deformation gradient, advanced nonlocal constitutive
models have been proposed. Most of these nonlocal formulations are derived based on the concept of the geometrically necessary dislocation (GND) density tensor as introduced by Nye [25]. Therefore, the nonlocal constitutive model should be coupled with RVE in order to precisely capture the deformation gradient especially along the grain boundaries.
2. RVE construction using dynamic microstructure generator
For this study, quasi-2D RVEs are generated using the dynamic microstructure generator (DMG) method [16]. This procedure couples particle simulations with the radical Voronoi tessellation to create the targeted microstructure. In general, the radical Voronoi tessellation makes use of seed positions along with corresponding diameters of spheres as input parameters. To obtain grain with the size, which is comparable to the sphere, the sphere distribution should be closely packed [16].
There are several methods to distribute spheres into space, including the random sequential addition (RSA) [26]. However, implementing RSA usually results in nonclose packed distribution of spheres in the microstructural space, which leads to some deviations from the targeted grain size distribution. Hence, a particle simulation method is implemented using LAMMPS [27, 28] to optimize the volume of the grains. With this simulation technique, size and a number of sphere are predefined in the way that a given grain size distribution can be mimicked [12]. In the first step, spheres are randomly distributed into a finite volume, which is larger than the intended final RVE. In the second step, the volume is compressed while the spheres move freely under a repulsive potential to avoid their overlapping. In the third step, the updated sphere positions and diameter assigned to each sphere are then transferred to the radical Voronoi tessellation package from open-source software Voro++ [29] for selected time steps. The resulting grain size distributions are then compared to the predefined size distribution of spheres and the RVE generated with
minimum difference is selected. It must be noted that the shape of the RVE created using DMG is rugged to leave grain intact, which improves the mesh quality and retains the periodicity [16]. Fitting the grains into periodic boxes with straight edges, require the cutting of the grains and typically results in partitions that are difficult to mesh.
For this study, the target grain size distribution is defined via a log-normal distribution with average grain diameter of 5 ^m and a standard deviation of 1 ^m. According to the prescribed distribution parameters, RVE of 100 grains is created as shown in Fig. 1, a. Comparison between discretized size distribution of seed spheres and resulting grain size distribution can also be found in Fig. 1, b. For this feasibility study, we have chosen a two-dimensional (2D) RVE; but we note that the generation of three-dimensional (3D) ones is possible with the same procedure.
To finish the RVE for micromechanical simulations, the geometry of 2D RVE is extruded for 0.2 ^m and meshed with 8 nodes linear brick elements (C3D8) using CUBIT [30]. This RVE consists of 8437 elements. It must be noted that the generated RVE from DMG is periodic and in order to apply the periodic boundary conditions, opposite surfaces must have exactly the same surface mesh, which requires a well-controlled meshing process. Finally, a random orientation distribution is assigned to the grains of the RVE.
3. Nonlocal crystal plasticity finite element method
Within the scope of this work, a crystal plasticity model is applied to describe plastic deformation within the individual grains of a polycrystal. Note that each grain is dis-cretized into several finite elements such that each element represents a part of a single crystalline region with a specific crystal orientation. The applied crystal plasticity method has already been described in detail in the literature [31], hence, here only an overview of the formulation is provided to define the used parameters and to make this work reproducible.
With respect to the kinematics of finite deformation, the deformation gradient F can be multiplicatively decom-
Seed
EIIU Grain
Seeds/grain diameter, Lim
Fig. 1. Generated RVE using DMG method; meshed RVE (a) and comparison of diameter distribution between defined seed spheres and constructed grain (b)
posed into the elastic deformation gradient Fe and plastic deformation gradient Fp:
F = Fe F p. (1)
The elastic deformations are calculated in the usual way, giving rise to stresses according to Hooke's law. The aniso-tropic elastic constants used in this work are provided in Table 1. The plastic deformation is characterized by the plastic velocity gradient Ip, as a function of the plastic deformation gradient Fp and its rate as
Lp = Fp Fp-1. (2)
In this work, we only consider plastic deformation by crystallographic slip of dislocations, hence II is described as sum of the plastic shear rate of all slip systems, yielding
(3)
ns
LL = £Y a Ma,
a=1
where ya is the plastic shear rate and M%a = da ® na is the Schmid tensor for slip system a, which is defined by the slip direction da and the slip plane normal na. The number of slip systems considered is Ns. The symbol ® denotes the dyadic product of two vectors resulting in a second rank tensor. For small elastic strains, the resolved shear stress Ta for each slip system can be calculated from the stress S in the intermediate configuration as
Ta= S • Ma. (4)
Following Ma and Hartmaier [31], the flow rule and strain hardening law are expressed as
GNDfc p
Y a=Y 0
T™ + T™
T a+T
GNDi
a
sign(Tc
. GND/K , + Ta )'
Ns
T„ = £ ^oXapl1 j 1 Y|
P 1
(5)
(6)
Here, yo is the reference shear rate, p1 is the inverse value of the strain rate sensitivity, h0 is the initial hardening rate, and X«p is the cross hardening matrix, which is assigned as 1.0 for coplanar slip systems and 1.4 otherwise, Tf is the saturation slip resistance, and p2 is a fitting parameter. The initial value for the initial value of critical resolved shear stress T0 is provided in Table 1. The flow rule in Eq. (5) contains furthermore two back stresses TGNDi and tGNDi that describe the additional hardening caused by GNDs resulting from strain gradients [31]. The nonlocal constitutive model applied here is formulated based on the concept of densities of geometrically necessary superdislocations. This formulation allows us to incorporate the pronounced plastic strain gradient at grain boundaries occurring in our RVEs. According to this model, the strain gradients at grain boundaries are considered as pile-ups of GNDs, which result in additional strain hardening. This additional strain hardening is split up into an isotropic hardening contribution t gnd and a kinematic hardening con-
tribution T
.gndi
For the sake of completeness, the basic
provided in the following, for details the reader is kindly referred to [31].
The second rank dislocation density tensor G in the reference configuration is calculated from the curl of Fp as defined by Nye [25]
Gy =-( Fp xV)r (7)
Since the dislocation density tensor only contains information about averages of dislocations and residual Burgers vectors, it is not possible to reconstruct meaningful crystallographic dislocation populations in a unique way. However, by projecting the dislocation density tensor to the global Cartesian coordinates of the system, a unique definition of geometrically necessary superdislocations is achieved and the far field stress of the crystallographic GND population can be described with good accuracy [31]. In this way, the GND density tensor can be separated into nine independent parts pa by evaluating
£ Padx® ta= 1G, (8)
a=1 __b
where da and ta are permutations of the Cartesian unit vectors as specified in [31] and b is the norm of the crystallographic Burgers vector. Please note that the super GND densities for a = 1, ..., 3 represent screw-type superdislocations, whereas the remaining components represent edge-type superdislocations, which is important when evaluating the internal stress fields resulting from the super GNDs.
The first contribution of these super GNDs to the work hardening is simply expressed by a Taylor-type work hardening formulation in the form
iGNDi
= c^l £ xGND I Pp I,
(9)
P=1
Table 1
Summary of material parameters for BCC and FCC steels
equations of the formulation of this nonlocal model are
BCC FCC
Phenomenological C11, MPa 231000 242000
C12, MPa 134700 146500
C44, MPa 116400 122000
t0, MPa 40 150
y 0 0.001
P1 45.5 200.0
h0, MPa 500 345
Tf, MPa 400 480
P2 5.00 2.25
GND L, nm 1
c1 0.01
xgnd /lap 1.00
Non-Schmid «1 «2 «3 «4 a5 a6
0.61 0.23 0.55 0.11 0.09 -0.20
where c1 is the Taylor hardening coefficient or a geometrical factor [32] and ^ = C44 is the shear modulus, X^® is the cross hardening matrix between crystallographic mobile dislocations and geometrically necessary superdislocations (super GNDs) [31]. This contribution causes an additional isotropic hardening for dislocation slip caused by plastic strain gradients.
A second contribution to work hardening results from the far-ranging internal stresses caused by GNDs in dislocation pile-ups. This contribution is taken into account by evaluating the second order gradient of Fp, which results in a super GND gradient paI in the form [31]
Pa,I = IbGjk ,i da® Ta. (10)
By evaluating these gradients within small volumes of dimension I3, it is possible to calculate internal stresses SGND in the intermediate configuration caused by dislocation pile-ups at grain boundaries, as described in [31]. From these internal stresses, finally, the second hardening contribution is quantified as the back stress tgnd k % GND
Ta =S ' Ma (11)
giving rise to additional kinematic work hardening.
For fcc steel, within the scope of this study, the above-mentioned constitutive laws are directly implemented and the applied material parameters are summarized from [31] and shown in Table 1. We considers dislocation slip on the usual crystallographic (111}(110) slip systems.
Furthermore, these equations can also be applied for bcc steels if the corresponding bcc slip systems are considered. In this work we take into account that the slip resistance for the glide of screw dislocations with nonplanar dislocation core in bcc systems does not only depend on the resolved shear stress, but also on further components of the stress tensor, see for example [33]. Here, we describe the influence of these so-called non-Schmid stresses in the form developed by Koster et al. [34], where the influence of the other shear components as well as the normal stress components on screw dislocation mobility are
-BCC
fourth and fifth terms consider the influence of tension and compression perpendicular to the glide direction while the sixth term is added to exclude the influence of hydrostatic stresses. It is noted that all parameters at have been calculated by atomistic methods in [33] and are applied in this work.
For this study, the material parameters of bcc steel are taken from [34] which are also given in Table 1. As bcc slip systems, only the {110}(111) systems are taken into account, which results in Ns = 12 slip systems.
4. Virtual mechanical testing
The nonlocal crystal plasticity model described above is implemented onto a user defined material subroutine (UMAT) [35] and applied in finite element simulations with the commercial software ABAQUS [36] to assess the mechanical properties of RVEs in so-called virtual mechanical tests. In these simulations, periodic boundary conditions are applied to fully maintain the periodicity during deformation [12]. The key principle of the periodic boundary condition is that the degrees of freedom of nodes on one side of the model are coupled to those of their counterparts on the opposite side. However, care must be taken to allow lateral shrinkage due to cross-contraction effects. The implemented formulations of periodic boundary condition applied in this work are based on the work of Smit et al. [37]. Details of the implementation are described in [16].
In order to calibrate the initial yield function, the mechanical response of RVE is tested under various multiaxial loading conditions, as shown in Fig. 2. In order to quantify the macroscopic loading condition of RVE, the mechanical boundary conditions of the finite element model must be transferred into macroscopic stress and strain tensors. The implementation of periodic boundary conditions, requires that the deformation state is applied to the four reference nodes Vx, V2, V4, and H1 as marked in Fig. 3. The RVE boundary nodes follow the kinetics of these reference nodes. Thus, the macroscopic strain tensor can eas-
node
taken into account as a further back stress Ta within the ily be calculated from the nodal displacements uf of
flow rule in the form
:Y o
_ - _GND& Ta + Ta
T +T GND
-T
BCC
sign( T
■ _GNDk \
(12)
these nodes [16]. The mathematical expression of macroscopic strain tensor £RVE is written as
This back stress is expressed as
BCC
: aiG : da ® naS + a2 a :(na x da ) ® na +
+ a3a ! (na,i X da
)® na ,1 + a4 a :na® na +
+ a5 a: (na X da)® (nax da) + a6a :da ® da. (13) The first term of the right-hand-side of Eq. (13) takes into account the influence of twinning-antitwinning asymmetry. The second and third terms integrate the effect of shear stresses perpendicular to the slip direction, na 1 is the specific {110} plane normal vector, which includes an angle of -60° with the reference slip plane defined by na. The
Fig. 2. All the combinations of applied strain in directions x andy
Fig. 3. Schematic diagram illustrating the vertex nodes Vj, V2, V4, and H
£rve -
uj2 Ax
/ v4 V2 Ul + u2l
Ay Ax
v
^ Hj V
uj 1
V4 V2 Ul + U22
Ay Ax
V4
u24
Ay
H
Hj
uj 1
+
V2
u32
Ax
) uH^+uV4'
Az Ay
V
y
H
.(14)
+ U3 j u2[ + U3 u3
Az Ax 2 Az Ay Az
V / V /
where Ax, Ay, and Az are the dimensions of the periodic
box in the global Cartesian coordinate system. The macroscopic stress tensor can be formulated from the reaction force vectors Fnode at the four reference nodes and the current nodal position vectors xnode of the reference nodes
which is written as [16] j
a = -
V
-sym[( xV
V
) ®FV +
RVE
+ (xV2 -xVj)®FV2 + (xH
j) ® F
H
(15)
where the symmetrization function is defined as sym = = j/2[A + AT] for tensor A and its transposed AT. It is noted that all RVE models have a thickness of only one element, such that the simulations are considered to be quasi-2D in nature under effective plane stress conditions.
5. Yield function determination and micromechanical RVE analysis
From the virtual mechanical testing of RVEs, macroscopic stress and strain tensors at different loading conditions are obtained. Here, we use the stresses at the onset of plastic yielding to parameterize the yield function Barlat YLD 2004-18p [8]. This yield function is selected, because it shows a great flexibility in describing anisotropic yield behavior. Thus, we can demonstrate that such advanced yield functions are an effective way to homogenize the results of micromechanical models. This chapter aims at briefly introducing this yield function along with the least square optimization method to determine the parameters.
5.1. Barlat YLD 2004-18p
Barlat YLD 2004-18p [8] is one of the most complex yield functions available, which is derived based on two
linear transformations of the deviatoric stress tensor. This yield function can be expressed as
«ms'-sy + \s1'-s2'\n + \s1'-£J'r +
+ \ s2' - s'\n + \ s2' - S2'\n + \ s2' - S3'|n +
+ \ S3' - s' \ nn + \ S3' - ¿2" r + \ ¿3' - S3' \ nn = 4an, (16)
where ay is the equivalent stress, n is the exponent which is either 6 or 8 for the case of bcc and fcc material respectively [37], St and St (i = 1, 2, 3) are the principal values of tensors s' and s" determining from two linear transformations of the deviatoric part a of the Cauchy stress tensor:
s = C: a, (17)
s" = D: a, (18)
C and D are the fourth-order tensors which are mainly responsible for characterizing the anisotropy of material or shape of the yield function. Each of this tensor contains of 9 material parameters which can be expressed with Voigt notation as
C =
D =
0
- C21
-C31 0 0 0 0
- d 21 -d31 0 0 0
- cl2 0
-c32 0 0 0
- dl2 0
-d32 0 0 0
23
- cl3
- c 0 0 0 0 -d
j3
- d
23
0 0 0
c44 0 0
0 0 0
d 44 0 0
0 0 0 0
c55 0
0 0 0 0 0
C66
(19)
0 0 0 0
d55 0
0 0 0 0 0
d66
(20)
and d66 define the out-of-
The parameters c55, c66, d55 plane properties while other 14 parameters describe the inplane properties of the sheet. If all parameters are set to a value of 1 and the exponent n in Eq. (16) is an even number, the Barlat yield function corresponds to the von Mises yield function.
5.2. Least square optimization
To calibrate 18 anisotropy parameters from simulations of RVE within the scope of this study, the nonlinear least square method for multivariable problems is used [38]. The error function E to be minimized is the sum of the quadratic relative difference between the yield stress aYLD obtained by evaluating the Barlat yield function and the yield stress aRVE obtained from the RVE simulation as E (C, D) -J , (a YLD/ a RVE)2. (21)
The minimization of this error function is done by varying all 18 components of the tensors C and D with the trust region reflective algorithm [32]. Convergence is assumed, when the error remains constant during several iterations.
Table 2
Summary of identified parameters for Barlat YLD 2004-18p of bcc and fcc steel
bcc fcc
c12 1.013 d12 0.774 c12 0.9990 d12 0.437
c13 1.091 d13 0.973 c13 0.9660 d13 0.941
c21 1.011 d21 0.903 c21 0.9010 d21 1.184
c23 1.020 d23 1.026 c23 1.0750 d23 0.954
c31 0.756 d31 0.981 c31 0.9990 d31 1.000
c32 0.894 d32 1.124 c32 1.0000 d32 0.999
c44 0.00001 d44 -1.510 c44 -0.00001 d44 0.371
c55 1.000 d55 1.000 c55 -1.5710 d55 0.999
c66 1.000 d66 0.999 c66 0.8240 d66 1.000
The RVE yield stress grve is defined as the macroscopic equivalent stress (calculated from the stress tensor defined by Eq. (15)) at the macroscopic equivalent plastic strain of 0.2 % (according to Eq. (14)).
5.3. Fitted initial yield function
The residual values of the error function for bcc and fcc steels reach the minimum values of 0.0014 and 0.0013, respectively. Thus it can be stated that the Barlat YLD 2004-18p, due to its 18 material parameters, allows a sufficient degree of flexibility for fitting even the complex shapes of initial yield surfaces obtained from RVE simulations with good precision. The identified parameters are summarized in Table 2. It can be seen that several parameters stay close to the value of 1, which indicates that they do not contribute to distortions of the yield surface with respect to the isotropic case. It is noted, here that in this study the uniqueness of the results has not been investigated, because the feasibility of the description of the numerical yield points stemming from RVE simulations with an analytical yield function was in the focus. Now that this feasibility has been confirmed, future work will address uniqueness and robustness of the parameter identification.
In Fig. 4, finally, the fitted initial yield functions together with the homogenized macroscopic yield stress are
plotted for both bcc and fcc steels. It can be seen that both initial yield functions fit the data of the RVE simulations very well and moreover provide a meaningful interpolation between the discrete data points. Thus, it is concluded that the mechanical behavior of RVE with respect to the onset of plastic yielding can successfully be described and homogenized with the Barlat yield function.
5.4. Micromechanical analysis of RVE
To illustrate the effect of applying nonlocal crystal plasticity for RVE simulations, the internal stresses and strains at different loading stages for uniaxial tension are further analyzed. The evolution of the von Mises stress distribution in RVE is shown in Fig. 5. Noticeably, the von Mises stress is always highest in the grain boundary region. This can be explained by the nonlocal crystal plasticity model, taking dislocation pile-ups at grain boundaries and their contribution to the work hardening behavior into account. As a consequence, the internal strains as shown in Fig. 6 are distributed more homogeneously within the grain, because further plastic deformation of the grain boundary regions is reduced by work hardening. Such behavior stands in contrast to simulations with local crystal plasticity models, where typically the plastic strains are localized in the grain boundary regions, due to the strong deformation incompatibilities there, whereas the stresses are rather homogeneous within one grain.
6. Conclusions
In this study, we demonstrate the feasibility of describing the plastic yield behavior of a micromechanical material model with the help of macroscopic yield functions. The micromechanical model is defined as a representative volume element that possesses a polycrystalline microstructure with the realistic grain size distribution. The mechanical behavior of RVE is described with the nonlocal crystal plasticity model that takes deformation gradients and the resulting dislocation pile-ups at grain boundaries into account. This crystal plasticity model has been parameterized for two cases, describing the generic behavior of single phase bcc and fcc steels, respectively.
Gyy, MPa.
- Fitted Yld2004-18p • Virtual experiment
-100-
-200 -100
100 <jxx, MPa
Gyy, MPa
- Fitted Yld2004-18p \b_ • Virtual experiment
300 g^, MPa
Fig. 4. Initial yield functions according to the Barlat Yld2004-18p yield function and numerical data obtain from RVE simulations for bcc (a) and fcc steel (b)
a V, MPa
1 600
450
z
300
1 150
aV, MPa
E 600
450
■ ..300
■ 150
*m
*
'M
«11: :■■■: : ■ .. i e
SiB
m
■
Pill
0
600 450
300 150
i
600 450 300 150
I
< EK
Fig. 5. Evolution of von Mises stress aV distribution during uniaxial tension simulation of RVE at various macroscopic equivalent plastic strains: 1.35 (a), 2.72 (b), 4.01 (c) and 5.41 % (d)
The mechanical behavior of this micromechanical model is characterized by finite element simulations under various uniaxial and multiaxial external loads. The onset of plastic yielding is defined by an engineering approach as the stress at which an equivalent plastic strain of 0.2 % is reached. In our work, we demonstrate that these yield stresses can be described very accurately by the Barlat YLD
2004-18p [8] yield function. Thus, a new way to homogenize micromechanical models is demonstrated that, furthermore, also offers the possibility to link the parameters of the macroscopic yield function to microstructural parameters or to specific features of crystal plasticity models.
In future work, the proposed micromechanical modeling strategy shall be further extended for fitting the evolu-
PEEQ 0.10
PEEQ
S 0.10
0.07 0.04 0.01
aV, MPa
Fig. 6. Evolution of equivalent plastic strain PEEQ distribution during uniaxial tension simulation of RVE at various macroscopic equivalent plastic strains: 1.35 (a), 2.72 (b), 4.01 (c) and 5.41 % (d)
tion of the yield surface during plastic flow, i.e. to describe also the micromechanical work hardening behavior correctly by a macroscopic work hardening law. Furthermore, the uniqueness and robustness of this homogenization method will be studied.
References
1. Zhang H., Diehl M., Roters F., Raabe D. A virtual laboratory using high resolution crystal plasticity simulations to determine the initial yield surface for sheet metal forming operations // Int. J. Plasticity. -2016. - V. 80. - P. 111-138.
2. von Mises R. Mechanik der plastischen Formanderung von Kristallen // Z. Angew. Math. Mech. - 1928. - V. 8. - P. 161-185.
3. Hill R. A theory of the yielding and plastic flow of anisotropic metals // Proc. R. Soc. A. Math. Phys. Eng. Sci. - 1948. - V. 193. - P. 281297.
4. Hosford W.F. Limitations of Non-Quadratic Anisotropic Yield Criteria
and Their Use in Analysis of Sheet Forming // Proc. of 15th Congress of Int. Deep Drawing Research Group. - Materials Park: ASM Int., 1988. - P. 163-170.
5. HillR. Theoretical plasticity of textured aggregates // Math. Proc. Cam-
bridge Philos. Soc. - 1979. - V. 85. - No. 1. - P. 179-191.
6. Hosford W.F. Comments on anisotropic yield criteria // Int. J. Mech. Sci. - 1985. - V. 27. - No. 7-8. - P. 423-427.
7. Barlat F., Lian K. Plastic behavior and stretchability of sheet metals. Part I: A yield function for orthotropic sheets under plane stress conditions // Int. J. Plasticity. - 1989. - V. 5. - No. 1. - P. 51-66.
8. Barlat F., Aretz H., Yoon J.W, Karabin M.E., Brem J.C., Dick R.E. Linear transformation-based anisotropic yield functions // Int. J. Plasticity. - 2005. - V. 21. - No. 5. - P. 1009-1039.
9. Karafillis A.P., Boyce M.C. A general anisotropic anisotropic yield criterion using bounds and a transformation weighting tensor // J. Mech. Phys. Solids. - 1993. - V. 41. - No. 12. - P. 1859-1886.
10. Banabic D., Kuwabara T., Balan T., Comsa D.S., Julean D. Non-quadratic yield criterion for orthotropic sheet metals under plane stress conditions // Int. J. Mech. Sci. - 2003. - V. 45. - No. 5. -P. 797-811.
11. Lademo O.G., Hopperstad O.S., Langseth M. An evaluation of yield criteria and flow rules for aluminum alloys // Int. J. Plasticity. - 1999. -V. 15. - No. 2. - P. 191-208.
12. Vajragupta N., Wechsuwanmanee P., Lian J., Sharaf M., Munstermann S., Ma A., Hartmaier A., Bleck W. The modelling scheme to evaluate the influence of microstructure features on microcrack formation in DP-steel: The artificial microstructure model and its applications to predict the strain hardening behaviour // Comp. Mater. Sci. -2014. - V. 94. - P. 198-213.
13. Regener B., Krempaszky C., Werner E., Stockinger M. Modelling the micromorphology of heat treated Ti6Al4V forgings by means of spatial tessellations feasible for FEM analyses of microscale residual stresses // Comp. Mater. Sci. - 2012. - V. 52. - P. 77-81.
14. Wu Y., Zhou W., Wang B., Yang F. Modeling and characterization of two-phase composites by Voronoi diagram in the Laguerre geometry based on random close packing of spheres // Comp. Mater. Sci. -2010. - V. 47. - No. 4. - P. 951-961.
15. Okabe A., Boots B., Sugihara K., Chiu S.N. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. - Chichester: John Wiley & Sons Ltd., 2000.
16. Boeff M. Micromechanical Modelling of Fatigue Crack Initiation and Growth: PhD Thesis. - Bochum: Ruhr-Universitat Bochum, 2016.
17. Roters F. Overview of constitutive laws kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications // Acta Mater. - 2010. - V. 58. -P. 1152-1211.
18. Tikhovsky I., Raabe D., Roters F. Simulation of earing of a 17% Cr stainless steel considering texture gradients // Mat. Sci. Eng. A. Struct. - 2008. - V. 488. - P. 482-490.
19. Wang Y Orientation dependence of nanoindentation pile-up patterns and of nanoindentation microtextures in copper single crystals // Acta Mater. - 2004. - V. 52. - P. 2229-2238.
20. Sachtleber M., Zhao Z., Raabe D. Experimental investigation ofplastic grain interaction // Mat. Sci. Eng. A. - 2002. - V. 336. - No. 1-2. -P. 81-87.
21. Stoelken J., Evans A. A microbend test method for measuring the plasticity length scale // Acta Mater. - 1998. - V. 46. - P. 5109-5115.
22. Suzuki K., Matsuki Y., Masaki K., Sato M., Kuroda M. Tensile and microbend tests of pure aluminum foils with different thicknesses // Mat. Sci. Eng. A. - 2009. - V. 513. - P. 77-82.
23. Hayashi I., Sato M., Kuroda M. Stain hardening in bent copper foils // J. Mech. Phys. Solids. - 2011. - V. 59. - P. 1731-1751.
24. FleckN.A., Muller G.M., Ashby M.F., Hutchinson J.W. Strain gradient plasticity: Theory and experiment // Acta Metall. Mater. - 1994. -V. 42. - P. 475-487.
25. Nye J.F. Some geometrical relations in dislocated crystals // Acta Metall. - 1953. - V. 1. - P. 153-162.
26. Tarjus G., Schaaf P., Talbot J. Random sequential addition: A distribution function approach // JSP. - 1991. - V. 63. - P. 167-202.
27. Plimpton S. Fast parallel algorithms for short-range molecular dynamics // J. Comput. Phys. - 1995. - V. 117. - No. 1. - P. 1-19.
28. Plimpton S. LAMMPS Molecular Dynamics Simulator. - http:// lammps.sandia.gov. Accessed: 2016-02-01.
29. Rycroft C.H., Grest G.S., Landry J.W., Bazant M.Z. Analysis of granular flow in a pebble-bed nuclear reactor // Phys. Rev. E. - 2006. -V. 74. - P. 021306.
30. Sandia National Laboratories "CUBIT 13.2.". - https://cubit. sandia.gov/. Accessed: 2016-05-04.
31. Ma A., Hartmaier A. On the influence of isotropic and kinematic hardening caused by strain gradients on the deformation behaviour of polycrystals // Philos. Mag. - 2014. - V. 94. - No. 2. - P. 125-140.
32. Newville M, Ingargiola A., Stensitzki T., Allen D.B. LMFT: NonLinear Least-Square Minimization and Curve-Fitting for Python. -Zenodo, 2014.
33. Vitek V., Mrovec M., Bassani J.L. Influence of non-glide stresses on plastic flow: From atomistic to continuum modeling // Mater. Sci. Eng. A. - 2004. - V. 365. - No. 1. - P. 31-37.
34. Koester A., Ma A., Hartmaier A. Atomistically informed crystal plasticity model for body-centered cubic iron // Acta Mater. - 2012. -V. 60. - P. 3894-3901.
35. ABAQUS User Subroutine Reference Manual Version 6.7. Dassault Systems: 2007.
36. ABAQUS/Analysis User's Manual. Version 6.12. ABAQUS Inc.
37. Smit R., Brekelmans W., Meijer H. Prediction of the mechanical behavior of nonlinear heterogeneous systems by multi-level finite element modeling // Comput. Method. Appl. Mech. Eng. - 1998. -V. 155.- No. 1-2. - P. 181-192.
38. Hosford W.F. A generalized isotropic yield criterion // Int. J. Appl. Mech. - 1972. - V. 39. - No. 2. - P. 607-609.
Поступила в редакцию 01.02.2017 г.
Сведения об авторах
Napat Vajragupta, Dr.-Ing., Res. Group Leader, ICAMS, Ruhr-Universität Bochum, Germany, [email protected] Shabaz Ahmed, Stud. Ass., ICAMS, Ruhr-Universität Bochum, [email protected] Martin Boeff, Dr.-Ing., ICAMS, Ruhr-Universität Bochum, [email protected]
Alexander Hartmaier, Prof., Dr., Director, ICAMS, Ruhr-Universität Bochum, [email protected]
