图书馆订阅: Guest
国际多尺度计算工程期刊

每年出版 6 

ISSN 打印: 1543-1649

ISSN 在线: 1940-4352

The Impact Factor measures the average number of citations received in a particular year by papers published in the journal during the two preceding years. 2017 Journal Citation Reports (Clarivate Analytics, 2018) IF: 1.4 To calculate the five year Impact Factor, citations are counted in 2017 to the previous five years and divided by the source items published in the previous five years. 2017 Journal Citation Reports (Clarivate Analytics, 2018) 5-Year IF: 1.3 The Immediacy Index is the average number of times an article is cited in the year it is published. The journal Immediacy Index indicates how quickly articles in a journal are cited. Immediacy Index: 2.2 The Eigenfactor score, developed by Jevin West and Carl Bergstrom at the University of Washington, is a rating of the total importance of a scientific journal. Journals are rated according to the number of incoming citations, with citations from highly ranked journals weighted to make a larger contribution to the eigenfactor than those from poorly ranked journals. Eigenfactor: 0.00034 The Journal Citation Indicator (JCI) is a single measurement of the field-normalized citation impact of journals in the Web of Science Core Collection across disciplines. The key words here are that the metric is normalized and cross-disciplinary. JCI: 0.46 SJR: 0.333 SNIP: 0.606 CiteScore™:: 3.1 H-Index: 31

Indexed in

Adaptive Model Selection Procedure for Concurrent Multiscale Problems

Mohan A. Nuggehally

Scientific Computation Research Center, Rensselaer Polytechnic Institute, 110 8th St., Troy, NY 12180, USA

Mark S. Shephard

Scientific Computation Research Center, Rensselaer Polytechnic Institute, 110 8th St., Troy, NY 12180, USA

Catalin R. Picu

Scientific Computation Research Center, Rensselaer Polytechnic Institute, 110 8th St., Troy, NY 12180, USA

Jacob Fish

Scientific Computation Research Center, Rensselaer Polytechnic Institute, 110 8th St., Troy, NY 12180, USA

Abstract

An adaptive method for the selection of models in a concurrent multiscale approach is presented. Different models from a hierarchy are chosen in different subdomains of the problem domain adaptively in an automated problem simulation. A concurrent atomistic to continuum (AtC) coupling method [27], based on a blend of the continuum stress and the atomistic force, is adopted for the problem formulation. Two error indicators are used for the hierarchy of models consisting of a linear elastic model, a nonlinear elastic model, and an embedded atom method (EAM) based atomistic model. A nonlinear indicator , which is based on the relative error in the energy between the nonlinear model and the linear model, is used to select or deselect the nonlinear model subdomain. An atomistic indicator is a stress-gradient-based criterion to predict dislocation nucleation, which was developed by Miller and Acharya [6]. A material-specific critical value associated with the dislocation nucleation criterion is used in selecting and deselecting the atomistic subdomain during an automated simulation. An adaptive strategy uses limit values of the two indicators to adaptively modify the subdomains of the three different models. Example results are illustrated to demonstrate the adaptive method.

Key words:

adaptive model selection, concurrent multiscale, model error indicators, continuum-atomistic coupling


Video materials:
Concurrent model configuration at 60th load step ( Concurrent model configuration at 60th load step ( indentation displacement). Colors represent the subdomains in which various models are used ) results obtained with the three-model hierarchy ) results obtained with the two-model hierarch indentation displacement).
The concurrent model and dislocation loops observed at load step 18 (trace of the strain = 0.0855).


Introduction

Key to the effective application of multiscale simulations is ensuring the accuracy of the simulation. Ultimately, this accuracy must be quantified against physical measurements. However, in the case of multiscale simulations with multiple interacting models and computational methods to solve the models, it is also critical to employ methods to ensure that the desired accuracy is obtained relative to the highest fidelity model in the multiscale hierarchy. Then, the process of validation against physical measurements can be focused on validation of the highest fidelity model. The goal of adaptive multiscale methods is to provide the mechanisms needed to estimate and control the model hierarchy errors and computational methods errors. More specifically, adaptive multiscale methods must control (i) the models applied at the relevant scales over various portions of the space/time domain of the problem, (ii) the numerical methods used to solve the problem for the different models applied, and (iii) the methods used to couple the models at different scales.

Our approach to the adaptive model selection is a strict and straightforward identification of breakdown of models in a hierarchy of models chosen a priori Model hierarchy for multiscale problems typically consists of models representing different length and time scales. This paper presents an adaptive concurrent multiscale procedure that considers a hierarchy of three models. Two of them are continuum models, and the third is an atomistic model. After a brief review of some of the adaptive multiscale methods, Section 3 describes the hierarchical concurrent model that is adopted here along with a discussion of the problem formulation. Section 4 presents error indicators and an adaptive strategy for model refinement and model coarsening. The relative error in the energy is used as an indicator between the two continuum models. A criterion based on the curl of the stress is employed as an atomistic indicator. Section 5 demonstrates numerical examples for adaptive model refinement and adaptive model coarsening.

Adaptive Multiscale Approaches

Modeling error estimation and adaptive model selection has become a research topic of interest in recent years. The goal of these methods is to assess the errors in quantities of interest produced by the currently applied model as it is used to solve the given problem compared to that of a model that is known to provide a higher level of accuracy for the quantities of interest. Knowledge of these error estimates can then be used to determine which portions of the space/time domain require better models. The potential to obtain quantitative error estimates is greatest in those cases where the set of models applied is based on a well qualified-hierarchy of model improvements. For many multiscale problems of interest, such a hierarchy of models is not readily available. However, for many of these problems there is a known hierarchy of models that is known to provide the ability to capture more physics and/or account for the influence of important physical properties on finer scales. In such cases, it is often possible to define error indicators that correctly indicate when a particular coarse model is not adequate, thus indicating the next model in the hierarchy needs to be used. Some of the adaptive multiscale techniques that exist in the literature are reviewed below.

Local/Nonlocal Criterion of the Quasicontinuum Method

The quasicontinuum method (QCM) [1] is a way of solving problems of interest at the atomistic scale by dramatically reducing the cost of computation by a coarse-graining procedure. The QCM consists of coarse-graining an atomistic domain by selecting a small subset of the total number of atoms, called representative atoms. In certain critical areas of the problem domain, all the atoms are selected as representative atoms (referred to as nonlocal atoms), whereas far away from the critical regions where the deformation is more uniform at the fine scale, a group of atoms are enslaved to move under the influence of the representative atoms at the nodes of a tetrahedral element (referred to as local atoms) enclosing the group of atoms. The energy calculation for the local atoms is approximated by applying the Cauchy-Born rule [2]. A criterion based on variation in the deformation gradient is used as an error indicator to adaptively distinguish the nonlocal atoms from the local atoms [3]. The criterion evaluated for the case of linear interpolation for displacements (constant deformation gradient F) is given by [3]

(1)

where is the kth eigenvalue of the right stretch tensor in element a. Indices a, b run over all the elements in the vicinity of a node that are within three times the cutoff distance of the interatomic potential. is the transpose of the deformation gradient tensor in element a. In addition to the local/nonlocal criterion, a mesh refinement criterion based on the ZZ error indicator [4,5] is used to add and remove representative atoms adaptively.

Recently, Miller and Acharya [6] developed a stress-gradient-based dislocation nucleation criterion that could be used in QCM adaptivity. The stress-gradient-based criterion is a mathematically rigorous derivation from the continuum thermomechanical equations of the dislocation field theory [6,7]. This is the criterion that we adopt as the basis for the atomistic indicator in Section 4.2.

Modeling Error and Adaptivity

The existence of a fine model that can very closely represent the physical processes of interest enables estimation of the error of coarser models with respect to the fine model. A brief summary of the modeling error estimation and adaptivity techniques presented by Oden et al. [8-10], Larsson and Runesson [11,12], and Stein and Ohnimus [13,14] is given here.

A fine-scale mathematical model representing a physical process of interest is formulated as follows:

Given appropriate boundary conditions, find such that

(2)

where and are function spaces of admissible functions. is a residual defined in a weighted integral form, is a semilinear form that may be nonlinear in while strictly linear in . is a linear functional that represents forcing function. For the cases where it is too expensive and unreasonable (even impossible) to solve Eq. (2), a coarse-scale model with a different semilinear operator is constructed as follows:

Given the boundary conditions, find such that -->

(3)

A finite dimensional discretization of coarse model 3 results in the following:

Given the boundary conditions, find such that

(4)

with the assumption that and . The residual associated with the finite dimensional solution of the coarse-scale model (3) is obtained by substituting in Eq. (2) as follows:

(5)

where

(6)

is the finite dimensional discretization error for the coarse model considered and

(7)

is the modeling error for the finite dimensional solution of the coarse model (3) with respect to the fine-scale model (2). Procedures to control the discretization error [Eq. (6)] is well developed and discussed in detail in [15,16]. The focus here is estimation of the modeling error term.

Equation (7) for the modeling error can only be computed approximately. The first approximation is with regard to , which might not be explicitly available or might be too expensive to compute even if available. Based on the existence of a model hierarchy, a model coarser than but finer than is replaced in Eq. (7) as

(8)

Clearly, the choice of a sufficiently ``good'' in Eq. (8) is crucial to obtain an accurate . Equation (8) is further approximated by substituting for .

An automatic detection of subdomains in which the mathematical model needs to be updated is presented in [13,14], where the hierarchy of models consists of continuum models of different complexity from 3D elastoplastic models to 2D plate models. The procedures for estimating the discretization error has been extended to include the error of a current continuum model, such as a 2D plate model, with respect to a master model, such as a 3D elastoplastic model. An expansion strategy of model adaptivity is employed, where they start with a simpler model and move up the hierarchy toward more complicated models. The discretization error is distinguished from the modeling error as follows for the special case of bilinear operators instead of semilinear operators [14]. The discretization error estimation procedure starts with computation of improved equilibrated boundary traction at the element interfaces by local variational problems on an element patch. Improved displacements are obtained by again solving local variational problems with Neumann boundary conditions using an extended test function space and trial function space as follows:

Find such that

(8)

where is an element of the finite element mesh, and are locally enhanced trial and test function spaces that contain functions of different polynomial orders in different directions such that

(9)

The discretization error estimator over is given by

(10)

To compute the model error estimator, further improved displacements are obtained by solving a local variational problem using the same enhanced test and trial function spaces given in Eq. (10) but with a different bilinear form as follows:

Find such that

(11)

where is the bilinear form corresponding to the master model in the hierarchy with a more sophisticated constitutive equation than the . The model error estimator for an element is given by

(12)

Further details regarding this work of Stein and Ohnimus is available in [13,14].

The works of Oden and Prudhomme [8,17], Larsson and Runesson [11,12], Braack and Ern [18] on modeling error estimation and model adaptivity are built on the goal-oriented error estimation methods originally developed in the context of finite element discretization error estimation [15,16]. In these works, the mathematical theory of discretization error control has been extended to assess modeling error among a hierarchy of computational models of different complexity.

Often the quantity of interest is not the solution itself but a function of the solution represented as called the goal function. A goal-oriented error measure of the finite dimensional solution of the coarse model (3) is introduced as

(13)

Dual formulation for the primal problem of Eq. (2) based on the quantity of interest (14) is discussed in [8,11]. If is the solution of the dual problem the error measure (14) is computed as

(14)

In addition to the approximations mentioned for Eq. (8) in computing the modeling error , approximation is introduced in the linearization of the dual problem since it depends on the unknown exact solution and the unavailable/expensive fine model of the primal problem (2). Approximate forms are introduced in the form of a tractable model , with the corresponding solution and goal function . Another approximation is in solving for the dual solution in a finite dimensional approximation space as [11].

Goal-oriented model error estimation and model adaptivity has been used with continuum model hierarchies in [9,11,12,17,18]. Recently, this method has been extended to concurrent multiscale problems with coupled continuum-atomistic models by Oden [10,19], in which application to the quasicontinuum method (QCM) [3] and bridging scale method (BSM) of Wagner and Liu [2] are discussed. For both applications, a fully atomistic model is the fine-scale model. QCM and BSM are used as coarse models to solve the molecular statics and the molecular dynamics problems respectively. Details of the problem formulation of primal, dual, and a goal function are given in [10]. As mentioned in [10], the error estimate is the estimate of the total error and does not distinguish between the discretization and modeling errors. In fact, the error estimate computed in [10] is used to refine the QCM mesh by replacing its mesh refinement criterion based on the ZZ error indicator [4,5]. The local/nonlocal criterion of Eq. (1) is still employed to adaptively distinguish the local and nonlocal atoms. It seems like the local/nonlocal criterion is the one that controls the model error while the QCM mesh refinement criterion controls the discretization error. For the BSM application, a goal-oriented error estimate is used to turn coarse-grained regions of the subdomain into a fully atomistic subdomain.

Hierarchical Concurrent Model

A model hierarchy consisting of an atomistic model and two continuum models is adopted here. The finest model in the hierarchy is the atomistic model with the embedded atom method (EAM) [21] based interatomic potential. The other two coarser-scale models are a nonlinear elastic continuum model, which has terms up to third-order elastic constitutive parameters, and a linear elastic continuum model. These continuum models fail to capture inhomogeneity due to the nucleation of defects unless mixed continuum-discrete models [22] are included in the hierarchy.

The nucleation of defects and their interactions at the atomistic level provide a fundamental understanding of the mechanical response of materials at the nanoscale that influences the macroscale response. Crystalline material behavior is elastic in the absence of dislocation nucleation and motion. A homogeneous elastic deformation prior to the nucleation of defects can be closely captured by a nonlinear elastic continuum model. The third-order elastic constants are found to be essential for a proper treatment of the finite deformation of crystals [23] and are also used extensively in the study of elastic waves in crystals. A brief overview of the models in the hierarchy is given below followed by the concurrent problem formulation.

EAM Atomistic Model

For the EAM interatomic potential, the total energy of a system of atoms obtained as the sum of energies of individual atoms is given as

(15)

where is

(16)

(17)

where is the total electron density at atom , is the embedding energy function. is the distance between the atoms and . is the pair potential term and is the electron density function, which have a cutoff distance in terms of as defined by the interatomic potential. Thus, the summation in Eqs. (17) and (18) is over the atoms in a neighborhood of the atom denoted by . Greek letters denote atoms and are used as subscripts to represent the quantities related to atoms. Also, there is no summation convention on the Greek subscripts.

Note, that the specific choice of EAM is not a limitation for the adaptive model selection discussed here. The EAM has been used to model crystal defects, surfaces of solids and liquids, etc. [24], and is also well suited for the purpose of accounting for the strain energy of some metallic crystals.

Nonlinear and Linear Elastic Continuum Models

The constitutive relationship for the nonlinear and linear elastic continuum models are obtained by the hypothesis of the existence of a stored internal energy density as a function of the Green-Lagrange strain . For an elastic solid, the following constitutive relationship can be derived by defining the Cauchy stress as the work conjugate of [25] as

(18)

where is the deformation gradient tensor. Uppercase subscripts denote components of a tensor in the reference or undeformed configuration and lowercase subscripts denote components of a tensor in the current or deformed configuration. To obtain a complete nonlinear elastic continuum theory, is expanded as a power series in as follows:

(19)

Thus, the third-order elastic constants are coefficients of the cubic terms in the energy expansion Eq. (20). is the second-order elastic tensor used in the linear theory. are chosen to be zero for a stress-free initial state.

Truncating the power series of Eq. (20) by retaining up to third-order terms, we obtain the constitutive relation for the nonlinear elastic model as

(20)

Retaining only up to second-order terms results in the constitutive equation for the linear elastic model as

(21)

The constitutive parameters and are chosen to be consistent with the EAM potential as discussed in [26].

Concurrent Model Problem Formulation

The problem domain shown in Fig. 1 is composed of subdomains in which various models from the hierarchy are used. The linear elastic description is used in , the nonlinear description is defined for , the EAM atomistic description is defined for and is the interphase subdomain, which has a blended description of the nonlinear elastic continuum and the atomistic models. is the continuum subdomain. The details of the problem formulation are given in [27]. The governing equations are standard equilibrium equations in the continuum region, the equilibrium of forces for each atom in the atomistic region, and a blend of the continuum stress and the atomistic force in the interphase region. Key points are given here.

Hybrid concurrent model domain

Figure 1. Hybrid concurrent model domain

The weak form of the equilibrium equation is stated as follows:

Given , , ,

Find displacements and such that

(24)

and a weak displacement compatibility is satisfied as follows:

(25)

and are continuum function spaces defined as follows:

(26)

(27)

where in Eqs. (25) and (26) denote Hilbert spaces [28]. and belong to the discrete phase space of the atomistic system given by

(28)

(29)

in Eq. (24) is the phase space of atoms in the interphase .

The terms appearing in Eq. (23) are explained below. and are the body forces for the continuum and atomistics, respectively. and are the essential and natural boundary conditions on essential boundary and natural boundary , respectively, where the domain boundary and . The internal force acting on atom due to the atom is given by

(29)

where is the displacement of the atom . Blend functions and are defined as follows:

(30)

The continuum blend function on , on and on , which is evaluated based on the proximity of the point to the boundaries and . By defining as a normalized distance in the physical domain from to , can be approximated to be a function of the scalar parameter . is the Dirac delta function whose value at is infinity and integral over is .

Discretized Constraint and Equilibrium Equations

A finite element discretization of the problem domain is denoted by . The discrete compatibility equation is constructed by choosing in Eq. (24) using piecewise constant shape functions defined to be constant over the finite element domains as follows:

(31)

where

(32)

Substituting Eqs. (31) and (32) into Eq. (24), we obtain

(33)

Requiring the arbitrariness of yields the following discrete compatibility equation for every element , being the number of atoms in an element :

(34)

Equation (34) yields number of constraint equations equal to the number of spatial dimension for each finite element . From Eq. (34), the degrees of freedom of one atom in can be expressed in terms of the degrees of freedom of the finite element and the degrees of freedom of the remaining atoms in the element. Note that at least one atom has to be positioned with an element in the interphase.

The continuum displacement and test functions defined over are discretized using continuous finite element shape functions. The discretized displacement is denoted by and the discretized test function is denoted by . The spaces and are given by

(35)

(36)

where are the finite element shape functions associated with the finite element nodes , are the nodal degrees of freedom and are the nodal multipliers corresponding to test functions. Summation convention over repeated index is employed. The discretized system of equations shown below is obtained by using Eqs. (35) and (36) in the equilibrium Eq. (23)

(37)

where is the number of independent atoms. Equation (37) is a nonlinear system of equations in continuum degrees of freedom and independent atomistic degrees of freedom. Assuming no follower forces and no body forces, Eq. (37) can be written in terms of residuals as

(38)

is expressed in terms of the continuum degrees of freedom according to appropriate constitutive equations (for example, Section 3.2) and is expressed in terms of the independent atomistic degrees of freedom [27]. The nonlinear system of Eq. (38) can then be solved either by Newton method or by conjugate gradient minimization of the residuals. Several options for the blend function are discussed in [27].

Adaptive model selection procedure

In order to use the hierarchy of models outlined in Section 3 effectively in a concurrent multiscale problem setup, we need an automated technique to identify subdomains of the problem domain to properly apply different computational models. An adaptive model selection strategy is to be devised such that we capture the relevant physics in the essential areas with a fine-scale model and use coarser models elsewhere in the problem domain.

Two error indicators are adopted to facilitate the adaptive refinement and adaptive coarsening of the models in the hierarchy. Relative error in the energy computed by the two continuum models over an element of a finite element discretization is used as an indicator between the nonlinear model and the linear model, which is discussed in Section 4.1. A stress-gradient-based dislocation nucleation criterion [6] is adopted as the atomistic indicator and is reviewed in Section 4.2. An adaptive strategy discussed in Section 4.3 employs the error indicators for adaptive model section and deselection.

Nonlinear Error Indicator

The strain energy computed over a subdomain by the linear elastic model and the nonlinear model (Section 3.2) are used to construct an indicator to distinguish the boundary of the nonlinear and linear continuum model subdomains. Specifically, the relative error in energy between the two models is the linear to nonlinear indicator, which is given by

(39)

where and are, respectively,

(40)

(41)

is an element of a finite element discretization of the problem domain . Note that is approximated with a minimal computational effort by using the strain field of the linear elastic model and , which is calibrated [26] based on the interatomic potential used for the finest model.

In comparison to the global modeling error estimation techniques developed by Oden and Prudhomme [8] and Oden [9], the elemental error indicator of the linear elastic model with respect to the nonlinear model is given by

(42)

which is related to the global modeling error by [9]

(43)

where the summation is over all the elements in the mesh . Comparing Eqs. (39) and (42), we can say that and are related.

Computing the Nonlinear Error Indicator Limit Value

A limit value of , denoted by , is computed to facilitate the selection of the nonlinear model during automated adaptive simulation. The relative error in the energy density of the linear elastic model with respect to the atomistic model is the basis for computing ; a procedure for this is explained below.

Relative error in the energy density versus load step

Figure 2. Relative error in the energy density versus load step

A rectangular block is subjected to uniaxial tension. The load is applied quasi-statically through Dirichlet boundary conditions. The constant energy density is computed for each load step, by choosing different models for the problem domain. The atomistic energy density is , where is given by Eq. (17) and is the atomic volume associated with the unit cell of atom . Figure 2 shows a plot of the relative error in the energy density of the linear model and the nonlinear model with respect to the atomistic model. is computed at shown in Fig. 2.

(44)

Figure 2 also shows the effectiveness of the nonlinear model in reducing the error in the energy with respect to the atomistic model. In other words, using the nonlinear model in the hierarchy reduces the size of the atomistic domain for a given error tolerance. This is illustrated in Section 5.1 for the nanoindentation example. is employed in the adaptive strategy (Section 4.3) to automatically select subdomains that are above this limit value and assign them the nonlinear model.

Atomistic Error Indicator

Dislocation motion and close-range interaction cannot be captured by continuum models. However, these phenomena are crucial for understanding incipient plasticity in materials. The atomistic model is used to capture dislocation nucleation and the initial stages of their interaction. A stress-gradient-based dislocation nucleation criterion [6] is adopted as the basis of the atomistic indicator. The criterion is derived from the thermomechanical considerations of a crystalline solid undergoing deformation. It is developed as a "driving force" or "work conjugate" of the dislocation density tensor and hence expresses the thermodynamic driving force, for an increase of the dislocation density. In the absence of this driving force the nucleation event ceases to exist.

The derivation of the criterion is based on an analysis of dissipation terms that identify the driving forces corresponding to the dissipative mechanisms according to the second law of thermodynamics. The mechanics and derivation details were described by Acharya in references [7,28]. Dislocation nucleation and dislocation motion are the two dissipative mechanisms under consideration according to the dislocation field theory presented in [6], which says that the rate of work done by the stress can be partitioned into a stored free energy rate and dissipated energy rate due to dislocation motion and dislocation nucleation as

(45)

where is the velocity of the matter ( ), is the driving force for dislocation motion with a velocity field . is derived from and is the nucleation rate tensor written in terms of the density of the Burger's vector per unit area per unit time of a particular type of dislocation as

(46)

where is the direction of Burger's vector and is the dislocation line direction. Thus, the last term in Eq. (45) reduces to

(47)

The quantity is the driving force for the variation of the dislocation density . It is shown in [6] that can be approximated by . The dislocation nucleation criterion is based on , which is given by

(48)

In order to predict the dislocation nucleation, theoretically, the in Eq. (48) needs to be evaluated for all possible values of the Burger's vector and the dislocation line direction at each point in the domain. The maximum value indicates the location and type of the dislocation expected to nucleate. In practice, a finite number of candidates for and are considered for a given crystal lattice type (for example, FCC lattices have 12 slip systems). As for the spatial location, they are evaluated at Gauss integration points of the elements of a finite element mesh. For 2D problems, specific values of and are known .

Atomistic Indicator Tolerance

The existence of a material-specific critical value for the dislocation nucleation criterion has been hypothesized in [6]. The criterion will predict nucleation of a specific dislocation characterized by and at the point when . To test the prediction capability of the criterion, needs to be evaluated at all points in the domain and for all possible combinations of and at each point. The maximum value of evaluated just before the nucleation would quantify . However, from a computationally practical perspective, the knowledge of slip systems of a crystal as the most likely candidates for the dislocation nucleation is utilized in [6] to quantify . Its prediction capability and material-property-like behavior is verified with atomistic model examples in [6].

The accuracy of is another factor affecting the prediction capability of the criterion. It is crucial to control the discretization error to the desired accuracy to obtain an accurate solution for . Because the finite element mesh of the problem domain is currently constructed and the mesh elements are large compared to the atomic spacing, we adopt a fraction of denoted by as the atomistic indicator tolerance. The specific value of the tolerance value used is given in the examples Section 5. Future work will be required to better address this process so that a less conservative value of can be used.

is similar to the local/nonlocal criterion of the quasicontinuum method in the sense that both of these are nonlocal in nature. [Eq. (48)] accounts for the nonlocality in terms of , where as the local/nonlocal criterion [Eq. (1)] of the quasicontinuum method discussed in Section 2.1 accounts for nonlocality in terms of variation in the deformation gradient of the elements in the neighborhood of a node. In the goal-oriented error estimation approach of Oden applied to multiscale problems [10,19], the total error estimate (discretization + modeling) of a goal function is used to refine a coarse mesh down to atomistic size in critical regions. However the local/nonlocal criterion [Eq. (1)] is still used to distinguish the local atoms from the nonlocal atoms [10,19].

Automated Adaptive Procedure for Model Refinement and Coarsening

An adaptive procedure that incorporates the error indicators for the models of the hierarchy is presented here. The adaptive procedure automatically assigns different models to different subdomains of a problem domain based on the error indicators. Tolerance values corresponding to and corresponding to the atomistic indicator are chosen as discussed in Sections 4.1.1 and 4.2.1. The sequence of steps in the adaptive procedure are as follows:

  1. Start the load step of a quasi-statically loaded problem with a coarse model (Note: Scale the load to a value where more than just a linear elastic model is needed to solve the problem).

  2. lncrement the load step .

  3. Solve the problem with the concurrent model problem formulation of Section 3.3.

  4. Compute appropriate error indicators for the model hierarchy chosen.

  5. Select or deselect the models from the hierarchy for subdomains of the problem domain according to the error indicators and corresponding tolerance values.

  6. If models change, go to step 3, else continue with step 2 till the end of load steps.

Example results

Examples illustrating the adaptive model refinement and adaptive model coarsening are presented here. The constitutive parameters of the continuum nonlinear and linear elastic models are consistent with the EAM potential of aluminum [29]. The cubic elastic constants for aluminum are , , and in the principal crystallographic direction. The nonlinear indicator tolerance is chosen as as discussed in Section 4.1.1, and the atomistic indicator tolerance is chosen as . The specific value of is chosen based on numerical experiments. For the solution accuracy of meshes used in the examples, the value of was calibrated to predict dislocation nucleation correctly with respect to a fully atomistic model. This value of works for the different examples in this section in accordance with the claim of Miller and Acharya [6] that the is a material specific value and does not depend on the loading.

The problem domain is discretized with tetrahedral elements such that the mesh is finer in the areas of high stress gradients. Adaptive discretization error control through a mesh modification procedure will be added in the near future. Quadratic shape functions are used to discretize displacement and displacement test function of Eqs. (35) and (36), respectively. The quadratic approximation of the displacement yields a piecewise linear stress field, which is the minimum polynomial order expected to compute stress gradient in the dislocation nucleation criterion [Eq. (48)].

The two finer models (Section 3) are selected or deselected for each element of the finite element mesh based on the error indicators (Sections 4.1 and 4.2) and the corresponding tolerance values. While computing for 3D problems, there are 12 slip systems for FCC aluminum crystal that needs to be considered. The dislocation line direction is searched within each slip system to obtain maximum value of .

The first example is a nanoindentation example that demonstrates the effectiveness of the three-model hierarchy over the two-model hierarchy. The second example is a nanovoid subjected to hydrostatic tension up to a certain load, after which it is unloaded steadily back to the initial state. The unloading stage demonstrates adaptive coarsening in which a fine-scale model subdomain is transformed into a coarse-scale model subdomain during the course of simulation. The last example is the case of multiple nanovoids (four in this case) subjected to cavitation.

Nanoindentation

A comparison of the three-model hierarchy, consisting of the atomistic model, the nonlinear model and the linear elastic model, with that of the two-model hierarchy, consisting of the atomistic model and the linear elastic model, is given here. The problem domain is a film of thickness placed on a rigid substrate. The indenter is rectangular in shape and wide. The indenter as well as the film are considered to be infinite in the direction (out of the plane); thus, a plane strain condition exists in the plane (Fig. 3). Homogeneous Dirichlet boundary conditions are imposed in the direction on the bottom face, in the direction on the left and right face of the continuum domain. The indenter load is applied quasi-statically through Dirichlet boundary condition by moving the indenter by for each load step. Periodic boundary conditions in the direction are imposed on the atomistic model to maintain the 3D lattice structure. The crystallographic orientation chosen (Fig. 3) is such that the dislocations generated from the corners of the indenter move in the negative direction. Thus, and are used to evaluate .

Concurrent model configuration at 60th load step ( indentation displacement). Colors represent the subdomains in which various models are used ) results obtained with the three-model hierarchy ) results obtained with the two-model hierarch

( Concurrent model configuration at 60th load step ( indentation displacement). Colors represent the subdomains in which various models are used ) results obtained with the three-model hierarchy ) results obtained with the two-model hierarch )

Figure 3. Concurrent model configuration at 60th load step ( Concurrent model configuration at 60th load step ( indentation displacement). Colors represent the subdomains in which various models are used ) results obtained with the three-model hierarchy ) results obtained with the two-model hierarch indentation displacement). Colors represent the subdomains in which various models are used a) results obtained with the three-model hierarchy b) results obtained with the two-model hierarch

Energy comparison between atomistic reference solution and concurrent models with different hierarchy

Figure 4. Energy comparison between atomistic reference solution and concurrent models with different hierarchy

Figure 3 shows the concurrent models at the 60th load step ( A indenter displacement) for the two cases of simulation. The one on left corresponds to a simulation with the three-model hierarchy, and the one on right is a simulation with the two-model hierarchy. The concurrent model with the three-model hierarchy has atoms, whereas the concurrent model with the two-model hierarchy has atoms. That is less degrees of freedoms to solve for and that many less residuals to compute. The neighbor list building operation for the extra atoms is also saved. Figure 4 shows the variation of the total strain energy in the model during the loading period. The discontinuities in the curve correspond to the nucleation of partial dislocations from the indenter corners. The two models predict an identical defect nucleation sequence and almost identical total energy, which demonstrates the usefulness of employing the three-model hierarchy.

Nanovoid Subjected to Cavitation

Adaptive model refinement and model coarsening is demonstrated here for the case of a nanovoid subjected to a cycle of hydrostatic tensile loading-unloading. The problem domain is a cube of side A with a spherical void of A diameter at the center. The load is applied quasi-statically in small increments of A on each face of the cube through Dirichlet boundary conditions. The simulation of a nanovoid subjected to hydrostatic load has been used to study nanovoid growth and cavitation in [30]. During the first 23 load steps, the problem domain is gradually loaded hydrostatically, whereas in the next 23 load steps it is unloaded gradually as per the adaptive procedure discussed in Section 4.3.

Figure 5 shows the concurrent model at load step 18 (trace of the strain = 0.0855) and the dislocation loops that just nucleate from the void surface. The atomistic model consists of 20,190 atoms.

The concurrent model and dislocation loops observed at load step 18 (trace of the strain = 0.0855)

( The concurrent model and dislocation loops observed at load step 18 (trace of the strain = 0.0855) )

Figure 5. The concurrent model and dislocation loops observed at load step 18 (trace of the strain = 0.0855)

The concurrent model and stacking fault tetrahedra around the void at load step 23 (trace of the strain = 0.1035)

Figure 6. The concurrent model and stacking fault tetrahedra around the void at load step 23 (trace of the strain = 0.1035)

These dislocation loops grow and react to form Lomer-Cottrell junctions [31] as seen in the right side of Fig. 6. Stacking fault tetrahedra form around the void. Atoms along the edges of these tetrahedra are seen on the right side of Fig. 6 (energetic atoms located in the core of dislocation), whereas the left side shows the concurrent model at that load step. The atomistic model at this load step consists of 69,070 atoms. The symmetry of the resulting dislocation configuration is due to the crystal symmetry. The crystallographic orientation is also shown in Fig. 6. The model is unloaded gradually back to the initial stage. The reaction products of the dislocation structures formed during the unloading stage remain in the model as debris and can be seen on the right side of Fig. 7. Also note that the atomistic model has shrunk in size compared to that at maximum load. The left portion of Fig. 7 shows the concurrent model at the end of the simulation and consists of 55,763 atoms.

The concurrent model and reaction products of stacking fault tetrahedra after unloading-load step 46 (trace of the strain = 0.0)

Figure 7. The concurrent model and reaction products of stacking fault tetrahedra after unloading-load step 46 (trace of the strain = 0.0)

Comparison of the energy between the atomistic reference solution and the concurrent model solution for the nanovoid simulation

Figure 8. Comparison of the energy between the atomistic reference solution and the concurrent model solution for the nanovoid simulation

Figure 8 shows a comparison of the energy plots between the fully atomistic reference solution and the concurrent model solution. The energy of the atoms within the interatomic cutoff distance from a free surface is subtracted from the model energy to eliminate the energy fluctuation due to surface relaxation. The linear elastic strain energy of the model is also subtracted from the total energy of the model so that the fluctuation in the energy due to nucleation and growth of dislocation loops is clearly distinguishable. Different events noted during the simulation are marked. The bottom portion of the curve corresponds to the loading cycle, and the top portion corresponds to that of the unloading cycle. The higher energy at the end of the simulation is due to the reaction products of dislocation structure that are left behind during unloading.

In another simulation of the same problem and loading conditions, loading cycle was stopped at load step 19 and unloaded for the next 19 load steps to bring it back to the initial unloaded state. It was observed that the dislocation structure that forms in this simulation is fully reversible. At the end of the simulation, the problem domain was free of dislocations and, hence, the adaptive model selection algorithm mandates a linear elastic model through out the problem domain. The two finer models were deselected gradually during the unloading part of the simulation. The full elimination of dislocation from a model that was not loaded up to the formation of a perfect stacking fault was verified using a fully atomistic model.

Multivoid Subjected to Hydrostatic Load

A block of material with four spherical nanovoids is subjected to hydrostatic tensile load. This is relevant to the study of porous materials at nanoscale. Each void is A in diameter, and the block is A cube. The load is applied quasi-statically through Dirichlet boundary conditions in increments of A of each face of the cube. The simulation is run adaptively as discussed in Section 4.3.

The concurrent model and the dislocation structures around voids-load step 20 (trace of strain = 0.0945)

Figure 9. The concurrent model and the dislocation structures around voids-load step 20 (trace of strain = 0.0945)

The left portion of Fig. 9 shows a concurrent model with different model subdomains during an adaptive simulation, while the right portion shows the high energy atoms that form stacking fault tetrahedra around voids.

Results shown in Fig. 9 reiterates dislocation loop formation around the voids in the presence of hydrostatic load, their growth, and reaction to form Lomer-Cottrell junctions and formation of stacking fault tetrahedra around the voids. The system goes through the same stages as in Section 5.2 as long as the dislocation structures do not interact. At larger loads, these start interacting with each other, new dislocations nucleate from the voids, and the overall dislocation configuration becomes much more complicated, as shown in Fig. 10.

Advanced dislocation structure around voids-load step 28 (trace of strain = 0.126)

Figure 10. Advanced dislocation structure around voids-load step 28 (trace of strain = 0.126)

This example illustrates the preliminary results of the problem of multiple interacting voids. The motivation here is to demonstrate the utility of the adaptive method to study bigger problems that are of interest for real-life applications. One interesting application ensuing the simulation of the multivoid problem is the study of the effect of the size of nanovoids on the change of porosity of the material. When the nanovoids are close to each other, they begin to interact faster and tend to coalesce to make the material more porous. Detailed study of such effects can be analyzed with the help of the adaptive method presented here.

Closing remarks

An adaptive multiscale model selection procedure is demonstrated with a model hierarchy consisting of a linear elastic model, a nonlinear elastic model, and an EAM-based atomistic model. The nonlinear error indicator based on the relative error in the energy between the nonlinear model and the linear model is used for the selection and deselection of nonlinear model subdomain. The atomistic indicator is based on the stress gradient criterion for dislocation nucleation and is used in the selection and deselection of the atomistic model subdomain. An adaptive strategy is devised to enable an automated simulation. The nanoindentation example illustrates the reduction in the size of the atomistic domain (translates to reduction in the computational cost) by using the three-model hierarchy instead of the two-model hierarchy. It was noted that the accuracy of the solution remained acceptable. The nanovoid examples show the adaptive refinement and coarsening and begin to demonstrate the ability of the method to support the study of important materials problems.

Acknowledgments

We acknowledge the support of NSF for Grant No. 0303902 for ``NIRT: Modeling and Simulation Framework at the Nanoscale: Application to Process Simulation, Nanodevices and Nanostructured composites'' and for Grant No. 0310596 for ``Multiscale Systems Engineering for Nanocomposites.'' We also acknowledge the support of the DOE program ``A Mathematical Analysis of Atomistic to Continuum Coupling Methods, '' Grant No. DE-FG01-05ER05-16. The computations performed were done on parallel computers provided by the NSF Grant No. 0420703 entitled ``MRI: Acquisition of Infrastructure for Research in Grid Computing and Multiscale Systems Computation.''

REFERENCES

  1. Tadmor, E.В B., Ortiz, M., and Phillips, R., Quasicontinuum analysis of defects in solids. Philos. Mag. A. 73(6):1529-1563, 1996.

  2. Ericksen, J.В L., The cauchy and born hypotheses for crystals. In Phase Transformations and Material Instabilities in Solids. Gurtin, M.В E. Ed., Mathematics Research Center, The University of Wisconsin-Madison, October 1983.

  3. Miller, R.В E., and Tadmor, E.В B., The quasicontinuum method: Overview, applications and current directions. J. Comput.-Aided Mater. Des. 9(3):203-239, 2002.

  4. Zienkiewicz, O.В C., and Zhu, J.В Z., The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique. Int. J. Numer. Methods Eng. 33:1331-1364, 1992.

  5. Zienkiewicz, O.В C., and Zhu, J.В Z., The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity. Int. J. Numer. Methods Eng. 33:1365-1382, 1992.

  6. Miller, R.В E., and Acharya, A., A stress-gradient based criterion for dislocation nucleation in crystals. J. Mech. Phys. Solids. 52:1507-1525, 2004.

  7. Acharya, A., Constitutive analysis of finite deformation field dislocation mechanics. J. Mech. Phys. Solids. 52:301-316, 2004.

  8. Oden, J. T., and Prudhomme, S., Estimation of modeling error in computational mechanics. J. Comput. Phys. 182:496-515, 2002.

  9. Oden, J. T., Prudhomme, S., Hammerand,В D.В C., and Kuczma, M. S., Modeling error and adaptivity in nonlinear continuum mechanics. Comput. Methods Appl. Mech. Eng. 190:6663-6684, 2001.

  10. Oden, J. T., Prudhomme, S., Romkes, A., and Bauman, P., Multi-scale modeling of physical phenomena: Adaptive control of models. Technical ReportВ No. 20, Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, May, 2005.

  11. Larsson, F., and Runesson, K., Modeling and discretization error in hyperelasto-(visco-) plasticity with a view to hierarchical modeling. Comput. Methods Appl. Mech. Eng. 193:5283-5300, 2004.

  12. Larsson, F., and Runesson, K., Adaptive computational meso-macro-scale modeling of elastic composites. Comput. Methods Appl. Mech. Eng. 195:324-338, 2006.

  13. Stein, E., and Ohnimus, S., Coupled model- and solution-adaptivity in the finite-element method. Comput. Methods Appl. Mech. Eng. 150:327-350, 1997.

  14. Stein, E., and Ohnimus, S., Anisotropic discretization and model-error estimation in solid mechanics by local neumann problems. Comput. Methods Appl. Mech. Eng. 176:363-385, 1999.

  15. Ainsworth, M., and Oden, J. T., A Posteriori Error Estimation in Finite Element Analysis. Wiley, 2000.

  16. Becker, R., and Rannacher, R., An optimal control approach to a posteriori error estimation in finite element methods. Acta Numer. 10:1-23, 2001.

  17. Oden, J. T., and Vemaganti, K. S., Estimation of local modeling error and goal-oriented adaptive modeling of heterogeneous materials. J. Comput. Phys. 164:22-47, 2000.

  18. Braack, M., and Ern, A., A posteriori control of modeling errors and discretization errors. Multiscale Model. Simul. 1(2):221-232, 2003.

  19. Oden, J. T., Prudhomme, S., and Bauman, P., On the extension of goal-oriented error estimation and hierarchical modeling to discrete lattice models. Comput. Methods Appl. Mech. Eng. 194:3668-3688, 2005.

  20. Wagner, G. J., and Liu, W. K., Coupling of atomistic and continuum simulations using a bridging scale decomposition. J. Comput. Phys. 190:249-274, 2003.

  21. Daw., M. S., and Baskes, M.В I., Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals. Phys. Rev. B. 29(12):6443-6453, 1984.

  22. Miller, R. E., Shilkrot, L.В E., and Curtin,В W.В A., A coupled atomistics and discrete dislocation plasticity simulation of nanoindentation into single crystal thin films. Acta Mater. 52:271-284, 2004.

  23. Musgrave, M.В J.В P., Crystal Acoustics: Introduction to the Study of Elastic Waves and Vibrations in Crystals. Holden-Day, San Francisco, 1970.

  24. Liu, X.-Y., Ohotnicky, P.В P., Adams, J.В B., Rohrer,В C.В L., and Hyland, RВ . WВ ., Jr., Anisotropic surface segregation in al-mg alloys. Surf. Sci. 373:357-370, 1997.

  25. Eringen, C. A., Mechanics of Continua. Chaps. 4 and 5, Wiley, New York, 1967.

  26. Datta, D. K., Picu, C. R., and Shephard, M. S., Composite grid atomistic continuum method: An adaptive approach to bridge continuum with atomistic analysis. Int. J. Multiscale Comput. Eng. 2:401-420, 2004.

  27. Fish, J., Nuggehally, M. A., Shephard, M. S., Picu, C. R., Badia, S., Parks, M. L., and Gunzburger, M., Concurrent atc coupling based on a blend of the continuum stress and the atomistic force. Comput. Methods Appl. Mech. Eng. 196:4548-4560, 2007.

  28. Acharya, A., A model of crystal plasticity based on the theory of continuously distributed dislocations. J. Mech. Phys. Solids. 49:761-784, 2001.

  29. Ercolessi, F., and Adams, J.В B., Interatomic potentials from first-principles calculations: The force-matching method. Europhys. Lett. 26:583-588, 1994.

  30. Marian,В J., Knap,В J., and Ortiz,В M., Nano-void cavitation by dislocation emission in aluminum. Phys. Rev. Lett. 93(165503):1-4,2004.

  31. Hull, D., and Bacon, D. J., Introduction to Dislocations. Butterworth Heinemann, 2001.

Begell Digital Portal Begell 数字图书馆 电子图书 期刊 参考文献及会议录 研究收集 订购及政策 Begell House 联系我们 Language English 中文 Русский Português German French Spain