Hierarchical multiscale quantification of material uncertainty
Burigede Liu, Xingsheng Sun, Kaushik Bhattacharya, Michael Ortiz
HHierarchical multiscale quantification of material uncertainty
Burigede Liu ∗ † , Xingsheng Sun ∗ † , Kaushik Bhattacharya ∗ , Michael Ortiz ∗ February 2021
Abstract
The macroscopic behavior of many materials is complex and the end result of mechanisms thatoperate across a broad range of disparate scales. An imperfect knowledge of material behavior acrossscales is a source of epistemic uncertainty of the overall material behavior. However, assessing thisuncertainty is difficult due to the complex nature of material response and the prohibitive computationalcost of integral calculations. In this paper, we exploit the multiscale and hierarchical nature of materialresponse to develop an approach to quantify the overall uncertainty of material response without theneed for integral calculations. Specifically, we bound the uncertainty at each scale and then combinethe partial uncertainties in a way that provides a bound on the overall or integral uncertainty. Thebound provides a conservative estimate on the uncertainty. Importantly, this approach does not requireintegral calculations that are prohibitively expensive. We demonstrate the framework on the problem ofballistic impact of a polycrystalline magnesium plate. Magnesium and its alloys are of current interestas promising light-weight structural and protective materials. Finally, we remark that the approach canalso be used to study the sensitivity of the overall response to particular mechanisms at lower scales ina materials-by-design approach.
Keywords
Material uncertainty; Multiscale modeling; Rigorous uncertainty quantification; Materials-by-design
The macroscopic behavior of many materials is complex and the end result of mechanisms that operateacross a broad range of disparate scales [Phillips, 2001]. The mesoscopic scales both filter (average) andmodulate (set the boundary conditions or driving forces for) the mechanisms operating at lower scales,which establishes a functional hierarchy among mechanisms. The complexity of the material response isoften a main source of uncertainty in engineering applications, a source that is epistemic in nature andtraceable to imperfect knowledge of material behavior across scales. This epistemic uncertainty often rendersdeterministic analysis of limited value. Instead, integral uncertainties must be carefully quantified in orderto identify adequate design margins and meet design specifications with sufficient confidence. However, thedirect estimation of integral material uncertainties entails repeated calculations of integral material responseaimed at determining worst-case scenarios at all scales resulting in the largest deviations in microscopicbehavior. Such integral calculations are almost always prohibitive and well beyond the scope of present-daycomputers.The complexity of materials response also poses a challenge to the development of new material systemsand the optimization of properties through a ‘materials-by-design’ approach [Olson, 2000]. In this approach, ∗ Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA 91125, United States † These authors contribute equally. a r X i v : . [ phy s i c s . c o m p - ph ] F e b ne seeks to ‘design’ materials with desired properties by targeting individual mechanisms at lower scales.However, while one can affect individual mechanisms and assess the response at a particular scale, forexample by adding solutes to affect dislocation kinetics and assessing it using molecular dynamics [Kohleret al., 2004], it is extremely challenging to understand how these changes percolate through the hierarchy. Inother words, the complexity makes it extremely difficult to understand the sensitivity of the overall materialresponse to individual mechanisms. Once again, direct estimation leads to computational problems that arewell beyond the scope of present-day computers. While our current work concerns uncertainty, there is aclose connection between these issues.Conveniently, the very multiscale and hierarchical nature of material response itself can be exploited forpurposes of uncertainty quantification. We recall that multiscale modeling is fundamentally a ’divide-and-conquer’ paradigm whereby the entire range of material behaviors is divided into a hierarchy of length scales[Ortiz et al., 2001]. The relevant unit mechanisms are then identified, namely, physical mechanisms that areirreducible and operate roughly independently: two mechanisms that are tightly coupled should be consideredas a single unit mechanism. In this hierarchy, the unit mechanisms at one scale represent averages of unitmechanisms operating at the immediately lower length scale. This functional relation introduces a partialordering of mechanisms that defines a directed graph. The nodes of the graph are the unit mechanisms,the root represents the integral macroscopic behavior of the solid and the edges define the upward flow ofinformation from the leaves to the root of the graph.A representative hierarchy for modeling strength in metals is shown in Fig. 1 by way of example. A numberof fundamental properties, such as the equations of state describing the compressibility of the material, theelastic moduli and heat capacity, can be calculated from first principles at the quantum-mechanical level up tohigh pressures and temperatures. Some transport properties, such as viscosity and thermal conductivity, canalso be characterized from first principles, e. g., within the Green-Kubo formalism [Green, 1954, Kubo, 1957].At the nanoscale, the properties of individual lattice defects, such as vacancies and dislocations, come intofocus [Leibfried and Breuer, 2006]. Such properties include dislocation core energies, kink mobilities, latticefriction, short-range dislocation-dislocation interactions, vacancy formation and migration energies, grain-boundary energies, and others. The sub-micron scale is dominated by dislocation dynamics [Messerschmidt,2010], i. e., the cooperative behavior of large dislocation ensembles, and, in the case of ductile fracture,void nucleation and cavitation. The sub-grain scale is characterized by the formation of highly-structureddislocation patterns and, in the case of metals undergoing solid-solid phase transitions, martensitic structures[Bhattacharya, 2003]. This intermediate scale is important, as it underlies scaling properties such as Hall-Petch scaling, or the inverse relation between strength and the square-root of the grain size [Hall, 1951, Petch,1953]. At the microscale, ductile fracture is characterized by plastic void growth and coalescence. Finally,the macroscopic response of polycrystalline metals represents the effective behavior of large ensembles ofgrains (cf., e. g., [Hutchinson, 1970, Wei and Anand, 2004] for reviews).In this work we show how, in materials for which a multiscale hierarchy is well defined, the quantificationof integral uncertainties can be reduced to the analysis of each unit mechanism in turn and the propaga-tion of unit uncertainties across the multiscale hierarchy according to an appropriate measure of interactionbetween the unit mechanisms. In particular, no integral calculation is required at any stage of the analysis .We specifically follow the approach of Topcu et al. [2011], which supplies rigorous upper bounds of integraluncertainties for hierarchies of interconnected subsystems through a systematic computation of moduli ofcontinuity for each individual subsystem. The moduli of continuity supply just the right measure of interac-tion between the subsystems enabling the propagation of uncertainties across the hierarchy. The resultinguncertainty bounds are rigorous, i. e., they are sure to be conservative and result in safe designs; theybecome sharper with an increasing number of input variables (concentration-of-measure phenomenon [Ta-lagrand, 1996]); they do not require differentiability of the subsystem response functions and account forlarge deviations in the response; and only require knowledge of ranges of the input parameters and not theirfull probability distribution, as is the case for Bayesian methods [Dashti and Stuart, 2011]. In addition, thecomputation of the bounds is non-intrusive and can be carried out using existing deterministic models of thesubsystems and external scripts.We specifically aim to assess this hierarchical multiscale approach to uncertainty quantification (UQ) by2igure 1: Example of multiscale hierarchy for modeling strength in metals. Overlaid are examples of modelsor calculations of several of the unit mechanisms, from bottom to top: i) density functional theory (DFT)calculations of the equation of state (EoS) in Ta [Miljacic et al., 2011]; ii) molecular dynamics (MD) cal-culations of kink mobility in Ta [Wang et al., 2001, Segall et al., 2001, Kang et al., 2012]; iii) phase-fieldsimulations of dislocation dynamics and forest hardening [Koslowski et al., 2002]; iv) lamination constructionfor sub-grain dislocation structures [Ortiz and Repetto, 1999]; v) direct numerical simulation (DNS) of poly-crystals [Zhao et al., 2004]; vi) Optimal Transportation Meshfree (OTM) simulations of ballistic perforation[Li et al., 2010].means of an example concerned with the ballistic impact of an elastic-plastic magnesium plate struck by aheavy rigid ball. We adopt as quantity of interest the maximum backface deflection of the plate. We considertwo scales of material response: a macromechanical scale characterized by the Johnson-Cook constitutivemodel [Johnson, 1983]; and a micromechanical scale in which the polycrystalline structure of the materialand its behavior at the single-crystal level are taken into account. For simplicity, we compute polycrystallineaverages by means of Taylor averaging [Taylor, 1938] assuming an isotropic initial texture. The single-crystalplasticity extended to include twinning follows Chang and Kochmann [2015]. We estimate ranges of param-eters for the single-crystal plasticity model based on experimental data compiled from a number of sources.Calculations are carried out using the commercial finite-element package LS-DYNA [Hallquist et al., 2007]on a single converged mesh. The analysis mirrors conventional design testing for impact resistance [Mukaseyet al., 2008], wherein performance is evaluated relative to a targeted set of characterized impact condi-tions. The hierarchical multiscale UQ protocol is implemented using the DAKOTA Version 6 .
12 softwarepackage [Adams et al., 2020] of the Sandia National Laboratories.It bears emphasis that the material model used in the calculations is intended for purposes of demonstrationof the methodology and not as an accurate model of material behavior. This proviso notwithstanding, thecalculations show that the integral uncertainties determined by the hierarchical multiscale UQ approach are3ufficiently tight for use in engineering applications. The analysis also sheds light on the relative contribu-tions of the different unit mechanisms to the integral uncertainty and the dominant propagation paths foruncertainty across the model hierarchy.
We start with a brief review of the uncertainty quantification theory for hierarchical systems [Lucas et al.,2008, Topcu et al., 2011, Sun et al., 2020] that provides the basis for the application to multiscale materialmodelling pursued in this work.
We begin by considering a single-scale, or monolithic, system with uncertain inputs x = ( x , ..., x N ) and aperformance measure y ∈ R . We assume that the inputs are independent, though not necessarily identicallydistributed, random variables defined on a bounded set X ⊂ R N . Let F : X (cid:55)→ R be the response functionthat maps x to y . Define the diameter of F with respect to input x i , 1 ≤ i ≤ N , as d F,i = sup (ˆ x i ,x i ) , (ˆ x i ,x (cid:48) i ) ∈X (cid:12)(cid:12) F (ˆ x i , x i ) − F (ˆ x i , x (cid:48) i ) (cid:12)(cid:12) , (1)where we write ˆ x i = ( x , . . . , x i − , x i +1 , . . . , x N ) . (2)Given a performance threshold y c ∈ R , with failure occurring when y ≥ y c , the McDiarmid’s inequality[Doob, 1940] gives the following upper bound on the probability of failure P [ y ≥ y c ] ≤ exp (cid:18) − M U (cid:19) , (3)where M = max(0 , y c − E [ y ]) is the design margin and U = d F = (cid:18) N (cid:88) i =1 d F,i (cid:19) / (4)is the system uncertainty. Thus, the system is certified if, for a given probability-of-failure tolerance (cid:15) > (cid:18) − M U (cid:19) ≤ (cid:15), (5)or equivalently, CF = MU ≥ (cid:115) log (cid:114) (cid:15) , (6)where CF is a confidence factor in the design [Sharp and Wood-Schultz, 2003].We note that Eq. (3) provides a rigorous probability-of-failure upper bound for the material and, therefore,sets forth a conservative certification criterion. Importantly, no a-priori assumption on the probabilitydistribution of the input variables, or prior , is required. The bound is uniquely determined by the designmargin M and the system uncertainty U , as unambiguously quantified by the system diameter d F , whichare both prior-free. 4igure 2: Graph representation of the input-output relations between the subsystems of a modular system.The nodes of the graph represent the subsystems of the system. An arrow indicates that the outputs of thesubsystem at the beginning of the arrow are among the inputs of the subsystem at the end of the arrow.The subsystems represented by circular boxes do not take input from other systems and are referred to as fundamental . The system has one single root subsystem, i. e., a subsystem that does not feed into any othersubsystem and through which the system output takes place. As noted in the introduction, the macroscopic behavior of many material systems is the result of mechanismsover multiple scales whose functional dependencies define a graph. Conveniently, this graph structure can beexploited to divide the material response into interconnected unit mechanisms, or subsystems, and estimateintegral uncertainties of the entire system from a quantification of uncertainties for each subsystem and anappropriate measure of interaction between the subsystems. We specifically follow the approach of Topcuet al. [2011], which we briefly summarize next.We consider hierarchical or modular systems representable by an oriented graph G ( V, E ) whose nodes V are the subsystems and whose edges E are the interfaces between the subsystems, Fig. 2. Specifically, thegraph contains an oriented edge from b to a if the state of the subsystem a depends on the state of thesubsystem b , i. e., if the outputs of subsystem b are contained among the inputs of system a . We then saythat a is an ancestor of a node b , denoted a ≺ b , and that b is a descendant of a , denoted b (cid:31) a . In order toavoid circular dependencies, we assume that G ( V, E ) is acyclic , i. e., it contains no closed-loop paths. The fundamental subsystems are those that take input from no other subsystems. We assume that the systemcontains a single integral subsystem that does not feed into any other subsystem. Thus, the fundamentalsubsystems and the integral subsystem are the leaves V L and the root R of the graph G ( V, E ), respectively.Suppose that the response of subsystem a ∈ V is characterized by a function F a : X a → Y a that mapsinput parameters X a ∈ X a to outputs Y a ∈ Y a . If a is an ancestor of b in the graph G ( V, E ) then Y b is asubspace of X a . The space of inputs of the system X is the Cartesian product of the input spaces of thefundamental subsystems, i. e., X = (cid:81) a ∈ V L X a . Likewise, the space of outputs of the system is Y = Y R . Forall subsystems other than the fundamental ones, a (cid:54)∈ V L , we have the relation X a = (cid:81) b (cid:31) a Y b , i. e., the inputspace of subsystems i is the Cartesian product of the output spaces of all its descendants. We note thatnon-fundamental subsystems could in principle have inputs of their own, not provided by any descendantsubsystem. We accommodate such cases within the present framework simply by adding a fundamentalsubsystem whose response function is the identity mapping and which supplies the requisite additionalinputs. 5 lgorithm 1 Hierarchical System Evaluation
Require:
Graph G ( V, E ); leaves V L ; root R ; response functions F a : X a → Y a , a ∈ V ; input X ∈ X .i) Initialize: V = V L , k = 0.ii) Reset: V k +1 = { a ∈ V : b ∈ ∪ kl =0 V l , ∀ b (cid:31) a } . for all a ∈ V k +1 do Compute X a = (cid:0) F b ( X b ) , b (cid:31) a (cid:1) . end forif V k +1 = { R } then return Y = F R ( X R ), exit . else k ← k + 1, goto ii). end if The function F : X → Y that describes the response of the integrated system can be evaluated recursively bymeans of Algorithm 1. The algorithm sets forth an “information wave” through the graph that propagatesinformation from the leaves to the root. Thus, in the example of Fig. 2, the sequence of active subsystemsis V = { , , , , , , , , , , , } , V = { , , , } , V = { , } , V = { } and V = { } .Correspondingly, the sequence of subsystem outputs is ( Y , Y , Y , Y ), ( Y , Y ), Y and Y = Y .Let { V k , k = 0 , . . . , N } be the sequence of nodal sets generated during the evaluation of the integratedresponse function F ( X ), with V = V L and V N = { R } . Define X k = (cid:81) a ∈ V k X a and Y k = (cid:81) a ∈ V k +1 X a ,i. e., the combined sets of inputs and outputs for each iteration. We note that X = X , Y k = X k +1 .Without loss of generality we may assume that dim Y N = dim Y = 1. Define further F k : X k → Y k as F k ( X k ) = ( F a ( X a ) , a ∈ V k ), X k ∈ X k , i. e., as the forward map for iteration k . By these definitions andreorganization of the data, we have the composition rule F = F N ◦ · · · ◦ F . (7)Conveniently, the propagation of uncertainty under composition is controlled by the moduli of continuityof the response functions. Recall that, given a function f : R n → R m , a real number δ > A ⊂ R n , the modulus of continuity ω ij ( f, δ, A ) of f i ( x ) with respect to x j over A is defined as [Efimov, 2001,Steffens, 2006] ω ij ( f, δ, A ) = sup {| f i ( x ) − f i ( x (cid:48) ) | : x, x (cid:48) ∈ A, x k = x (cid:48) k for k (cid:54) = j, | x j − x (cid:48) j | ≤ δ } . (8)Thus, ω ij ( f, δ, A ) measures the variation of the function f i ( x ) over A when the variable x j is allowed todeviate by less than δ . We note that this component-wise definition of the modulus of continuity does notrequire the range or image of the function f to be a normed space. This is important in practice, sincethe inputs and outputs of subsystems often comprise variables measured in different units which belong tovector spaces with no natural norm. Consider now two functions f : A ⊂ R n → R m and g : B ⊂ R m → R p ,with B a hyper-rectangle such that f ( A ) ⊂ B , and let δ >
0. Let g ◦ f : A ⊂ R n → R p be the compositionof the f and g , i. e., ( g ◦ f )( x ) = g ( f ( x )). Then, we have [Topcu et al., 2011] ω ij ( g ◦ f, δ, A ) ≤ m (cid:88) k =1 ω ik ( g, ω kj ( f, δ, A ) , B ) . (9)This inequality shows that the moduli of continuity of a composite function g ◦ f can be estimated conser-vatively from the moduli of continuity of the individual functions.6 lgorithm 2 Hierarchical Uncertainty Quantification
Require:
Graph G ( V, E ); leaves V L ; root R ; response functions F a : X a → Y a , a ∈ V ; set A ⊂ X = X . for all j = 1 , . . . , dim X do Compute: D j = sup {| x j − x (cid:48) j | : x, x (cid:48) ∈ A, x k = x (cid:48) k , for k (cid:54) = j } .Compute: D (0) ij = ω ij ( F , D j , A ) , i = 1 , . . . , dim Y . end forfor all k = 1 , . . . , N do Find hyper-rectangles B k ⊂ X k containing ( F k − ◦ · · · ◦ F )( A ).Compute: D ( k ) ij = (cid:80) dim X k l =1 ω il ( F k , D ( k − lj , B k ) , i = 1 , . . . , dim Y k . end forreturn { D F,i = D ( N )1 i , i = 1 , . . . , dim X } .This property of the moduli of continuity in turn enables uncertainties of the integral system to be boundedonce the uncertainties of the subsystems and their interfaces are known. The systematic application ofestimate (9) to the graph G ( V, E ) is described in Algorithm 2 and results in a set of approximate diameters { D F,i , i = 1 , . . . , dim
X } . The fundamental theorem proven by Topcu et al. [2011] is that the approximatediameters { D F,i } bound above the exact integral diameters { d F,i } , Eq. (1), i. e., d F,i ≤ D F,i , i = 1 , . . . , dim X . (10)By the monotonicity of McDiarmid’s inequality, with respect to the diameters, it follows that replacing { d F,i } by { D F,i } in (4) and (5) results in probabilities of integral outcomes and, by extension, in conservativecertification criteria.It bears emphasis that every step of Algorithm 2 requires the execution of subsystem tests only, and thatat no time during the analysis an integral test is required. The sequence D ( k ) ij generated by the algorithmmay be regarded as a measure of uncertainty in the i th output variable due to the variability of the j thinput variable after k levels of operation of the system. The algorithm propagates uncertainties in the inputvariables associated with a leaf i through possibly multiple paths connecting the node i with the root R . Thealgorithm additionally identifies the path responsible for most of the uncertainty of the integral outcomewith respect to the variable i , namely, the path with the highest flow of uncertainty. As just shown in the foregoing, the integral uncertainty for a hierarchical multi-scale system can be evaluatedfrom the moduli of continuity of the sub-system maps without integral testing. We proceed to demonstratethe feasibility of the approach by means of an example of application concerning the ballistic impact of amagnesium plate. For purposes of this demonstration, we restrict attention to three length scales as depictedin Fig. 3: i) the microscale, where the behavior of the magnesium, including slip and twinning, is modelledat the single crystal level; ii) the mesoscale, where the polycrystalline response is computed using Tayloraveraging and the single-crystal model; and iii) the macroscopic impact problem, where the material behavioris approximated by the Johnson-Cook constitutive model and the ballistic performance of the magnesiumplate is simulated using finite elements.For purposes of illustration of the UQ methodology, we assume that all parameters are uncertainty-freesave those shown in red in Fig. 3, namely the micromechanical critical resolved shear stresses (CRSS) forthe basal slip, prismatic slip, pyramidal slip and twinning mechanisms, the mesomechanical Johnson-Cookparameters and the integral quantity of interest. These assumptions lead to the following simplified two-levelhierarchical system: X F −→ Y = X F −→ Y R , (11)7 𝐴, 𝐵, 𝑛, 𝐶)(𝐴, 𝐵, 𝑛, 𝐶) , Figure 3: Graph representation of the hierarchical multiscale model of ballistic impact. Highlighted in redare the micro and mesomechanical parameters assumed to be uncertain in the UQ analysis.cf. Fig. 3. Thus, the graph V of the system contains three nodes, denoted { micro ≡ , meso ≡ , macro ≡ R } ,and X = Uncertain single-crystal model parameters: (CRSS) for the basal slip, prismatic slip, pyramidal slipand twinning mechanisms. F = Mapping that returns an optimal set Y of Johnson-Cook parameters ( A, B, n, C ) for a given realizationof X , obtained by performing regression on stress paths computed using both the Johnson-Cook modeland the single crystal model combined with Taylor averaging along selected strain paths. Y = X = Optimal Johnson-Cook parameters. F = Finite element model of ballistic impact of magnesium plate using the Johnson-Cook model withoptimal parameters X . Y R = Quantity of interest extracted from the results of the finite element calculations.A brief description of the different subsystems is given next, followed by a specification of the maps F and F . We specifically model single-crystal magnesium within the framework of finite-deformation single crystalplasticity extended to include twinning following Chang and Kochmann [2015]. The reader is referred to theoriginal publication for details of the model. The model accounts for basal slip, prismatic slip, pyramidalslip and tensile twinning, cf. Fig. 4. The entire collection of parameters of the model is listed in Table 1,including the values of the parameters determined by Chang and Kochmann [2015]. For each slip system α ,8igure 4: Schematic view of the slip and twin systems of magnesium considered in the present work. Blueplanes represent the slip/twin planes, while red arrows represent the corresponding Burgers vector/twinningshear direction. σ ∞ α is the ultimate strength, h α is the self-hardening modulus, m α is the slip rate sensitivity exponent, ˙ γ ,α isa reference slip rate and h ij are off-diagonals latent hardening moduli. For twinning, h β is the self-hardeningmodulus, m β is the rate-sensitivity exponent, ˙ λ ,β is a reference twin volume-fraction rate, γ t denotes thetwin strain and k ij are interaction moduli. In addition, the elasticity is assumed to be isotropic with Lam´econstants λ e and G . Table 1: Parameters of the single-crystal magnesium model.Parameter Value UnitBasal (cid:104) a (cid:105) h α σ ∞ α h ij m α γ ,α − Prismatic (cid:104) a (cid:105) h α
40 GPa σ ∞ α
170 MPa h ij
20 MPa m α γ ,α − Pyramidal (cid:104) c + a (cid:105) h α
30 GPa σ ∞ α
200 MPa h ij
25 MPa m α γ ,α − Tensile twin h β k ij
40 GPa m β λ ,β − γ t λ e
24 GPa G
25 GPaWe assume that the target plate is made of polycrystalline magnesium and compute the polycrystallineresponse of the polycrystal from the single crystal model just outline by means of Taylor averaging [Taylor,1938]. The computational cost of Taylor averaging is relatively small compared to other averaging schemessuch as periodic boundary conditions [Geers et al., 2010]. In addition, Chang and Kochmann [2015] havedemonstrated good agreement between Taylor averaging and experimental measurements when the numberof the grains is sufficiently large ( > SO (3) [Haar, 1933]. A total of 128 grains areused in the current study for the mesoscale model. 9 .2 Johnson-Cook model The combination of a single-crystal model and Taylor averaging just described, together with suitable initialconditions, enables the calculation of Cauchy stress histories σ ( t ) at a material point from known deformationgradient histories F ( t ). However, such calculations are computationally intensive and cannot be carried outon-the-fly as part of large-scale finite-element calculations. Instead, in said calculations we choose to use asimpler and faster surrogate or mesomechanical model.For definiteness, we specifically choose the Johnson-Cook model [Johnson, 1983] σ M (cid:0) (cid:15) p , ˙ (cid:15) p ) = (cid:2) A + B(cid:15) np (cid:3)(cid:2) C ln ˙ (cid:15) ∗ p (cid:3) , (12)as mesomechanical model, where σ M = (cid:112) / s · s is the Mises stress and s = σ − / tr( σ ) I denotes thedeviatoric part of the Cauchy stress σ , (cid:15) p denotes the equivalent plastic strain, ˙ (cid:15) p is the plastic strain rateand ˙ (cid:15) ∗ p = ˙ (cid:15) p / ˙ (cid:15) p is a normalized plastic strain rate using reference ˙ (cid:15) p . The model parameters are thestrength A , the hardening modulus B , the strain-hardening exponent n , and the rate-sensitivity modulus C ,or X = ( A, B, n, C ). Figure 5: Schematic illustration of the plain-strain cavity expansion problem used to generate representativestress and strain histories for purposes of regression.We now describe the function F that returns the optimal value Y of the mesoscopic model parametersfrom given micromechanical model parameters X . We understand optimality in the sense of achieving theclosest agreement between the mesomechanical and the micromechanical models over selected deformationhistories.In order to generate histories for purposes of regression, we exercise the micromechanical and mesomechan-ical Johnson-Cook models under conditions corresponding to a plane-strain cavity expansion model, Fig. 5.Expanding cavity solutions have been commonly used to approximate conditions that arise in ballistic pen-etration [Bishop et al., 1945, Hanagud and Ross, 1971, Forrestal et al., 1988]. We specifically consider theexpansion of a solid cylinder of radius b at time t = 0 to a hollow cylinder of internal radius a ( t ) = ct at time t at constant rate c . Following Bishop et al. [1945] we assume that the material is incompressible, whereuponthe logarithmic strains in the cylindrical coordinate system ( r, θ ) follow as ε rr ( r, t ) = − ε θθ ( r, t ) = ln r (cid:112) r + a ( t ) and ε zz ( r, t ) = 0 , (13)and the logarithmic strain rate as ˙ ε rr ( r, t ) = − car + a ( t ) . (14)10 straightforward calculation gives the Mises equivalent stress from the Johnson-Cook solution as σ JCeq ( r, t ) = (cid:113) E ln( a ( t )+ r r ) , σ eq ≤ σ ( r, t ) , (cid:16) A + B ( (cid:113) ln( a ( t )+ r r ) (cid:17) n (cid:16) C ln( ca √ ε p ( r + a ( t )) ) (cid:17) , otherwise , (15)where σ ( r, t ) = A (cid:16) C ln( 2 ca √ ε p ( r + a ( t )) ) (cid:17) , (16)is the Mises effective stress at first yield and E is the Young’s modulus. We then measure the discrepancybetween the micromechanical and Johnson-Cook models in the root-mean square (RMS) senseError( Y trial0 |X ) = (cid:16) (cid:88) c (cid:90) b (cid:90) T | σ eq ( r, t ) − σ JCeq ( r, t ) | dt dr (cid:17) / , (17)where the sum is carried out over a representative sample of expansion rates c , σ JCeq ( r, t ) is the Mises equivalentstress computed using Johnson-Cook with parameters Y trial0 and σ eq ( r, t ) is the target Mises equivalent stresscomputed using the micromechanical model with parameters X . The optimal Johnson-Cook parameters Y then follow by minimization of Error( ·|X ), i. e., Y = argmin Error( ·|X ) . (18)This step concludes the regression algorithm and the definition of the sought micro-to-mesomechanical map-ping F : X → Y . r / a H oop s t r ess ( M P a )
25 s(Taylor)25 s(JC)50 s(Taylor)50 s(JC)100 s(Taylor)100 s(JC)
Figure 6: Comparison between the radial distributions of hoop stress σ θθ computed from the micromechanical(Taylor) and the Johnson-Cook (JC) calculations.In calculations, we employ generic algorithms (GA) [Mitchell, 1998] to solve the minimization problem. Theintegrals in (17) are computed by discretizing the spatial and temporal spaces. We specifically consider acylinder radius b = 0 . c = 100 . T = 100 . µ sand set the Young’s modulus to 27 . (cid:15) p = 1 . − . Fig. 6 illustratesthe good agreement that is achieved between the micromechanical and the Johnson-Cook model with A =20 .
98 MPa, B = 161 .
84 MPa, n = 0 .
346 and C = 0 . .
35 and the maximum logarithmic strain rate is 5 , . − , which provide adequatecoverage of the conditions that arise in the ballistic calculations.11 .4 Meso-macromechanical forward solver x x x x x x y R (a) (b) Figure 7: Schematic illustration of magnesium plate struck by a spherical steel projectile at a sub-ballisticspeed. (a) Initial setup. (b) Dynamic indentation process with maximum back surface deflection labeledas Y R . In each subfigure, the top figure shows a perspective view of the projectile/plate system, and thebottom figure shows the cross-sectional view.For purposes of illustration, we consider the problem of assessing the sub-ballistic performance of a magne-sium plate struck by a spherical steel projectile. The problem is solved using the explicit dynamics solveravailable within the commercial software LS-DYNA [Hallquist et al., 2007] and the Johnson-Cook materialmodel. For definiteness, we choose the maximum back surface deflection as the outcome quantity of interest Y R . For every sample Johson-Cook parameter set X the calculations return one value of Y R and thus definea mapping Y R = F ( X ) , (19)which completes the hierarchical model (11).Table 2: Fixed material parameters used in the LS-DYNA simulation. The Gruneisen parameters of themagnesium are provided by Feng et al. [2017].Parameter Value UnitTarget (magnesium) Mass density 1 .
77 g/cm Young’s modulus 27 . .
35 -Specific heat 1 .
04 J/(K · g)Gruneisen intercept 4520 . .
54 -Gruneisen slope 1 .
242 -Projectile (steel) Mass density 7 .
83 g/cm Young’s modulus 210 . .
30 -A schematic of the finite-element model is shown in Fig. 7. The diameter of the projectile is 1 .
12 cm,and the size of the plate is 10 × × .
35 cm. The attack velocity is 200 m/s with normal impact. Thebackface nodes of the target near the edges are fully constrained to prevent displacement in all directions.The projectile is resolved using 864 elements, while the number of elements for the plate is 70 , . µ s before termination. This simulationduration is sufficiently long to allow for the rebound and separation of the projectile from the plate in allthe calculations. The calculations are adiabatic with the initial temperature set at room temperature. Theequation-of-state, which controls the volumetric response of the material, is assumed to be of the Gruneisentype. For simplicity, the projectile is assumed to be rigid and uncertainty-free. All other material parametersare fixed and listed in Table 2. 12 .5 Implementation of Hierarchical UQ Table 3: Experimentally reported CRSS at room temperature for slip and twin systems in pure magnesiumsingle crystals. The unit is MPa. CRSS SourceBasal (cid:104) a (cid:105) [0 .
44 0 .
58] Yoshinaga and Horiuchi [1964]0 .
48 Akhtar and Teghtsoonian [1969]0 .
52 Herring and Meshii [1973]0 .
76 Burke and Hibbard [1952]0 .
81 Schmid [1931]4 .
16 Chapuis and Driver [2011]Prismatic (cid:104) a (cid:105) . . . (cid:104) c + a (cid:105) [27 .
76 49 .
22] Obara et al. [1973]44 . .
87 58 .
52] Ando et al. [2007][55 .
48 86 .
0] Kitahara et al. [2007]Tensile twin 1 .
86 Roberts [1960][3 . , .
7] Chapuis and Driver [2011]Table 4: Lower and upper bounds for the CRSS in different slip and twin systems. The unit is MPa.Basal (cid:104) a (cid:105) Prismatic (cid:104) a (cid:105) Pyramidal (cid:104) c + a (cid:105) Tensile twin[0 .
44 4 .
16] [17 . .
0] [27 .
76 86 .
0] [1 .
86 11 . (cid:104) a (cid:105) , prismatic (cid:104) a (cid:105) , pyramidal (cid:104) c + a (cid:105) and tensile twin cf. Fig. 4. These CRSSs are allowedto vary over a certain range in order to cover the experimental data. Table 3 lists a compilation of CRSSvalues reported in the literature for the slip and twin systems in pure magnesium single crystals at roomtemperature. Based on these data, the ranges of CRSSs used as the input in the multiscale UQ analysisare listed in Table 4. All other material parameters in the crystal-plasticity model are fixed and listed inTable 1.The calculation of the diameters D ( k ) ij requires the solution of a global, constrained optimization prob-lem, while the calculation of hyper-rectangles B k entails the solutions of two global optimization problemswith regard to each output. In order to solve these optimization problems efficiently, we have developed anon-intrusive, high-performance computational framework based on DAKOTA Version 6 .
12 software pack-age [Adams et al., 2020] of the Sandia National Laboratories. We employ genetic algorithms (GA) to solveall the optimization problems involved in the UQ analysis [Mitchell, 1998, Sun et al., 2020]. GA, as a globaland derivative-free optimization method, offers great flexibility in applications, such as considered here, tohighly non-linear non-convex problems without the availability of gradients. Another advantage of GA is itshigh degree of concurrency, since individuals in each iteration can be evaluated independently across multipleprocessors. In all the GA calculations, we choose throughout a fixed population size of 64. The crossoverrate and mutation rate are fixed at 0 . . We recall that the modular UQ analysis is a divide-and-conquer approach that enables each subsystem to beassessed independently. In keeping with this paradigm, we analyze the micro-meso and meso-macro maps13ndependently and combine the results to determine integral uncertainty bounds for the entire system.
Basal Pris. Pyra. Twin
CRSS M odu l u s o f c on t i nu i t y ( M P a ) (a) Basal Pris. Pyra. Twin
CRSS M odu l u s o f c on t i nu i t y ( M P a ) (b) Basal Pris. Pyra. Twin
CRSS M odu l u s o f c on t i nu i t y (c) Basal Pris. Pyra. Twin
CRSS M odu l u s o f c on t i nu i t y (d) Figure 8: Moduli of continuity from single-crystal properties to Johnson-Cook (JC) parameters. (a) JCparameter A . (b) JC parameter B . (c) JC parameter n . (d) JC parameter C .The moduli of continuity for the micro-meso mapping, relating uncertainties in the basal, prismatic, pyra-midal and twin CRSS to uncertainties in the Johnson-Cook parameters are shown in Fig. 8. Clearly themicro-mechanical parameters contribute to different degrees to the uncertainties in the Johnson-Cook param-eters. Remarkably, twinning uncertainty contributes the least. The parameters A and n are most sensitiveto the pyramidal CRSS, while C is most sensitive to the prismatic CRSS. By contrast, the parameter B isequally sensitive to all three slip mechanisms.The integral moduli of continuity relating uncertainties in the CRSSs to the maximum back surface deflectionof the plate are shown in Fig. 9. An important property of the moduli of continuity is that, since they aredimensionally homogeneous, they can be compared and rank-ordered, which in turn provides a quantitativemetric of the relative contributions of the input parameters to the overall uncertainty. In the present case,the ranking is C > B > A > n for basal CRSS,
C > A > n > B for both prismatic and pyramidal CRSSs,and
C > n > B > A for twin CRSS. Remarkably, the slip uncertainty flow through the Johnson-Cookparameter C is significantly larger than other paths, whereas the twin uncertainty is transmitted nearlyuniformly by all Johnson-Cook parameters.Numerical results for the modular upper bounds are collected in Fig. 10. The rank-ordering of the CRSSs tothe overall uncertainty in ballistic performance is found to be pyramidal > prismatic > basal > twin, with14 asal Pris. Pyra. Twin CRSS M odu l u s o f c on t i nu i t y ( c m ) A B n C
Figure 9: Moduli of continuity from single crystal properties to ballistic performance through differentJohnson-Cook parameters.
Basal Pris. Pyra. Twin
CRSS M odu l a r upp e r bound ( c m ) Figure 10: Sub-diameter upper bounds computed using modular approach.the pyramidal and prismatic CRSSs contributing the most, the twin CRSS the least and the basal CRSS inbetween. We also note that the modular upper bounds lie above the system diameters, cf. Eq. (10). Thus,the modular upper bounds provide a conservative estimate of the probability of departure from the meanwhen inserted into McDiarmid’s inequality Eq. (3) in place of the system diameter.15
Concluding remarks
We have presented a framework to assess the uncertainties of a hierarchical multi-scale material system.The hierarchical structure of the multi-scale systems can be viewed as a directed graph, and we exploit thisstructure by bounding the uncertainty at each scale and then combining the partial uncertainties in a waythat provides a bound on the overall or integral uncertainty. The bound provides a conservative estimate onthe integral uncertainty. Importantly, this approach does not require integral calculations, which often areprohibitively expensive.We have demonstrated the framework on the problem of ballistic impact of a polycrystalline magnesiumplate. Magnesium and its alloys are of current interest as promising light-weight structural and protectivematerials. We start at the microscopic sub-grain scale where the behavior of various slip and twinningsystems is described using extended crystal plasticity, pass to the mesoscopic scale of a representative volumeinvolving multiple grains where the behavior is described using a Johnson-Cook constitutive model and studya macroscopic problem of ballistic impact of a plate. We study the uncertainty of the ballistic response dueto uncertainties in the strength of individual slip and twin system.We close with a discussion of how our approach also informs the ‘materials-by-design’ approach. In thisapproach, we seek to ‘design’ a material with desired properties by affecting some underlying mechanismat a fine scale. As a concrete example, consider the improvement of ballistic performance (measured bymaximum deflection) of a magnesium plate. A potential way to doing so is to change the CRSS of a slip ortwin system by adding solutes or precipitates. We can estimate the change in CRSS of a particular systemto a solute by conducting molecular dynamics simulations. The question then, is how does this change inCRSS manifest itself in a change of ballistic performance. In the notation of (11), we seek to understandhow the change of X affects Y R , or the modulus of continuity of this composite map. The method we havepresented provides an outer bound (optimistic in the context of design) on this modulus. Importantly we cancompute this bound by studying individual scales ( F and F ) without the need for prohibitively expensiveintegral calculations. Acknowledgement
We are grateful to Dennis Kochmann for making available to us the single-crystaland Taylor averaging codes used in this work to generate data. This research was sponsored by the ArmyResearch Laboratory and was accomplished under Cooperative Agreement Number W911NF-12-2-0022. Theviews and conclusions contained in this document are those of the authors and should not be interpretedas representing the official policies, either expressed or implied, of the Army Research Laboratory or theU.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmentpurposes notwithstanding any copyright notation herein.
References
B. M. Adams, W. J. Bohnhoff, K. R. Dalbey, M. S. Ebeida, J. P. Eddy, M. S. Eldred, R. W. Hooper,P. D. Hough, K. T. Hu, J. D. Jakeman, M. Khalil, K. A. Maupin, J. A. Monschke, E. M. Ridgway,A. A. Rushdi, D. T. Seidl, J. A. Stephens, L. P. Swiler, and J. G. Winokur. Dakota, a multilevel parallelobject-oriented framework for design optimization, parameter estimation, uncertainty quantification, andsensitivity analysis: version 6.12 user’s manual.
Sandia National Laboratories, Tech. Rep. SAND2020-5001 ,2020.A. Akhtar and E. Teghtsoonian. Solid solution strengthening of magnesium single crystals—i alloying be-haviour in basal slip.
Acta Metallurgica , 17(11):1339–1349, 1969.S. Ando, N. Harada, M. Tsushida, H. Kitahara, and H. Tonda. Temperature dependence of deformationbehavior in magnesium and magnesium alloy single crystals. In
Key Engineering Materials , volume 345,pages 101–104. Trans Tech Publications Ltd, 2007.16. Bhattacharya.
Microstructure of martensite: why it forms and how it gives rise to the shape-memoryeffect , volume 2. Oxford University Press, 2003.R. Bishop, R. Hill, and N. Mott. The theory of indentation and hardness tests.
Proceedings of the PhysicalSociety , 57(3):147, 1945.E. Burke and W. Hibbard. Plastic deformation of magnesium single crystals.
Jom , 4(3):295–303, 1952.C. M. Byer, B. Li, B. Cao, and K. Ramesh. Microcompression of single-crystal magnesium.
Scripta Materi-alia , 62(8):536–539, 2010.Y. Chang and D. M. Kochmann. A variational constitutive model for slip-twinning interactions in hcp metals:application to single-and polycrystalline magnesium.
International Journal of Plasticity , 73:39–61, 2015.A. Chapuis and J. H. Driver. Temperature dependency of slip and twinning in plane strain compressedmagnesium single crystals.
Acta Materialia , 59(5):1986–1994, 2011.M. Dashti and A. M. Stuart. Uncertainty quantification and weak approximation of an elliptic inverseproblem.
SIAM Journal on Numerical Analysis , 49(6):2524–2542, 2011.J. L. Doob. Regularity properties of certain families of chance variables.
Transactions of the AmericanMathematical Society , 47(3):455–486, 1940.A. Efimov. Modulus of continuity, encyclopaedia of mathematics, 2001.J. Feng, P. Chen, Q. Zhou, K. Dai, E. An, and Y. Yuan. Numerical simulation of explosive welding usingsmoothed particle hydrodynamics method.
International Journal of Multiphysics , 11(3), 2017.P. W. Flynn, J. Mote, and J. E. Dorn. On the thermally activated mechanism of prismatic slip in magnesiumsingle crystals.
Transactions of the Metallurgical Society of AIME , 221(6):1148–1154, 1961.M. J. Forrestal, K. Okajima, and V. K. Luk. Penetration of 6061–t651 aluminum targets with rigid longrods.
Journal of Applied Mechanics , 55(4):755–760, 1988.M. G. Geers, V. G. Kouznetsova, and W. Brekelmans. Multi-scale computational homogenization: Trendsand challenges.
Journal of computational and applied mathematics , 234(7):2175–2182, 2010.M. S. Green. Markoff random processes and the statistical mechanics of time-dependent phenomena. ii.irreversible processes in fluids.
The Journal of Chemical Physics , 22(3):398–413, 1954.A. Haar. Der massbegriff in der theorie der kontinuierlichen gruppen.
Annals of mathematics , pages 147–169,1933.E. Hall. The deformation and ageing of mild steel: Iii discussion of results.
Proceedings of the PhysicalSociety. Section B , 64(9):747, 1951.J. O. Hallquist et al. Ls-dyna keyword users manual.
Livermore Software Technology Corporation , 970:299–800, 2007.S. Hanagud and B. Ross. Large deformation, deep penetration theory for a compressible strain-hardeningtarget material.
AIAA Journal , 9:905–911, 1971.R. Herring and M. Meshii. Thermally activated deformation of gold single crystals.
Metallurgical Transac-tions , 4(9):2109–2114, 1973.J. Hutchinson. Elastic-plastic behaviour of polycrystalline metals and composites.
Proceedings of the RoyalSociety of London. A. Mathematical and Physical Sciences , 319(1537):247–272, 1970.G. R. Johnson. A constitutive model and data for materials subjected to large strains, high strain rates, andhigh temperatures.
Proc. 7th Inf. Sympo. Ballistics , pages 541–547, 1983.17. Kang, V. V. Bulatov, and W. Cai. Singular orientations and faceted motion of dislocations in body-centered cubic crystals.
Proceedings of the National Academy of Sciences , 109(38):15174–15178, 2012.T. Kitahara, S. Ando, M. Tsushida, H. Kitahara, and H. Tonda. Deformation behavior of magnesium singlecrystals in c-axis compression. In
Key Engineering Materials , volume 345, pages 129–132. Trans TechPubl, 2007.C. Kohler, P. Kizler, and S. Schmauder. Atomistic simulation of precipitation hardening in α -iron: influ-ence of precipitate shape and chemical composition. Modelling and Simulation in Materials Science andEngineering , 13(1):35, 2004.M. Koslowski, A. Cuitino, and M. Ortiz. A phase-field theory of dislocation dynamics, strain hardening andhysteresis in ductile single crystals.
Journal of the Mechanics and Physics of Solids , 50:2597–2635, 2002.R. Kubo. Statistical-mechanical theory of irreversible processes. i. general theory and simple applications tomagnetic and conduction problems.
Journal of the Physical Society of Japan , 12(6):570–586, 1957.G. Leibfried and N. Breuer.
Point defects in metals I: introduction to the theory , volume 81. Springer, 2006.B. Li, F. Habbal, and M. Ortiz. Optimal transportation meshfree approximation schemes for fluid and plasticflows.
International Journal for Numerical Methods in Engineering , 83(12):1541–1579, 2010.L. J. Lucas, H. Owhadi, and M. Ortiz. Rigorous verification, validation, uncertainty quantification andcertification through concentration-of-measure inequalities.
Computer Methods in Applied Mechanics andEngineering , 197(51-52):4591–4609, 2008.U. Messerschmidt.
Dislocation dynamics during plastic deformation , volume 129. Springer Science & BusinessMedia, 2010.L. Miljacic, S. Demers, and A. van de Walle. Equation of state and viscosity of tantalum and iron from firstprinciples.
American Physical Society , APS March Meeting 2011:21–25, March 2011.M. Mitchell.
An introduction to genetic algorithms . MIT press, 1998.M. Mukasey, J. L. Sedgwick, and D. Hagy. Ballistic resistance of body armor, nij standard-0101.06. , 2008.T. Obara, H. Yoshinga, and S. Morozumi. { }(cid:104) (cid:105) slip system in magnesium. Acta Metallurgica , 21(7):845–853, 1973.G. B. Olson. Designing a new material world.
Science , 288(5468):993–998, 2000.M. Ortiz and E. Repetto. Non-convex energy minimization and dislocation structures in ductile singlecrystals.
Journal of the Mechanics and Physics of Solids , 47:397–462, 1999.M. Ortiz, A. M. Cuiti˜no, J. Knap, and M. Koslowski. Mixed atomistic-continuum models of materialbehavior: The art of transcending atomistics and informing continua.
MRS Bulletin , 26(3):216–221, 2001.N. Petch. The cleavage strength of polycrystals.
Journal of the Iron and Steel Institute , 174:25–28, 1953.R. Phillips.
Crystals, defects and microstructures: Modeling across scales . Cambridge University Press, 2001.R. E. Reed-Hill and W. D. Robertson. Deformation of magnesium single crystals by nonbasal slip.
Jom , 9(4):496–502, 1957.C. S. Roberts.
Magnesium and its Alloys . Wiley, 1960.E. Schmid. Beitr¨age zur physik und metallographie des magnesiums.
Zeitschrift f¨ur Elektrochemie undangewandte physikalische Chemie , 37(8-9):447–459, 1931.D. E. Segall, T. A. Arias, A. Strachan, and W. A. Goddard. Accurate calculations of the peierls stress insmall periodic cells.
Journal of Computer-Aided Materials Design , 8(2):161–172, 2001.18. H. Sharp and M. M. Wood-Schultz. Qmu and nuclear weapons certification-what’s under the hood?
LosAlamos Science , 28:47–53, 2003.K.-G. Steffens. Constructive function theory: Kharkiv.
The History of Approximation Theory: From Eulerto Bernstein , pages 167–190, 2006.X. Sun, T. Kirchdoerfer, and M. Ortiz. Rigorous uncertainty quantification and design with uncertainmaterial models.
International Journal of Impact Engineering , 136:103418, 2020.M. Talagrand. A new look at independence.
The Annals of probability , pages 1–34, 1996.G. Taylor. Plastic strain in metals.
Journal of the Institute of Metals , 62:307–324, 1938.U. Topcu, L. J. Lucas, H. Owhadi, and M. Ortiz. Rigorous uncertainty quantification without integraltesting.
Reliability Engineering & System Safety , 96(9):1085–1091, 2011.G. Wang, A. Strachan, T. Cagin, and W. A. Goddard. Molecular dynamics simulations of 1/2 a¡1 1 1¿ screwdislocation in ta.
Materials Science and Engineering: A , 309-310:133 – 137, 2001. Dislocations 2000: AnInternational Conference on the Fundamentals of Plastic Deformation.Y. Wei and L. Anand. Grain-boundary sliding and separation in polycrystalline metals: application tonanocrystalline fcc metals.
Journal of the Mechanics and Physics of Solids , 52(11):2587–2616, 2004.H. Yoshinaga and R. Horiuchi. On the nonbasal slip in magnesium crystals.
Transactions of the JapanInstitute of Metals , 5(1):14–21, 1964.Z. Zhao, R. Radovitzky, and A. Cuitino. A study of surface roughening in fcc metals using direct numericalsimulation.