FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 17, N o 2, 2019, pp. 243 - 254 https://doi.org/10.22190/FUME190403028R © 2019 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper MICROSTRUCTURE-BASED SIMULATIONS OF QUASISTATIC DEFORMATION USING AN EXPLICIT DYNAMIC APPROACH Varvara Romanova 1 , Ruslan Balokhonov 1 , Evgeniya Emelianova 1 , Olga Zinovieva 2 , Aleksandr Zinoviev 2 1 Institute of Strength Physics and Materials Science, SB RAS, Tomsk, Russia 2 Airbus Endowed Chair for Integrative Simulation and Engineering of Materials and Processes, University of Bremen, Bremen, Germany Abstract. Microstructure-based simulations of the deformation processes require substantial computational resources due to the necessity of using detailed meshes with a large number of elements. An approach that considerably reduces the computational costs implies simulation of quasistatic deformation within a dynamic approach involving a solution of the motion equations rather than the equilibrium equations. It enables a transition from implicit to explicit time integration providing a significant gain in the computational capacity. In this paper, we show that the explicit dynamic approach can be successfully used in the microstructure-based simulations of quasistatic deformation, considerably reducing the computational costs without losing the information and solution accuracy. The following conditions have to be met to ensure a close agreement between the dynamic and static solutions: (i) the load velocity in the dynamic calculations must be smoothly increased to its amplitude value and then kept constant to minimize the acceleration term appearing in the equation of motion and (ii) the constitutive model employed must describe a quasi-rate- independent response. An examination of the mesh convergence and the strain-rate dependence for a polycrystalline aluminum model has supported this conclusion. Key Words: Microstructure-based Simulations, Quasistatic Deformation, Explicit Dynamic Approach, Crystal Plasticity Received April 03, 2019 / Accepted June 28, 2019 Corresponding author: Varvara Romanova Institute of Strength Physics and Materials Science SB RAS Akademicheskii prospect 2/4, 634055 Tomsk, Russia E-mail: varvara@ispms.tsc.ru 244 V. ROMANOVA, R. BALOKHONOV, E. EMELIANOVA, O. ZINOVIEVA, A. ZINOVIEV 1. INTRODUCTION The microscale deformation mechanisms are known to be strongly affected by the material microstructure. The knowledge of deformation and fracture mechanisms operating in loaded materials at different scales is of critical importance since gradual accumulation of irreversible microdeformation and damage is commonly followed by macroscopic failure of the engineering structure. Along with the experimental methods, numerical simulations explicitly taking the material microstructure into account appear to be useful tools for studying the multiscale deformation processes. While considerable progress in this field has been made in the recent few decades, the microstructure-based numerical analysis in a 3D case remains to be a challenge for the researchers due to mathematical complexity of the 3D problem, difficulties in its numerical implementation and high computational demands on numerical solving the boundary-value problem (BVP). On the one hand, the microstructure model has to contain a sufficient number of structural elements for the micro- and mesoscale processes to be simulated as realistically as possible; the microstructure constituents and interface regions have to be approximated in sufficient detail to ensure a reasonable accuracy of the solution. This necessitates the use of detailed meshes with a large number of elements. In some cases, the high-resolution meshes require the memory and computational time so large that the microstructure-based numerical analyses become impractical. It is, therefore, challenging to reduce the computational costs without losing the information and solution accuracy. An approach that considerably reduces the memory, disk space and computational time requirements implies simulation of quasistatic deformation processes in a dynamic formulation where the equations of motion are solved instead of the equations of equilibrium [1-3]. This enables a transition from implicit to explicit time integration, which provides significant advantages from the viewpoint of computational capacity. The benefit of the dynamic calculations becomes much more significant for any kind of nonlinearity, e.g., nonlinear constitutive behavior, microstructural inhomogeneity, nonlinear loading history, and the like. Among the numerical problems, namely, those where the explicit dynamic approach has a distinct advantage over the static one, there are contact problems in which it is quite difficult to make the implicit solution converge [4]. The drawback of the dynamic simulations is the conditional stability of the numerical scheme, which places strong restrictions on the time step value so that too many computational steps would be necessary to achieve a reasonable degree of straining at quasistatic loading rates. To overcome this trouble the loading is artificially sped up, which under certain conditions would result in wave dynamics untypical for quasistatic deformation. Moreover, the material free surface and interfaces of different kinds (e.g., grain boundaries, matrix–particle or substrate–coating interfaces, etc.) are the sources of wave reflection, refraction and dissipation and as such they would affect the numerical solution to a greater or lesser extent. For this reason, the dynamic approach has limited applications in the microstructure-based simulations, where the material interfaces are treated explicitly. In this paper we discuss the dynamic approach applicability in microstructure-based simulations of quasistatic deformation phenomena. The conditions under which the static and dynamic solutions converge to a high degree of accuracy are analyzed for an aluminum polycrystal as an example. Microstructure-based Simulations of Quasistatic Deformation Using an Explicit Dynamic Approach 245 2. BOUNDARY-VALUE PROBLEMS IN STATICS AND DYNAMICS The mechanical boundary-value problems (BVPs) discussed in this Section are formulated in a rectangular Cartesian coordinate system in the absence of body forces. In the dynamic formulation the elastic-plastic BVP includes the equations of motion ,i ij j u  , (1) the strain rate-displacement relations , , 1 ( ) 2 ij i j j i u u   (2) and the constitutive equations in the rate form of the generalized Hooke’s law ( ) e p ij ijkl kl ijkl kl kl C C      . (3) Here ρ is the current density, ui are the components of the displacement vector, σij are the stress tensor components, Cijkl is the fourth-order tensor of elastic moduli, εij, ε e ij and ε p ij are the components of the total, elastic and plastic strain tensors; the upper dot denotes the time derivative. The kinematic boundary conditions take the form i iS u   , (4) and the traction boundary conditions are ij j i S n T    , (5) where υi are the velocity vector components prescribed on the Sυ surface, and Ti and ni are the force and normal vector components on the Sσ surface. The governing equations of a static problem are the equations of equilibrium , 0 ij j   (6) complemented by the strain-displacement relations , , 1 ( ) 2 ij i j j i u u   . (7) and the constitutive equations (3). The traction boundary conditions are given by Eq. (5) while the kinematic boundary conditions are given by surface displacements Ui prescribed on the SU surface U i iS u U . (8) Among the numerical methods for solving partial differential equations (PDEs) the finite-element (FE) method is considered to be the most universal and widely used one. It implies a transition from the strong formulation of the BVPs to a weak integral form using, e.g., the virtual work principle. The integral equations are approximated by a system of algebraic equations determined on a finite-element mesh. The finite-element 246 V. ROMANOVA, R. BALOKHONOV, E. EMELIANOVA, O. ZINOVIEVA, A. ZINOVIEV method is discussed in detail in a large number of papers (e.g., [5]). Let us briefly address the main features of the FE implementation using explicit and implicit approaches. Both in dynamic and static formulations the deformation process is simulated in a stepwise manner: the load applied to the boundaries is incremented and the stress, strain and displacement fields are updated at the end of each time increment to achieve a new state of the static or dynamic balance. The time integration methods in current use are based on implicit or explicit approaches. With an implicit FE solver the unknown quantities to be calculated at each time increment are expressed through the parameters that are also unknown at the beginning of this increment. Therefore, iteration algorithms are necessary to achieve a numerical solution. The static equilibrium Eq. (6) can be solved by implicit numerical methods alone. The general form of the global FE equations of the static BVP is    [ ] 0 K u f , (9) where {f} and {u} are the global vectors of the nodal forces and displacements, and [K] is the global stiffness matrix relating the forces and displacement vectors. In the implicit computational procedure the stiffness matrix has to be inverted, which requires substantial computational resources. While the implicit solution is assumed to be attained at relatively large time increments, the iteration procedure and calculations of the inverse stiffness matrix make the computations rather expensive in terms of memory, disk space, and computational time. Any kind of nonlinearity (e.g., nonlinear constitutive behavior, irregularly-shaped interfaces, complex geometry of the computational domain, nonlinear loading path, etc.) additionally increases the number of iterations necessary to distinguish the features of nonlinear phenomena to a proper accuracy. In contrast to the static equilibrium problem, hyperbolic equations of motion can be solved using an explicit scheme. The unknown quantities appearing in the PDEs are expressed through the parameters known from the previous time step. The components of acceleration and velocity vectors at the (n+1)-th time increment are expressed through the values known from the n-th and (n+1/2)-th time steps 1 1 1/ 2n n n n i i i u u t u      , (10) 1 1/ 2 1/ 2 2 n n n n n i i i t t u u u         . (11) The accelerations are calculated at the beginning of the time increment from the equation of motion written in the matrix form as follows       [ ]  M u f K u , (12) where [M] is the lumped mass matrix. The explicit solution requires neither iterations nor inversion of the stiffness matrix, which substantially reduces the computational costs for each time increment. The drawback of the explicit schemes is their conditional stability. The stability condition is associated with the velocity of the fastest process described by the PDEs. For the mechanical problem the time step has to be proportional to the smallest element size and inversely proportional to the velocity of mechanical wave propagation in the material (i.e., the sound velocity) with the coefficient <1. The stability condition Microstructure-based Simulations of Quasistatic Deformation Using an Explicit Dynamic Approach 247 ensures that the stress wave across each time increment does not cover the distance exceeding the smallest mesh step. While each increment of the explicit time integration is much less computationally consuming than that of the implicit calculations, the time step providing the stability of the explicit scheme is too small for the simulations of long-time processes to be practical. Particularly, too many computational steps would be necessary to achieve a reasonable degree of straining at quasistatic loading rates. To overcome these troubles, the load velocities in the explicit simulations of quasistatic processes are artificially increased. Another method for reducing the computational time in the dynamic calculations involves scaling the material density to speed up the stress wave propagation. 3. MICROSTRUCTURE-BASED SIMULATIONS 3.1 Generation of polycrystalline geometrical models The construction of a geometrical microstructure model, which is the starting point in the microstructure-based numerical analysis, is a sophisticated problem in a 3D case. A rigorous approach is based on processing a series of experimental microstructural images obtained by means of specimen sectioning, X-ray tomography and other time and money consuming techniques. Alternative methods rely on the computer-aided design of synthetic models with microstructural features similar to those of real materials. Earlier [6] we have proposed a semi-analytical method of step-by-step packing (SSP), which enables generating 3D microstructures with a wide variety of geometrical features. In the general case, the SPP-procedure includes the following steps. A 3D volume is discretized with a regular or irregular mesh. Certain mesh elements are selected to be seeds of the microstructure elements. Each kind of seeds is associated with a certain analytical function according to which the volumes surrounding the seeds are grown in a stepwise manner. The mesh elements, whose coordinates fall within any of the incremented seed volumes, are added to this microstructure phase. The SSP-procedure is repeated until the growing phases reach the preset volume content. The main parameters controlling the resulting microstructure geometry are the number of seeds, their types and spatial distributions, and the laws of their growth. The equations of ellipsoids, spheres, cylinders and the like are the basic analytical functions enabling us to construct microstructures typical of many materials. In this paper the SSP-method has been employed to generate three-dimensional polycrystalline models on regular meshes with different resolutions. In order to obtain the same polycrystalline structure on different meshes a set of grain seeds was once randomly selected and then applied in all SSP-simulations. All grains were grown by the equation of a sphere at the same growth rate. The SSP generation was terminated when the entire computational domain was packed with grains. This algorithm provides polycrystalline aggregates with convex polyhedral grains characterized by plane faces and straight-line edges much similar to those constructed analytically by a Voronoi tessellation. An advantage of the SSP-models generated on a mesh by default is that they can be directly imported into the finite-element or finite-difference computations. 248 V. ROMANOVA, R. BALOKHONOV, E. EMELIANOVA, O. ZINOVIEVA, A. ZINOVIEV Two polycrystalline models generated on the coarsest and finest meshes are shown in Fig. 1 together with a selected grain presented for 100×20×100, 160×32×160, 200×40×200, 250×50×250, 320×64×320 and 400×80×400 meshes. The mesh resolution was found to have but a minor effect on the generated polycrystalline structures. While finer meshes describe smoother grain interfaces, the shape, size and spatial arrangement of the grains are the same as those for coarser meshes. a) b) Fig. 1 Polycrystalline models generated on 100×20×100 (a) and 400×80×400 meshes (b) and a selected grain meshed with a step of 20, 12.5, 10, 6.25 and 5 µm (left-to-right) 3.2 Constitutive behavior of aluminum grains The constitutive response of aluminum grains is described in the framework of the crystal plasticity theory [5]. A polycrystal is treated as an aggregate of single crystals of varying crystallographic orientations with respect to the specimen coordinate system (XYZ-axes in Fig. 1). The generalized Hooke’s law (3) is written with regard to a local coordinate system associated with the crystal axes. The plastic strains appearing in Eq. (3) are thought to result from dislocation gliding in the active slip systems ( ) ( )p ij ij       , with ( ) ( )1 ( ) 2 ij i j j i s m s m      (13) where Si (α) and mi (α) are the components of slip direction and slip plane normal vectors for a slip system (SS) α. Equations (13) provide a relation between specimen strains and dislocation slip in the directions prescribed by the orientation tensor Θij. The shear strain rate ( )  in Eq. (13) is the unknown quantity and has to be determined by a constitutive dependence on the resolved shear stress τ (α) =σijΘij (α) acting on the slip system α. Rate- independent models provide an ambiguous definition for a set of active slip systems [7]. In order to overcome this ambiguity, the rate-dependent models are commonly used in crystal plasticity simulations Microstructure-based Simulations of Quasistatic Deformation Using an Explicit Dynamic Approach 249 ( ) ( ) * ( ) ( ) sgn CRSS            , (14) where *  and ν are the parameters controlling the strain rate sensitivity. Of primary importance is the description of critical resolved shear stress (CRSS) τ (α) CRSS with proper account of the strengthening mechanisms for particular materials. In this paper, a simple phenomenological equation for τ (α) CRSS is used to describe the grain boundary strengthening and strain hardening ( ) ( ) 0 ( ) poly p CRSS eq f         , (15) where τ0 (α) is the CRSS of a single crystal, τ poly takes into account the CRSS value increasing due to the presence of grain boundaries. The third term of the sum describes the strain hardening as a function of the accumulated equivalent plastic strain ε p eq. In what follows, calculations are presented for an aluminum alloy which is characterized by face-centered cubic (FCC) crystalline structure. Twelve <111>{110} SSs are potentially active in FCC metals, all of them are activated at the same value of τ (α) CRSS. The strain hardening function for the aluminum alloy is determined as 2 2 / / 1 1 ( ) (1 ) (1 ) p p eq eqa bp eq f a e b e          , (16) where a1, a2, b1 and b2 are chosen to fit the experimental data [8]. The model parameters and material constants used in the simulations are C1111=108 GPa, C1122=61 GPa, C2323=28 GPa, τ0 (α) =2 MPa, τ poly =18 MPa, a1=73 MPa, a2=0.07, b1=16 MPa, b2=0.002. 3.3 Numerical implementation The polycrystalline constitutive models were imported into the finite-element software package ABAQUS/Standard and ABAQUS/Explicit through UMAT and VUMAT User Subroutines, respectively. The crystal plasticity FE-implementation in ABAQUS/Standard is reported in [5]. Let us consider briefly the explicit numerical procedure. The simultaneous solution to Eqs. (3), (13) and (14) at each time increment calls for an iterative procedure. We employed the method of simple iterations which provided a fast solution convergence with a reasonable accuracy. The solution is shown to converge for one or two iterations provided that the parameters of Eq. (14) are well-defined, with a purely elastic state being chosen as the initial approximation. In the ABAQUS/Explicit solver, the constitutive equations are formulated with respect to local orientations given by Θij. The tensors and vectors used in the calculations of constitutive behavior are automatically rotated with respect to the local coordinates before they are imported to VUMAT. Thus, we do not have to reformulate the constitutive Eqs. (3), (13)-(17) for each local coordinate system. The boundary conditions, formulated with respect to the specimen coordinate system, simulate uniaxial tension along the X-axis. The specimen top surface in all simulations is free of external forces, while the bottom surface is taken as a symmetry plane. The lateral surfaces parallel to the tensile axis are assumed to be free of external forces. 250 V. ROMANOVA, R. BALOKHONOV, E. EMELIANOVA, O. ZINOVIEVA, A. ZINOVIEV 4. COMPUTATIONAL RESULTS 4.1 Mesh convergence The mesh convergence of the numerical solution has been checked for six polycrystalline models of 0.2×0.04×0.2 cm consisting of 1600 grains generated on different meshes. Two models generated on the coarsest and finest meshes are shown in Fig. 1. The set of grain orientations was identical in all simulations. a) b) c) d) e) f) Fig. 2 Equivalent stress (a-c) and plastic strain fields (d-f) in the polycrystalline structure meshed with 100×20×100 (a, d), 200×40×200 (b, e) and 400×80×400 resolution (c, f) Explicit calculations of uniaxial tension have been performed for the models loaded at the same strain rate of 10 2 s -1 . All components of the stress and plastic strain tensors obtained for different meshes have been compared to analyze the effect of the mesh resolution on the numerical solution accuracy. For the sake of illustration, the equivalent stress and plastic strain fields are presented in Fig. 2 for three meshes, with the general conclusion being supported by the whole set of numerical data. The respective stress-strain curves and those of mesh dependence of the maximum stress and strain values are plotted in Fig. 3. For all mesh approximations the stress and strain distributions are found to be reproduced with a reasonable accuracy. Even the stress and strain fields calculated for the coarsest mesh (Fig. 2a, d) qualitatively reproduce the main features which become more detailed on the finer meshes (Fig. 2b-c, e-f). For the meshes finer than 160×32×160, qualitative differences between the stress-strain patterns become hardly distinguishable. Microstructure-based Simulations of Quasistatic Deformation Using an Explicit Dynamic Approach 251 100 200 300 400 340 360 380 400  eq Fitting curve m a x . e q u iv a le n t s tr e s s , M P a No. of mesh elements along X-axis 0,08 0,10 0,12 0,14  p eq m a x . e q u iv a le n t p la s tic s tra in 0,00 0,03 0,06 0,09 0 100 200 300 S tr e s s , M P a Strain Mesh size 100x20x100 160x32x160 200x40x200 320x64x320 400x80x400 × a) b) Fig. 3 Mesh dependence of the maximum values of equivalent stresses and plastic strains (ε = 0.029) (a) and averaged stress-strain curves for different mesh densities (b) The stress-strain curves for all mesh approximations mostly coincide except the curve for the coarsest mesh which lies somewhat lower (Fig. 3b). It is worth noting that the maximum values of stresses and strains also tend to converge upon mesh refinement (Fig. 3a), though at a slower rate than the averaged values do. This supports the conclusion made by Harewood and McHugh [1] that the rate-dependent models enable eliminating mesh sensitivity of the solution when plastic strain localizes in shear bands. 4.2 Explicit dynamic simulations of quasi-static deformation In quasi-static simulations using the dynamic approach it is critical to ensure that the inertia effects are insignificant. Apparently, the dynamic and static solutions converge if the inertia term appearing in the left-hand side of Eq. (1) vanishes. The acceleration value is non-zero in the cases where the velocity is changed and is neglected when it is kept constant. Thus, the load velocity in the initial deformation stage has to be increased smoothly to minimize the acceleration and, thus, to eliminate the dynamic effects involved. We have found that the wave effects become negligible if the time of the velocity increase up to the amplitude value is 3-4 times larger than that necessary for the elastic wave to propagate through the computational domain. Taking into account the ratio of the model longest length to the speed of elastic wave propagation in aluminum, the time of load increasing in the explicit simulations was chosen to be 3 µs. The solutions to the implicit static and explicit dynamic problems have been compared for the polycrystalline structure approximated by a mesh consisting of 1,600,000 elements. Note that the maximum number of elements, which might be calculated with ABAQUS/Explicit, is an order of magnitude larger than that in the implicit static simulations run in the same computer. The element-by-element comparison between static and dynamic stress and strain fields showed a coincidence to within 0.1% for the most part of the elements, with only few elements belonging to interface regions demonstrating a disagreement of 3-5% (Fig. 4a). 252 V. ROMANOVA, R. BALOKHONOV, E. EMELIANOVA, O. ZINOVIEVA, A. ZINOVIEV a) 0,0 0,5 3,5 4,0 0 20 40 60 N o f fi n it e e le m e n ts , % Discordance between dynamic and static solutions, % b) 0,00 0,04 0,08 0,12 0 100 200 300 Strain S tr e s s , M P a Strain rate, 1/s 10 2 10 3 10 4 Fig. 4 The element-by-element comparison between plastic strain fields obtained in dynamic and static calculations (a) and the stress-strain curves calculated using an explicit dynamic approach at different load velocities (b) In the dynamic simulations of quasi-static deformation, where load velocities are artificially increased, it is important to use rate-independent constitutive models. Thus, the strain-rate sensitivity parameters, appearing in the rate-dependent crystal plasticity Eq. (14), have to be chosen in order to eliminate the rate dependence. Harewood and McHugh [1] reported that a material becomes quasi-rate-independent for large ν values, which may have an adverse effect on the iteration convergence. In our simulations the ν value was equal to 10, which reduced the strain-rate sensitivity of the average material response for the strain rates up to 10 3 s -1 . The grain scale stress and strain fields, however, demonstrated conspicuous differences. With increasing the load velocity, the regions of plastic strain localization became wider, while the strain values in the shear bands decreased (cf., e.g. Fig. 5a and b). This is due to the fact that the local strain rates in the strain localization regions can be several orders of magnitude higher than the strain rate prescribed by the boundary conditions. If so, the plastic strain rates calculated by Eqs. (13)-(14) might not be high enough to achieve a static balance between the load and the material response. To overcome this drawback, we suggest substituting constant *  in Eq. (14) by the relationship * eq k  , (17) where eq  is the equivalent total strain rate and k is the constant value <1. In our simulations k was chosen to be 0.8. Larger values of k had an adverse effect on iteration convergence, while smaller values led to rate-dependent effects. The calculation results for *  given by Eq. (17) are plotted in Fig. 4b and 5 for the tensile strain rates varied by three orders of magnitude. A close agreement of the averaged stress-strain curves indicates the rate-independence of the material model at the macroscale. The corresponding plastic strain fields compared in Fig. 5 at two most different strain rates support this conclusion for the microscale as well. The differences in the plastic strain patterns are almost undistinguishable either qualitatively or quantitatively. Microstructure-based Simulations of Quasistatic Deformation Using an Explicit Dynamic Approach 253 a) b) Fig. 5 Equivalent plastic strain fields calculated with ABAQUS/Explicit for the strain rates 10 2 (a) and 10 4 s -1 (b). Tensile strain is 8% Finally, let us compare the computational costs needed for solving the same 3D quasistatic problem using the ABAQUS Implicit and Explicit solvers. The estimations of the element number, which can be accommodated in the memory, were performed only for C3D8R finite elements for the initial static and explicit dynamic steps. The maximum number of elements in the explicit calculations was found to be 10-12 times larger than that in the implicit computations. Table 1 Normalized computational time in static and dynamic calculations No. of CPUs 1 2 3 4 ABAQUS/Implicit runtime 1 0.64 0.57 0.52 ABAQUS/Explicit runtime 0.16 0.08 0.06 0.04 In order to compare the runtime values of the implicit and explicit calculations, quasistatic uniaxial tension up to the same strain value was solved with ABAQUS/Standard and Explicit. The tension velocity in the dynamic calculations was chosen to provide a close agreement between the quasistatic and dynamic solutions. The calculation time values normalized to those required for nonparallel static calculations are given in Table 1 for different numbers of CPUs. The time needed for the explicit calculations is much shorter than that for the implicit ones. The dynamic calculations become even more advantageous in the case of parallel computations. 5. CONCLUSION The numerical solution to the mechanical boundary-value problem with an explicit account of the material microstructure requires substantial computational resources due to a necessity of using detailed meshes with a large number of elements. An approach that considerably reduces the requirements for computer memory, disk space, and computational time implies the solution of quasistatic problems in a dynamic formulation, where the equation of motion is solved instead of static equilibrium equation. This enables a transition 254 V. ROMANOVA, R. BALOKHONOV, E. EMELIANOVA, O. ZINOVIEVA, A. ZINOVIEV from implicit to explicit calculations providing a significant improvement of computational capacity. In this paper we have shown that the explicit dynamic approach can be successfully used in the microstructure-based simulations of quasistatic deformation, substantially reducing the computational costs without losing the information and solution accuracy. The following conditions have to be met to ensure a close agreement between the dynamic and static solutions: (i) the load velocity in the dynamic calculations must be smoothly increased to its amplitude value and then kept constant to minimize the acceleration term in the initial loading stage and (ii) the constitutive model used must be quasi-rate-independent. An examination of the mesh convergence and strain-rate dependence of the polycrystalline aluminum model has supported this conclusion. Acknowledgements: This work is supported by the Russian Academy of Sciences within the Fundamental Research Program for 2013–2020. The constitutive model presented in Section 3 has been developed within the joint research project of Deutsche Forschungsgemeinschaft (Grant No. PL 584/4-1) and the Russian Foundation for Basic Research (Grant No. 18-501-12020). REFERENCES 1. Harewood, F.J., McHugh, P.E., 2007, Comparison of the implicit and explicit finite element methods using crystal plasticity, Computational Materials Science, 39(2), pp. 481-494. 2. Kutt, L.M., Pifko, F.B., Nardiello, J.A., Papazian, J.M., 1998, Slow-dynamic finite element simulation of manufacturing processes, Computers and Structures, 66(1), pp. 1-17. 3. Hu, X., Wagoner, R.H., Daehn, G.S., Ghosh, S., 1994, Comparison of explicit and implicit finite element methods in the quasistatic simulation of uniaxial tension, Communications in Numerical Methods in Engineering, 10(12), pp. 993-1003. 4. Dimaki, A.V., Shilko, E.V., Popov, V.L., Psakhie, S.G., 2018, Simulation of fracture using a mesh-dependent fracture criterion in a discrete element method, Facta Universitatis, Series: Mechanical Engineering, 16(1), pp. 41-50. 5. Roters, F., Eisenlohr, P., Bieler, T.R., Raabe, D., 2010, Crystal plasticity finite element methods: in materials science and engineering, Wiley‐ VCH Verlag GmbH & Co. KGaA. 6. Romanova, V.A., Balokhonov, R.R., 2009, Numerical simulation of surface and bulk deformation in three-dimensional polycrystals, Physical Mesomechanics, 12(3-4), pp. 130-140. 7. Busso, E.P., Cailletaud, G., 2005, On the selection of active slip systems in crystal plasticity, International Journal of Plasticity, 21(11), pp. 2212-2231. 8. Teplyakova, L.A., Bespalova, I.V., Lychagin, D.V., 2009, Spatial organization of deformation in aluminum [1 ī 2] single crystals in compression, Physical Mesomechanics, 3(12), pp. 166-174.