A Computationally Tractable Framework for Nonlinear Dynamic Multiscale Modeling of Membrane Fabric
Philip Avery, Daniel Z. Huang, Wanli He, Johanna Ehlers, Armen Derkevorkian, Charbel Farhat
RReceived: Added at production Revised: Added at production Accepted: Added at productionDOI: xxx/xxxx
RESEARCH ARTICLE
A computationally tractable framework for nonlinear dynamicmultiscale modeling of membrane fabric
Philip Avery | Daniel Z. Huang | Wanli He | Johanna Ehlers | Armen Derkevorkian | CharbelFarhat Department of Aeronautics andAstronautics, Stanford University, Stanford,CA 94305, USA Department of Mechanical Engineering,Stanford University, Stanford, CA 94305,USA Institute for Computational andMathematical Engineering, StanfordUniversity, Stanford, CA 94305, USA Jet Propulsion Laboratory, CaliforniaInstitute of Technology, Pasadena, CA91109, USA
Summary
A general-purpose computational homogenization framework is proposed for thenonlinear dynamic analysis of membranes exhibiting complex microscale and/ormesoscale heterogeneity characterized by in-plane periodicity that cannot be effec-tively treated by a conventional method, such as woven fabrics. The proposedframework is a generalization of the “Finite Element squared” (or FE ) method inwhich a localized portion of the periodic subscale structure – typically referred toas a Representative Volume Element (RVE) – is modeled using finite elements. Thenumerical solution of displacement-driven problems using this model furnishes amapping between the deformation gradient and the first Piola-Kirchhoff stress ten-sor. This unconventional material model can be readily applied in the context ofmembranes by using a variant of the approach proposed by Klinkel and Govindjee for using conventional, finite strain, three-dimensional material models in beam andshell elements. The approach involves the numerical enforcement of the plane stressconstraint, which is typically performed on out-of-plane components of the secondPiola-Kirchhoff stress tensor ( 𝑺 ). Observing the principal of frame invariance, theRVE solution is reinterpreted as a mapping between the right stretch tensor and thesymmetric Biot stress tensor conjugate pair. This facilitates the development of adrop-in replacement for any conventional finite strain plane stress material model for-mulated in terms of the in-plane components of the Green-Lagrange strain tensor and 𝑺 . Finally, computational tractability is achieved by introducing a regression-basedsurrogate model to avoid further solution of the RVE model when data sufficientto fit a model capable of delivering adequate approximations is available. For thispurpose, a physics-inspired training regimen involving the utilization of our gener-alized FE method to simulate a variety of numerical experiments – including butnot limited to uniaxial, biaxial and shear straining of a material coupon – is proposedas a practical method for data collection. The proposed framework is demonstratedfor a Mars landing application involving the supersonic inflation of an atmosphericaerodynamic decelerator system that includes a parachute canopy made of a wovenfabric. Several alternative surrogate models are evaluated including a neural network. KEYWORDS:
FE, regression, multiscale, woven fabrics, neural network a r X i v : . [ c s . C E ] J u l AVERY
ET AL
Nonlinear multiscale problems – defined here as those nonlinear problems that exhibit vastly different scale features that aresignificant to the macroscopic behavior – are ubiquitous in science and engineering. They arise, for example, in the modeling ofwoven fabrics (see Figure 1) used in body armor and inflatable structures such as vehicle air bags, parachutes and other atmo-spheric decelerators. Numerical methods that attempt to resolve all relevant scales can lead to massive discretized problems.However, recent developments using a variety of alternative surrogate modeling techniques – including nonlinear model orderreduction , kriging , and neural networks – to accelerate the solution of one or more scales within the context of a computa-tional homogenization framework present a coherent methodology by which a computationally tractable approximation can beattained without resorting to ad-hoc approximations. Notably, thin shell and membrane discretizations have not been consideredin this context prior to this work although several frameworks for multiscale modeling of shells without emphasis on computa-tional efficiency have been proposed . In particular we will address the case of a hybrid discretization in which plane stressmembrane elements are employed for the sake of convenience and numerical efficiency at the macroscopic scale, while three-dimensional solid elements are preferred for the sake of generality and in order to most precisely represent geometric featuresand deformation modes at the mesoscopic and/or microscopic scales. FIGURE 1
Optical microscope imaging of parachute membrane fabric yarns under tensile loading.We choose to build our framework on a macro-micro concept that generalizes to 𝑛 -levels, although without loss of generalitywe present only the two-level case here anticipating that this case can be sufficient for many problems of interest. The frameworkallows for the treatment of unilateral contact constraints at both macroscale and microscale. We propose using a microscalemodel discretized with solid elements to allow accurate representation of microscopic geometric features such as yarns andvoids. However, the proposed framework readily generalizes to alternative microscale discretizations such as shell elements. Forthe case of a macroscale model that is also discretized with solid elements, well-established localization/homogenization scalebridging strategies have been developed and typically provide a mapping between the three-dimensional deformation gradient( 𝑭 ) and first Piola-Kirchhoff stress ( 𝑷 ) tensors from which a constitutive relation is inferred. However, when the macroscalemodel is discretized with membrane elements and the microscale is discretized with solid elements then coupling between thetwo scales requires careful attention. The proposed treatment has two novel components. • First, we observe that due to the principal of material frame invariance the conventional three-dimensional 𝑭 – 𝑷 scalebridging approach can be reformulated using the polar decomposition of the deformation gradient to furnish a mapping VERY
ET AL between the right stretch tensor and its conjugate, the symmetric Biot stress tensor. Using this reformulation combinedwith some straightforward transformations, we will show that an unconventional FE material model can conveniently beused as a drop-in replacement for any conventional three-dimensional material model formulated as a mapping between thecommonly used Green-Lagrange strain ( 𝑬 ) and its conjugate, the second Piola-Kirchhoff stress tensor ( 𝑺 ). The relevanceof this development to the issue of multiscale membrane-solid coupling will be addressed in what follows. • Material models used in membrane elements are typically of the plane stress variety. In some cases, a plane stress variantof a three-dimensional material model can be derived for which the plane stress condition is enforced analytically. When ananalytical solution is not available, numerical enforcement of the plane stress condition is commonly used, for example inthe case of J2 elastoplasticity . This involves solving numerically – with a root-finding method such as Newton’s methodor the bisection method – a nonlinear equation to enforce the plane stress condition. Klinkel and Govindjee have shown how numerical enforcement of the plane stress condition can be used to construct an interface that in principle enablesany three-dimensional material model to be “converted” into a plane stress variant which can then be used in a shell orbeam element. However, the plane stress condition is typically expressed by constraints on the out-of-plane componentsof the second Piola-Kirchhoff stress tensor, so the method of Klinkel and Govindjee is presented by the authors in themost convenient way using the 𝑬 – 𝑺 conjugate pair. We simply observe that the method of Klinkel and Govindjee can betrivially adapted to membranes and furthermore can be conveniently used with any material model of the form 𝑺 = ̂ 𝑺 ( 𝑬 ) including, but not limited to, constitutive relations inferred from a FE computational homogenization formulated usingthe polar decomposition as described previously.Motivated by the fact that the proposed plane stress constitutive law is essentially a mapping is between two pairs of 3-dimensional vectors, we consider a lightweight alternative in which a regression based-model is used as a surrogate forconstitutive function evaluations that would otherwise require solution of a finite-element model of the microscale RVE. Threealternative surrogates each capable of achieving computational tractability are presented and evaluated: (1) the classical linearelastic model fitted to data using linear regression, (2) a quadratic model fitted to data using linear regression, and (3) an artificialneural network model fitted to data using the PyTorch library. In each case the data used to train and test the model is obtainedby exercising the proposed high-fidelity multiscale membrane model on a series of numerical experiments intended to mimicthe familiar physical experimental-based methodology typically used in the development of conventional material models.The remainder of this paper is organized as follows. Section 2 provides an overview of the proposed two-level multiscaleframework with a locally attached microscale, focusing on the context of large-deformation structural mechanics with macro-scopic discretization using membrane elements, microscale discretization using 3D solid elements, bridging between the scales,and the solution of the discrete coupled multiscale problem including treatment of contact at both scales. In Section 3, threeregression-based surrogate microscale models and their training methodologies are presented and compared. A numericalexample is provided in Section 4 to evaluate the proposed framework, involving a realistic simulation of the deployment of aDisk-Gap-Band (DGB) parachute in the Martian atmosphere. Finally, we offer conclusions in Section 5. In this section, a multiscale continuum mechanics formulation suitable for membranes based on the concept of a locally attachedmicrostructure is presented. As formulated, the stress-strain relationship for a heterogeneous membrane is not defined by aconventional plane stress constitutive law but rather by (a) the solution at each material point of one or more boundary valueproblems governing its microstructure, and (b) the numerical enforcement of the plane stress condition. Although the conceptgeneralizes naturally to three or more scales, for simplicity it is applied here to address problems that exhibit precisely twoseparate scales. Specifically, the stress-strain relationship at the coarse scale is defined by the solution of the boundary valueproblems at the fine scale, an appropriate scale transition method, and a constraint on the out-of-plane components of thehomogenized stress tensor. At the finest scale, where all heterogeneities can be adequately resolved and described by an availableconstitutive theory, it is defined by an analytical constitutive law. All considered length scales are assumed to be much largerthan the molecular dimension so that the continuum assumption holds. Furthermore, scale separation is assumed to looselycouple the various scales through localization from coarse to fine scales and homogenization from fine to coarse scales. Forfurther details, the reader is referred to for the concept of a locally attached microstructure. The approach adopted here
AVERY
ET AL can be interpreted as a kind of generalization and/or application to the case of membranes of the localization/homogenizationscale bridging strategy presented in . Consider a domain ℬ ⊂ ℝ that defines a highly heterogeneous membrane structure of interest. Assume that its boundary 𝜕 ℬ is subject to prescribed displacements on 𝜕 ℬ 𝒖 ⊂ 𝜕 ℬ , and tractions on 𝜕 ℬ ⧵ 𝜕 ℬ 𝒖 . Let 𝝋 𝑡 ∶ ℬ → ℬ 𝑡 denote the nonlineartransformation that maps a point in the reference configuration, 𝑿 ∈ ℬ , at time 𝑡 to a counterpart in the current configuration, 𝒙 ( 𝑿 ; 𝑡 ) = 𝝋 𝑡 ( 𝑿 ) ∈ ℬ 𝑡 . The current configuration of the membrane ℬ 𝑡 is assumed to be defined as ℬ 𝑡 = { 𝒙 ∈ ℝ | 𝒙 = 𝜙 ( 𝜉 (1) , 𝜉 (2) ) + 𝜉 (3) 𝒏 } where the map 𝜙 ∶ 𝒜 → ℝ defines the current position of the mid-surface of the membrane, ( 𝜉 (1) , 𝜉 (2) ) ∈ 𝒜 ⊂ ℝ arecoordinates parameterizing the mid-surface, 𝒏 ∈ ℝ is the unit normal to the mid-surface in the current configuration, 𝜉 (3) ∈ [ − ℎ ∕2 , ℎ ∕2 ] is a coordinate parameterizing the direction normal to the surface, and ℎ is the upper bound of the membranethickness. Similarly, the reference configuration of the membrane ℬ is defined as ℬ = { 𝑿 ∈ ℝ | 𝑿 = Φ ( 𝜉 (1) , 𝜉 (2) ) + 𝜉 (3) 𝑵 } where the the map Φ ∶ 𝒜 → ℝ defines the reference position of the mid-surface of the membrane, and 𝑵 ∈ ℝ is theunit normal to the mid-surface in the reference configuration. The deformation of this domain is governed by a reduction of thefinite deformation continuum equations to the mid-surface with a plane stress but otherwise unknown constitutive law due to theassumed highly heterogeneous fine scale structure. For this reason, generalizing the work described in , the deformationproblem is solved here by locally attaching an appropriately defined microstructure to each mid-surface point, computing thestress-strain relationship at each such point through the solution of a microstructure boundary value problem, bridging thescales via a localization and homogenization strategy, and numerically enforcing the plane stress constraint on the resultinghomogenized stress tensor. An appropriately defined microstructure in this context is one that represents only a minuscule“Representative Surface Element” (RSE) of the membrane within which the entire thickness of the membrane is represented inits entirety. Hence, the range of the in-plane coordinates 𝜉 (1) , 𝜉 (2) in the microscale domain should be much smaller than in themacroscale domain, while the ranges of the normal coordinate 𝜉 (3) should be identical in both domains. The separation of scalesand assumed periodicity in only two of the three spatial dimensions is a notable characteristic of the problem of interest and itsproposed treatment that distinguishes it from the ubiquitous alternative multiscale treatments devised for fully three-dimensionalscale bridging.Here and throughout the remainder of this paper, the subscripts and denote quantities associated with the coarse ( -th)and fine ( -st) scales, respectively. For simplicity, a 𝑘 -th scale is also referred to as scale 𝑘 or level 𝑘 , interchangeably. Thedeformation at both scales is governed by the finite deformation continuum equations, with the stress-strain relationship definedby the solution of a constrained boundary value problem formulated at a finer scale for level 𝑘 = 0 , or an assumed constitutivelaw at the fine scale designated by level 𝑘 = 1 . Let 𝝋 ∶ ℬ → ℬ ′1 denote the nonlinear transformation that maps a point in thefine scale reference configuration, 𝑿 ∈ ℬ , to a counterpart in the fine scale current configuration, 𝒙 ( 𝑿 ) = 𝝋 ( 𝑿 ) ∈ ℬ ′1 . Asin the formulation of the macroscale problem, we define 𝜕 ℬ as the boundary of ℬ and 𝜕 ℬ 𝒖 as its part where a displacementis prescribed.The boundary conditions at scale are defined by the physical problem of interest, while those at scale depend on thedeformations at the coarse scale. The constitutive law at the fine scale is chosen based on the expected response of this scale, whileat the coarse scale there is no preassigned constitutive law but rather a dependence on the microstructure response to evaluatethe constitutive function. Arbitrarily complex fine scale constitutive relationships involving nonlinearities and path-dependencyare allowed and described here by 𝑺 = ̂ 𝑺 ( 𝑬 , 𝚵 ) where 𝑺 , 𝑬 and 𝚵 denote the microscale second Piola-Kirchhoff stress tensor, Green-Lagrange strain tensor and historyvariables, respectively, and ̂ 𝑺 is the microscale constitutive function. At the coarse scale we wish to devise a multiscale, planestress constitutive function of the form 𝑺 𝑚 = ̂ 𝑺 𝑚 ( 𝑬 𝑚 , 𝚵 ) (1)where 𝑺 , 𝑬 and 𝚵 denote the macroscale second Piola-Kirchhoff stress tensor, Green-Lagrange strain tensor and historyvariables, respectively. The subscript 𝑚 applied to a tensor quantity (for example 𝑺 𝑚 and 𝑬 𝑚 ) denotes a restriction of the tensor VERY
ET AL to its in-plane membrane components. For example, the membrane part of 𝑺 is given by 𝑺 𝑚 = [ 𝑆 (11)0 𝑆 (12)0 𝑆 (21)0 𝑆 (22)0 ] The superscript 𝑚 applied to a constitutive function (for example ̂ 𝑺 𝑚 ) indicates that the function is a particular plane stress typeof constitutive relation that evaluates the in-plane membrane components of a stress tensor while constraining the out-of-planecomponents to be zero. A general numerical procedure used to construct this function will be defined subsequently. Following the work presented in , the boundary conditions on ℬ are defined so that the pointwise deformation gradient tensorat level , 𝑭 , is equal to the volumetric average of the deformation gradient tensor at level 𝑭 = 1 | ℬ | ∫ ℬ 𝑭 𝑑𝑉 This localization transmission condition can be conveniently enforced by prescribing a boundary deformation of the form 𝒙 || 𝜕 ℬ 𝒖 = 𝑿 || 𝜕 ℬ 𝒖 𝑭 + 𝒘 (2)subject to some conditions (see ), where 𝒘 represents the non-uniform part of the boundary deformation. Without loss ofgenerality, in this work we assume uniform essential boundary conditions defined by the condition 𝒘 = 0 .The pointwise first Piola-Kirchhoff stress tensor at level is defined as the volumetric average of the stress tensor at level 𝑷 = 1 | ℬ | ∫ ℬ 𝑷 𝑑𝑉 (3)This homogenization transmission condition can be conveniently determined from quantities defined solely on 𝜕 ℬ 𝒖 by applyinga Gauss-type identity to (3) 𝑷 = 1 | ℬ | ∫ 𝜕 ℬ 𝒖 𝑷 𝑵 ⊗ 𝑿 𝑑𝐴 (4a) = 1 | ℬ | 𝑿 𝑇 ||| 𝜕 ℬ 𝒖 𝒇 || 𝜕 ℬ 𝒖 (4b)where 𝒇 || 𝜕 ℬ 𝒖 is the vector of so-called reaction forces associated with the prescribed deformations (2) and the superscript 𝑇 designates the transpose operation.In this context, the microscale volume measure | ℬ | should be interpreted as the entire volume of a bounding box enclosingthe microscale volume (see Figure 2) including both regions of solid material and voids. The height of the bounding box ℎ + 𝜀 should be slightly larger than the minimum enclosing dimension ℎ in the 𝜉 (3) direction (i.e. 𝜀 > ) so that the microscale volume ℬ does not intersect with the box’s upper and lower faces. The magnitude of 𝜀 is otherwise arbitrary as the dependence ofthe homogenized stress tensor on this parameter will be subsequently canceled when evaluating the membrane stress resultant.Note that the boundary 𝜕 ℬ 𝒖 used to define the transmission conditions are entirely contained within the four side faces of thebounding box, i.e. the faces whose normals coincide with the 𝜉 (1) and 𝜉 (2) axes.Equations (2) and (4) constitute a relation of the form 𝑷 = ̂ 𝑷 ( 𝑭 , 𝚵 ) (5)which is evaluated in three steps:1. First, the microscale problem with prescribed boundary values given by (2) is solved.2. Next, the solution of the microscale problem is postprocessed to obtain the reaction forces. This can be done either by aboundary integral as shown in (4a), or alternatively by a volume integral over the region of the domain adjacent to theboundary .3. Finally, the reaction forces are combined and scaled according to (4b) to produce the homogenized first Piola-Kirchhoffstress tensor 𝑷 . AVERY
ET AL
FIGURE 2
Representative Surface Element of the membrane within which the entire thickness of the membrane is representedin its entirety.Unfortunately, (5) is not directly compatible with the stated application of interest, namely a plane stress relation of the form(1) expressed in terms of the in-plane components of the Green-Lagrange strain and second Piola-Kirchhoff stress tensors.To formally adapt the homogenization methodology to this setting, we first assume without loss of generality that relation (5)satisfies the principal of material frame invariance which can be stated as follows ̂ 𝑷 ( 𝑸𝑭 ) = 𝑸 ̂ 𝑷 ( 𝑭 ) ∀ 𝑸 ∈ 𝑆𝑂 (3) (6)where 𝑆𝑂 (3) is the group of special orthogonal transformations defined as 𝑆𝑂 (3) = { 𝑸 ∈ ℝ ∶ 𝑸 𝑇 𝑸 = 𝑸𝑸 𝑇 = 𝑰 , det( 𝑸 ) = 1} and 𝐈 is the identity matrix. Regarding the assumption of material frame indifference, under some conditions it can be shownthat (5) is objective and hence (6) holds. For cases in which (6) does not hold, the proposed alternative formulation whichfollows can be interpreted as imposing or restoring material frame invariance which is generally considered to be appropriatefor constitutive relations in solid mechanics.Introducing the polar decomposition of the deformation gradient 𝑭 = 𝑹 𝑼 where 𝑹 ∈ 𝑆𝑂 (3) is the rotation tensor and 𝑼 ∈ ℝ is the right stretch tensor which is positive definite and symmetric, andthe unsymmetric Biot stress tensor defined as 𝑩 = 𝑹 𝑇 𝑷 , it follows from (6) taking 𝑸 = 𝑹 𝑇 that the homogenized constitutive law (5) can equivalently be stated as a relation betweenthe right stretch tensor and the unsymmetric Biot stress using the same functional form, i.e. 𝑩 = ̂ 𝑷 ( 𝑼 , 𝚵 ) . (7)This can be interpreted simply as a variant of the standard transmission conditions (2, 4) in which the right stretch tensor is usedinstead of the deformation gradient to compute the microscale prescribed boundary deformations, and the homogenized stresstensor obtained by evaluating the constitutive function ̂ 𝑷 is identified as the Biot measure rather than the first Piola-Kirchhoff.Specifically, 𝒙 || 𝜕 ℬ 𝒖 = 𝑿 || 𝜕 ℬ 𝒖 𝑼 (8a) 𝑩 = 1 | ℬ | 𝑿 𝑇 ||| 𝜕 ℬ 𝒖 𝒇 || 𝜕 ℬ 𝒖 . (8b)A more convenient relation between the Green-Lagrange strain and the first Piola-Kirchhoff stress tensor can be obtainedfrom (7) by utilizing well-known transformations : VERY
ET AL • First, the right stretch tensor can be obtained from the Green-Lagrange strain using 𝑼 = ( 𝐂 ) = 𝑖 =3 ∑ 𝑖 =1 𝜆 𝑖 𝑁 𝑖 ⊗ 𝑁 𝑖 (9)where 𝐂 = 2 𝑬 + 𝐈 is the right Cauchy-Green deformation tensor, and 𝜆 𝑖 and 𝑁 𝑖 are the eigenvalues and eigenvectors,respectively, of 𝐂 . • Second, the second Piola-Kirchhoff stress can be obtained from the Biot stress using . ( 𝑺 𝑼 + 𝑼 𝑺 ) = 𝑻 (10)where 𝑻 is the symmetric part of the Biot stress tensor 𝑻 = 0 . ( 𝑩 + 𝑩 𝑇 ) . Note that (10) has the form of the Lyapunovequation whose solution is given by a linear system of equations, namely vec ( 𝑺 ) = [ 𝐈 ⊗ 𝑼 + 𝑼 𝑇 ⊗ 𝑰 ] −1 vec ( 𝑻 ) (11)where ⊗ denotes the Kronecker product, and vec ( ⋅ ) denotes vectorization. For example, vectorization of 𝑺 is given by vec ( 𝑺 ) = [ 𝑆 (11)0 𝑆 (21)0 𝑆 (31)0 𝑆 (12)0 𝑆 (22)0 𝑆 (32)0 𝑆 (13)0 𝑆 (23)0 𝑆 (33)0 ] 𝑇 . The dimension of Eq. (11) can be reduced further to six due to symmetry.Substituting (9) and (11) into (7) produces a constitutive function relating the macroscale Green-Lagrange strain and secondPiola-Kirchhoff stress tensor of the form 𝑺 = ̂ 𝑺 ( 𝑬 , 𝚵 ) (12)which is evaluated in five steps:1. First, the macroscale right stretch tensor 𝑼 is computed from the Green-Lagrange strain 𝑬 .2. Second, the microscale problem with prescribed boundary values given by (8a) is solved.3. Third, the solution of the microscale problem is postprocessed to obtain the reaction forces.4. Next, the reaction forces are combined and scaled according to (8b) to produce the homogenized unsymmetric Biot stresstensor 𝑩 .5. Finally, we take the symmetric part of the Biot stress tensor and solve the Lyapunov equation (11) get the homogenizedsecond Piola-Kirchhoff stress tensor 𝑺 . The three-dimensional constitutive law (12) can be adapted to plane stress (and hence membrane elements) using a variant ofthe method proposed by Klinkel and Govindjee for using finite strain three-dimensional material models in beam and shellelements. This method involves solving a local nonlinear equation using Newton’s method to enforce the plane stress condition.Specifically, the requirement that the out-of-plane components of the second Piola-Kirchhoff stress tensor are zero, i.e. 𝑺 𝑧 = [ 𝑆 (33)0 𝑆 (13)0 𝑆 (31)0 ] 𝑇 = 0 , (13)is enforced by iteratively solving for the corresponding out-of-plane components 𝑬 𝑧 of the Green-Lagrange strain tensor whichare treated as unknowns. Each Newton iteration incurs a single evaluation of the three-dimensional constitutive function (12)and its constitutive tangent. Solving the plane stress equation (13) for 𝑬 𝑧 given 𝑬 𝑚 , and then evaluating the in-plane componentsof the second Piola-Kirchhoff stress tensor 𝑷 𝑚 at the resulting configuration corresponds to the evaluation of a plane stressconstitutive relation of the form (1) which can be used as a drop-in replacement for a conventional finite strain plane stressconstitutive equation. This will be demonstrated in what follows using the general purpose finite element analyzer AERO-S.To complete the description of this multiscale material model we note that for a static analysis or a dynamic analysis usingan implicit time-stepping scheme, the consistent constitutive tangent of the plane stress constitutive law is typically required.This quantity, namely 𝜕 ̂ 𝑺 𝑚 ∕ 𝜕 𝑬 𝑚 is readily obtained using the constitutive tangent of the three-dimensional constitutive law; theprecise definition is given in . AVERY
ET AL
In this section, the discretized form of the equations governing the multiscale problem of interest are presented, notably includingcontact at both scales. Specifically, • At the macroscale we seek the solution of a dynamic contact problem. The deforming bodies are discretized in space usingmembrane finite elements, and in time using the explicit central difference time-integration scheme. The contact part ofthe problem is solved using an implicit approach . • At the microscale we seek the solution of static contact problems with nonhomogeneous prescribed boundary displacementand discretized in space using solid finite elements.With regards to notation, a distinction is made in this work between unconstrained degrees of freedom (dofs), i.e., dofs that arenot constrained by any essential boundary condition, and constrained dofs, i.e., dofs that are constrained by essential boundaryconditions. A matrix or vector defined over the set of unconstrained dofs is not designated by any specific symbol. However, avector of constrained dofs is designated by the ring symbol as in ̊ 𝐯 , and a vector defined over the entire set of constrained andunconstrained dofs is designated by the overline symbol as in ̄ 𝐯 . In other words ̄ 𝐯 = [ 𝐯 ̊ 𝐯 ] (14)We assume, without loss of generality, that the discrete form of the macroscale governing equations can be written as adifferential-algebraic inequality (DAI): 𝑴 ̈ 𝒖 ( 𝑛 +1)0 + 𝒇 𝑖𝑛𝑡 ( ̄ 𝒖 ( 𝑛 +1)0 ) + 𝑮 ( ̄ 𝒖 ( 𝑛 +2)0 ) 𝝀 ( 𝑛 +1)0 = 𝒇 𝑒𝑥𝑡 ( 𝑡 ( 𝑛 +1) ) (15a) 𝒈 ( ̄ 𝒖 ( 𝑛 +2)0 ) ≥ (15b) 𝝀 ( 𝑛 +1)0 ≤ (15c) 𝝀 ( 𝑛 +1)0 𝑇 𝒈 ( ̄ 𝒖 ( 𝑛 +2)0 ) = 0 (15d)where 𝑴 is the (diagonal) mass matrix, 𝒇 𝑖𝑛𝑡 and 𝒇 𝑒𝑥𝑡 are the internal and external force vectors, 𝒖 ( 𝑛 +1)0 and ̈ 𝒖 ( 𝑛 +1)0 are thedisplacements and accelerations at time 𝑡 ( 𝑛 +1) , 𝒈 is the gap, a vector-valued constraint function representing the discretizednon-penetration condition, 𝑮 is the transpose of the constraint Jacobian matrix 𝑮 = [ 𝜕 𝒈 𝜕 𝒖 ] 𝑇 and 𝝀 ( 𝑛 +1)0 is a vector of Lagrange multipliers at time 𝑡 ( 𝑛 +1) . The dependence of the internal force vector on the history variablesis acknowledged but not explicitly stated.Given some initial values 𝒖 ( 𝑛 )0 , ̇ 𝒖 ( 𝑛 )0 and ̈ 𝒖 ( 𝑛 )0 at time 𝑡 ( 𝑛 ) , the solution at time 𝑡 ( 𝑛 +1) is obtained using the following updatingprocedure:1. Update displacement 𝒖 ( 𝑛 +1)0 = 𝒖 ( 𝑛 )0 + Δ 𝑡 𝑛 ̇ 𝒖 ( 𝑛 )0 + 0 . 𝑡 𝑛 ̈ 𝒖 ( 𝑛 )0
2. Update acceleration and velocity using the predictor-corrector iterative method(a) predictor: 𝑘 = 0 ̈ 𝒖 ( 𝑛 +1) , = 𝑴 −10 [ 𝒇 𝑒𝑥𝑡 ( 𝑡 ( 𝑛 +1) ) − 𝒇 𝑖𝑛𝑡 ( ̄ 𝒖 ( 𝑛 +1)0 )] ̇ 𝒖 ( 𝑛 +1) , = ̇ 𝒖 ( 𝑛 )0 + 0 . 𝑡 𝑛 [ ̈ 𝒖 ( 𝑛 )0 + ̈ 𝒖 ( 𝑛 +1) , ] 𝒖 ( 𝑛 +2) , = 𝒖 ( 𝑛 +1)0 + Δ 𝑡 𝑛 +1 ̇ 𝒖 ( 𝑛 +1) , + 0 . 𝑡 𝑛 +1 ̈ 𝒖 ( 𝑛 +1) , (b) corrector iterations: 𝑘 = 1 , … ̈ 𝒖 ( 𝑛 +1) ,𝑘 = ̈ 𝒖 ( 𝑛 +1) ,𝑘 −10 + Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 ̇ 𝒖 ( 𝑛 +1) ,𝑘 = ̇ 𝒖 ( 𝑛 +1) ,𝑘 −10 + 0 . 𝑡 𝑛 Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 𝒖 ( 𝑛 +2) ,𝑘 = 𝒖 ( 𝑛 +2) ,𝑘 −10 + 0 . [ Δ 𝑡 𝑛 Δ 𝑡 𝑛 +1 + Δ 𝑡 𝑛 +1 ] Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 VERY
ET AL At each corrector iteration, the acceleration increment Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 is obtained by linearizing the gap function 𝒈 and solving thelinearized sub-problem 𝑴 Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 + 𝑮 ( ̄ 𝒖 ( 𝑛 +2) ,𝑘 −10 ) 𝝀 ( 𝑛 +1) ,𝑘 = − ̃ 𝒇 𝑘 −10 (16a) 𝑮 ( ̄ 𝒖 ( 𝑛 +2) ,𝑘 −10 ) 𝑇 Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 ≥ − ̃ 𝒈 𝑘 −10 (16b) 𝝀 ( 𝑛 +1) ,𝑘 ≤ (16c) 𝝀 ( 𝑛 +1) ,𝑘 𝑇 [ 𝑮 ( ̄ 𝒖 ( 𝑛 +2) ,𝑘 −10 ) 𝑇 Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 + ̃ 𝒈 𝑘 −10 ] = 0 (16d)where: ̃ 𝒇 𝑘 −10 = 𝑴 [ ̈ 𝒖 ( 𝑛 +1) ,𝑘 −10 − ̈ 𝒖 ( 𝑛 +1) , ] ̃ 𝒈 𝑘 −10 = 2Δ 𝑡 𝑛 Δ 𝑡 𝑛 +1 + Δ 𝑡 𝑛 +1 𝒈 ( ̄ 𝒖 ( 𝑛 +2) ,𝑘 −10 ) The corrector sub-problem (16) has the form of a quadratic program and can be solved by the primal-dual active set method
1. Initialize Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 , 𝝀 ( 𝑛 +1) ,𝑘
2. Iterate • Choose active set: = { 𝑖 ∶ [ 𝝀 ( 𝑛 +1) ,𝑘 ] 𝑖 > [ 𝑮 ( ̄ 𝒖 ( 𝑛 +2) ,𝑘 −10 )] 𝑇𝑖 Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 + [ ̃ 𝒈 𝑘 −10 ] 𝑖 < } • Set inactive Lagrange multipliers to zero: [ 𝝀 ( 𝑛 +1) ,𝑘 ] 𝑖 = 0 ∀ 𝑖 ∉ • Solve for Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 and active Lagrange multipliers: 𝑴 Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 + 𝑮 ( 𝒖 ( 𝑛 +2) ,𝑘 −10 ) 𝝀 ( 𝑛 +1) , ,𝑘 = − ̃ 𝒇 𝑘 −10 (17a) 𝑮 ( 𝒖 ( 𝑛 +2) ,𝑘 −10 ) 𝑇 Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 = − ̃ 𝒈 ,𝑘 −10 (17b)where the superscript applied to a vector indicates its restriction to the active set. Similarly, the superscript applied to a matrix indicates its column-wise restriction to the active set.The active set method sub-problem (17) is a linear saddle-point system. To solve for the active Lagrange multipliers, we firsteliminate Δ ̈ 𝒖 ( 𝑛 +1) ,𝑘 and then solve the remaining Schur complement system [ 𝑮 𝑇 𝑴 −10 𝑮 ] 𝝀 = ̃ 𝒈 − 𝑮 𝑇 𝑴 −10 ̃ 𝒇 . (18)To simplify notation, the superscripts denoting time-step index and predictor-corrector iteration have been omitted here but canbe inferred from (17). After solving (18) for the Lagrange multipliers, the acceleration increment can be obtained from (17a).If 𝑮 is rank-deficient then the active set iterations may not converge. In this case, a penalty parameter ( 𝜇 ) can be used toregularize the system, leading to the perturbed systems of the form [ 𝑮 𝑇 𝑴 −10 𝑮 + 1 𝜇 𝑰 ] 𝝀 = ̃ 𝒈 − 𝑮 𝑇 𝑴 −10 ̃ 𝒇 or equivalently, [ 𝑴 + 𝜇 𝑮 𝑮 𝑇 ] Δ ̈ 𝒖 = − ̃ 𝒇 − 𝜇 𝑮 ̃ 𝒈 . This completes the description of the macroscale discrete governing equation and its solution algorithm. Significantly, eachtime-step incurs only one evaluation of 𝒇 𝑖𝑛𝑡 which in the context of a multiscale simulation invariably dominates the compu-tational cost of the entire time-step. In order to evaluate this discrete vector of internal forces, the homogenized stress tensormust be computed at each Gauss point of the macroscale finite element model which in turn involves the iterative solution of theKlinkel-Govindjee plane stress equation with one solution of the discrete microscale governing equation required per iteration. AVERY
ET AL
In the presence of contact at the microscale – for example, non-penetration and sliding of yarns in a woven fabric – the discreteform of the microscale governing equation has a similar form to that of the macroscale (15) but without the time-dependenceand associated temporal discretization. The external force term is also identically zero and can be omitted; the problem is insteaddriven by prescribed values of the constrained dofs and can be described as follows 𝒇 𝑖𝑛𝑡 ( ̄ 𝒖 ) + 𝑮 ( ̄ 𝒖 ) 𝝀 = 0 (19a) 𝒈 ( ̄ 𝒖 ) ≥ (19b) 𝝀 ≤ (19c) 𝝀 𝑇 𝒈 ( ̄ 𝒖 ) = 0 (19d)The quantities 𝑭 , ̄ 𝒖 , 𝒈 , 𝑮 , and 𝝀 are all microscale counterparts of the corresponding macroscale quantities defined previ-ously. This problem can be tackled in similar fashion to that of the macroscale by solving a series of linearized sub-problems ofthe form 𝑲 𝑡𝑔𝑡 Δ 𝒖 𝑘 + 𝑮 ( ̄ 𝒖 𝑘 −11 ) 𝝀 𝑘 = − ̃ 𝒇 𝑘 −11 (20a) 𝑮 ( ̄ 𝒖 𝑘 −11 ) 𝑇 Δ 𝒖 𝑘 ≥ − ̃ 𝒈 𝑘 −11 (20b) 𝝀 𝑘 ≤ (20c) 𝝀 𝑘 𝑇 [ 𝑮 ( ̄ 𝒖 𝑘 −11 ) 𝑇 Δ 𝒖 𝑘 + ̃ 𝒈 𝑘 −11 ] = 0 (20d)where 𝑲 𝑡𝑔𝑡 is the microscale tangent stiffness matrix 𝑲 𝑡𝑔𝑡 = 𝜕 𝒇 𝑖𝑛𝑡 𝜕 𝒖 . Problem (20) can again be solved by the dual-primal active set method proposed for the corresponding macroscale problem(16), although numerous alternatives exist.The computational homogenization method described in Section 2 provides a very general framework for solving the problemof interest without resorting to any ad-hoc approximation. However, without introducing any further approximation, the frame-work – although amenable to parallel implementation – is impractical for all but the most modest of applications due to itscomputational complexity. For example, we estimate that to simulate the inflation of a parachute using a macroscale modelcomprising 182,554 nodes and 279,025 triangular membrane elements would require 49,604,444,444 constitutive function eval-uations and a total run time of approximately 479.95 years for a parallel execution using 96 processing units. Hence, to achievecomputational tractability we propose an alternative regression-based surrogate modeling methodology which is described andevaluated in what follows. We emphasize that this methodology relies exclusively on the general framework that has beenpresented in this section to obtain “training data” that can be used to construct a low-dimensional surrogate model.
Here, a methodology featuring a regression-based surrogate model is presented for dramatically accelerating the solution ofnonlinear dynamic multiscale problems modeled using the multiscale formulation based on the concept of a locally attachedmicrostructure overviewed above. The methodology features a novel training strategy based on the concept of a coupon testanalogy. Regression-based surrogate models can be loosely classified as follows:1. Models whose forms are determined à-priori and whose parameters are fitted to available data. Examples of such modelsare: • The St. Venant-Kirchhoff hyperelastic model, a 2-parameter model characterized by a linear relationship betweenthe second Piola-Kirchhoff stress and the Green-Lagrange strain. • Hyper-viscoelastic models incorporating a hyperelastic model such as that of St. Venant-Kirchhoff, combined witha viscoelastic component based on a Prony series.2. Models whose forms are not entirely predetermined but rather discovered at least in part by the regression/fitting process.An example of such a model is an artificial neural network (NN). In this case, certain characteristics of the model maystill be specified à-priori , such as the number of hidden layers and the functional form of the activation function.
VERY
ET AL A training strategy, i.e., a procedure for sampling the parameter space and collecting stress and strain data is proposed here forconstructing the regression-based surrogate models described in this section. For this purpose a model of a small coupon of themacroscale is employed; in this work we utilize a coupon comprising a single membrane element. We emphasize that due to theoverwhelming cost of an entire multiscale simulation using a high-dimensional macroscale model, it is not practical to collectdata specifically customized to the target application as is sometimes done to train reduced-order models. However, the rangeof strains to which the coupon model is subject to during the training can be customized to a certain extent for example to targetapplications with small, medium or large deformations. Due to the small size of the coupon macroscale model, it is feasible tocollect data that comprehensively samples in a regular grid the parameter space which is dimension three (recall the microscaleprescribed boundary displacements are obtained by a mapping from the in-plane components of the macroscale symmetricGreen-Lagrange strain tensor). The macroscale strain can be indirectly specified for the purposes of training by prescribingdisplacements on the boundary of the macroscale coupon. Figure 3 shows the deformed configurations and corresponding vonMises stress contours for several sampled points in the parameter space obtained during a training performed for the applicationdescribed in the following section.
FIGURE 3
Microscale training solution snapshots showing von Mises stress contours for selected points in the sampledparameter space: uniaxial tension (top left), biaxial tension (top right), uniaxial compression (bottom left), and shear (bottomright).
In this section, the following regression-based surrogate microscale models are considered: • linear [ 𝑆 (11)0 𝑆 (22)0 𝑆 (12)0 ] 𝑇 = 𝑪 [ 𝐸 (11)0 𝐸 (22)0 𝐸 (12)0 ] 𝑇 AVERY
ET AL • quadratic [ 𝑆 (11)0 𝑆 (22)0 𝑆 (12)0 ] 𝑇 = 𝑪 [ 𝐸 (11)0 𝐸 (22)0 𝐸 (12)0 𝐸 (11)0 2 𝐸 (22)0 2 (2 𝐸 (12)0 ) 𝐸 (12)0 𝐸 (22)0 𝐸 (12)0 𝐸 (22)0 𝐸 (11)0 𝐸 (22)0 ] 𝑇 • neural network (linear model with neural network correction) [ 𝑆 (11)0 𝑆 (22)0 𝑆 (12)0 ] 𝑇 = 𝑪 [ 𝐸 (11)0 𝐸 (22)0 𝐸 (12)0 ] 𝑇 + 𝑵 ([ 𝐸 (11)0 𝐸 (22)0 𝐸 (12)0 ] 𝑇 ) where 𝑪 ∈ ℝ and 𝑪 ∈ ℝ are coefficient matrices, and 𝑵 denotes a fully connected neural network which corrects thelinear model by mapping the strain to a stress correction. We will show the proposed neural network model out-performs othersurrogate models in terms of the training/test errors. Finally, the neural network model is applied to describe the behavior of aDGB parachute inflation during a supersonic Mars landing event. Strain-stress data pairs ( 𝑬 ( 𝑖 )0 , 𝑺 ( 𝑖 )0 ) , 𝑖 = 1 , … , 𝑁 are generated by performing a numerical coupon test 𝑁 times, where 𝑁 is thenumber of training data points. Each coupon test is depicted in Fig. 4, the right triangle fabric coupon piece is of length m.The displacements of the right angle node and all out-of-plane displacement are constrained to be zero, and prescribed in-planedisplacements are applied to the other nodes to generate a specified target strain field. The microscale model (see Section 2.2) issolved at the single Gaussian quadrature point located at the center of the right triangle, which delivers the homogenized strainand stress pairs. ̊ 𝒖 ̊ 𝒖 FIGURE 4
Schematic of the coupon test with prescribed displacements (strain fields).The strain field at the macroscale (fabric coupon) level [ 𝐸 (11)0 𝐸 (22)0 𝐸 (12)0 ] 𝑇 is sampled in a cube of extent [−0 . , .
25] ×[−0 . , .
25] × [−0 . , . . Here, the range of strains is customized to match the application of interest, specifically the inflationof a DGB parachute at typical Mars landing conditions. We sample uniformly in the cube with equidistant points in eachcomponent of strain, which accounts for a total of training data points . Each training data point requires solution of thediscrete governing equations (15) governing the multiscale coupon. To facilitate the sampling procedure, each data point istreated as a time-step of a single multiscale simulation in which the prescribed boundary conditions are varied in time along atrajectory (see Figure 5) that traverses all of the data points. Each time-step can be interpreted as an independent static simulation;alternatively each line segment of the trajectory can be interpreted as being associated with the numerical counterpart of a singlephysical coupon test in which two components of strain are held fixed while the third is varied. Crucially, the converged solutionof the microstructure at the previous data point is used as the initial guess for the iterative solution by Newton’s method at thesubsequent point. In total, the data generation procedure takes about 40 CPU hours. To validate the surrogate models, a further test data points were generated by shifting the trajectory. The data is available at https://github.com/Zhengyu-Huang/Fabric-Data.git
VERY
ET AL FIGURE 5
Trajectory used for sampling training data points.
The neural network considered in the present work consists of one hidden layer for efficient purpose. Both tanh and ReLU areused as activation functions. The loss function is defined as 𝑁 ∑ 𝑖 =1 || 𝑺 ( 𝑖 )0 − ( 𝑬 ( 𝑖 )0 ) || + 𝜆 || 𝜽 || , (21)where represents the surrogate model, 𝜽 denotes the hyperparameters in the surrogate model, and 𝐿 regularization is addedwith the regularization parameter 𝜆 = 10 −4 . However, for the fabric material, the shear stress 𝑆 (12)0 is several magnitudes smallerthan the axial stresses 𝑆 (11)0 and 𝑆 (22)0 (see Figures 6). To capture the shear effect better, an alternative weighted loss function isalso considered 𝑁 ∑ 𝑖 =1 ( 𝑆 (11) , ( 𝑖 )0 − (11) ( 𝑬 ( 𝑖 )0 ) ) + ( 𝑆 (22) , ( 𝑖 )0 − (22) ( 𝑬 ( 𝑖 )0 ) ) + 𝑤 ( 𝑆 (12) , ( 𝑖 )0 − (12) ( 𝑬 ( 𝑖 )0 ) ) + 𝜆 || 𝜽 || , (22)where 𝑤 is a weighting constant set to . In total, six regression-based surrogate models are considered: • linear model, • quadratic model, • NN-tanh model: neural network with six tanh neurons and trained with the unweighted loss function (21), • NN-ReLU model: neural network with six ReLU neurons and trained with the unweighted loss function (21), • NN-ReLU model (weighted 6): neural network with six ReLU neurons and trained with the weighted loss function (22),and AVERY
ET AL • NN-ReLU model (weighted 20): neural network with 20 ReLU neurons and trained with the weighted loss function (22).Both the linear and quadratic models are trained by least squares without regularization ( 𝜆 = 0 ). All neural network models aretrained with the limited-memory BFGS (L-BFGS-B) method with regularization ( 𝜆 = 10 −4 ). We use the line search routinein , which attempts to enforce the Wolfe conditions by a sequence of polynomial interpolations. Note the BFGS is applicablein our case since the data sets are typically small ; for large data set the stochastic gradient descent method is suggested.The total relative training/test errors (excluding the regularization term) and relative errors of each component are reported inTable 1. All neural network models lead to errors one magnitude smaller than the linear or quadratic models. The training/testdata and the predictions from all surrogate models are depicted in Figures 6 and 7 for each of the component-wise relations 𝐸 (11)0 − 𝑆 (11)0 , 𝐸 (22)0 − 𝑆 (22)0 , and 𝐸 (12)0 − 𝑆 (12)0 . The training data shows that the fabric material is flexible with respect to shearingand compression. In particular, the shear stresses are two orders of magnitude smaller than the axial stresses under similarstrains. Furthermore, the 𝐸 (11)0 − 𝑆 (11)0 and 𝐸 (22)0 − 𝑆 (22)0 curves are “flat” when the fabric is compressed; their slopes changesuddenly at zero and remain constant in the stretching regime. Due to these features (especially the discontinuity in the slope) inthe strain-stress relations, neural network models deliver better approximations and out-perform linear and quadratic regressionmodels. It is worth mentioning the shear stress is small but highly nonlinear. Neural networks trained with the unweighted lossfunction focus mainly on the axial stresses and fail to capture the nonlinearity in the shear stress (see Figure 6 and Figure 7).By introducing the weighted loss function, the accuracy in the shear stress prediction is improved. Moreover, increasing thenumber of neurons improves the prediction accuracy (see Table 1). However, it is worth noting that increasing the complexityof the neural network architecture increases computational cost and might lead to over-fitting, especially when the training datais inadequate. linear quadratic NN-tanh NN-ReLU NN-ReLU(weighted 6) NN-ReLU(weighted 20)training set 19.5% 10.4% 1.53% 1.57% 2.55% 1.47% 𝑆 𝑥𝑥 𝑆 𝑦𝑦 𝑆 𝑥𝑦 𝑆 𝑥𝑥 𝑆 𝑦𝑦 𝑆 𝑥𝑦 TABLE 1
Total relative error and relative errors of each component of regression-based surrogate models on both training/testdata sets.Regarding the computational cost, the number of operations for a single model evaluation are 𝑂 (15) , 𝑂 (51) , and 𝑂 (13 × neurons + 15) for the linear, quadratic, and neural network models, respectively. This indicates that for explicit time-integrationschemes whose cost is predominately accounted for by constitutive function evaluations, the simulation cost when a neuralnetwork model is utilized can be up to 6 (6 neurons) or 19 (20 neurons) times more expensive than when the linear modelis utilized. However, for implicit time-integration schemes, the additional cost incurred by the neural network model wouldtypically be substantially less than for explicit schemes due to the cost of the equation solver which contributes a significantportion of the computational cost regardless of which constitutive model is utilized. Finally, the trained NN-ReLU model is applied for describing the canopy behavior in the Fluid-Structure Interaction (FSI)simulation of the inflation dynamics of a DGB parachute system in the low-density, low-pressure, supersonic Martian atmo-sphere . While such a simulation is crucial to the understanding of the effects of the fabric material on the performance ofthe parachute during the deceleration process, its main purpose here is to validate the proposed multiscale fabric model usingflight data from the landing on Mars of NASA’s rover Curiosity.
VERY
ET AL Component Parameter Description ValueCanopy 𝐷 Diameter 15.447 m 𝑡 Thickness 7.6 × −5 m 𝐸 Microscale yarn Young’s modulus 3497 MPa 𝜈 Microscale yarn Poisson’s ratio 0.2 𝜌 𝐶 Density 1154.25 kg m −3 𝛼 Porosity 0.08Suspension lines 𝐿 Length 36.56 m 𝐷 Diameter 3.175 × −3 m 𝐸 Young’s modulus 29.5 GPa 𝜌 𝑆𝐿 Density 1154.25 kg m −3 TABLE 2
Geometric and material properties of a DGB parachute system .To this end, the DGB parachute system that successfully deployed in 2012 for the Mars landing of Curiosity is considered(see Figure 8-left). This aerodynamic decelerator system consists of three main components : • the canopy, which is made of nylon material (see Figure 1), • the suspension lines, which are made of Technora T221 braided cords, and • the reentry vehicle.Its geometric and material properties are listed in Table 2.The simulation discussed herein starts from the line stretch stage where the suspension line subsystem is deployed but thecanopy is folded (see Figure 8-right), and the entire system is prestressed by the folding pattern . The incoming supersonicflow is at the state defined by 𝑀 ∞ = 1 . , 𝜌 ∞ = 0 . kg m −3 , and 𝑝 ∞ = 260 Pa.Since the Martian atmosphere is mainly composed of carbon dioxide, the viscosity of this gas is modeled using Sutherland’sviscosity law with the constant 𝜇 = 1 .
57 × 10 −6 kg m −1 s −1 and the reference temperature 𝑇 = 240 K. The Reynolds numberbased on the canopy diameter is .
06 × 10 . Hence, the flow is assumed to have transitioned to the turbulent regime, which ismodeled here using Vreman’s eddy viscosity subgrid-scale model for turbulent shear flow with model constant 𝐶 𝑠 = 0 . .The in-house Eulerian computational framework with an immersed (embedded) boundary method — Finite Volume methodwith Exact two-phase or two-material Riemann problems (FIVER) is adopted in the present work, due to the largedeformation of the parachute system. This method has previously been successfully employed for the simulation of the failureanalysis of submerged structures subjected to explosions and implosions . It incorporates in the framework a parallel AdaptiveMesh Refinement (AMR) based on newest vertex bisection , which enables the capturing of various interactions betweenthe fluid subsystem, the nonlinear parachute subsystem, and the forebody.The canopy of the DGB parachute consists of band and disk gores. Here, these are discretized by 279,025 geometrically non-linear membrane elements. The suspension line subsystem contains 80 lines, each of which is discretized by 500 geometricallynonlinear beam elements. The reentry vehicle, it is modeled as a fixed rigid body that is embedded, together with the entireaerodynamic decelerator system, in the embedding computational fluid domain (see Figure 8).The aforementioned computational fluid domain is a box of size 200 m ×
160 m ×
160 m. It is initially discretized by a meshcomprising 2,778,867 nodes and 16,308,672 tetrahedra. During the FSI simulation reported below, AMR is applied to track andresolve the boundary layer and flow features. The specified characteristic mesh sizes near the reentry vehicle and canopy are . cm and cm, respectively. The specified characteristic mesh size in the wake and near the shock is cm.Since the canopy is made of nylon fabric with an void fraction, its permeability is modeled using an homogenized porouswall model . Due to the massive self-contact of the parachute canopy during its dynamic inflation, the explicit central dif-ference time-integration scheme is used to advance in time the semi-discrete state of the structural subsystem. A small amountof Rayleigh damping is applied to stabilize the system. The microscale yarn Young’s modulus is roughly estimated from the Young’s modulus of the macroscale nylon. AVERY
ET AL
Number of cores Wall-clock time (h)Training data generation 1 40Training 1 0.01Fluid solver 480 96.14NN-ReLUstructure solver 96 19.52NN-ReLU (weighted 20)structure solver 96 24.82St. Venant-Kirchhoffstructure solver 96 2.96FE -basedstructure solver 96 4204336.55 TABLE 3
Computational costs of each component in the proposed framework for the parachute inflation dynamics simulation.First, a quasi-steady state of the flow past the folded parachute configuration shown in Figure 8-right is computed assumingthat this configuration is rigid and fixed. Using this CFD solution and the aforementioned prestressed state of the structuralmodel of the parachute system as initial fluid and structural conditions, respectively, the FSI simulation of the inflation dynamicsof the DGB parachute is performed in the time-interval [0 , . s. The length of this time-interval is such that it covers theinflation process as well as a few breathing cycles of the DGB parachute system. As stated above, the explicit central differencetime-integrator is applied to advancing in time the semi-discrete structural subsystem. On the other hand, the implicit, 3-pointBDF scheme is applied to time-integrate the semi-discrete fluid state. The fluid and structural discretizations are coupled forthis simulation using the stability-preserving, second-order, time-accurate, implicit-explicit fluid-structure staggered solutionprocedure presented in . The fluid-structure coupling time-step is set to Δ 𝑡 𝐹 ∕ 𝑆 = 10 −5 s.Figure 9 graphically depicts the time-evolutions of the dynamic inflation of the DGB parachute and the flow Mach numberaround it. The parachute is fully inflated at approximately 𝑡 = 0 . s; after this time, it starts the breathing cycles expected froma violent, high-speed, dynamic, inflation process.Figure 10 reports the time-histories of the total drag force predicted by the FSI simulation described above. For reference, thisfigure also includes the time-history of the total drag generated by the parachute system of NASA’s rover Curiosity as measuredduring Mars landing , and the simulation with the classical St. Venant-Kirchhoff model in . The NN-ReLU models deliver astable result, which is in reasonably good agreement with the experimental data. However, the effects of the constitutive relationsof the nylon fabric is rather weak in terms of the drag performance. Figure 11 reports the time-histories of the maximum vonMises stresses, an indicator of material failure, predicted by the FSI simulation described above. The results delivered by NN-ReLU model and NN-ReLU (weighted 20) model, which is trained with the weighted loss function are similar. And this indicatesthe shear effect of the nylon fabric is not significant during the inflation. However the comparison with the result delivered by St.Venant-Kirchhoff model illustrates that flexibility with respect to shearing and compression in the multiscale fabric model leadsto lower von Mises stresses in the parachute “breathing” cycle after the full inflation. This disparity is also shown in Fig. 12,which depicts the time-evolutions of the von Mises stress fields. Although further (experimental) investigation is required toconclude which model is more reliable, this comparison illustrates the potential of multiscale constitutive model for improvingthe prediction of material failure.The computational costs of each component, including both the training procedure and the FSI simulation, are reported inTable 3. It is worth mentioning the the estimated simulation time for the direct FE simulation is also reported in Table 3, andthe speedup of the NN-based surrogate model based on this estimation is about (with training costs included). Thisdemonstrates the strength of NN-based surrogate in the constitutive modeling. The wall-clock time is estimated as the multiplication of the number of FE model evaluations on each CPU and the cost of a single evaluation. VERY
ET AL A general framework has recently been developed for computationally tractable, nonlinear multiscale modeling of membranefabrics. It is enabled by the coherent utilization of several established methodologies • computational homogenization known as Finite Element squared or FE based on the concept of a locally attachedmicrostructure, • numerical enforcement of the membraneâĂŹs plane stress condition, and • regression; a surrogate constitutive model is “discovered” using data generated by many multiscale numerical simulationsof a small fabric coupon.This framework encompasses a cascade of multiscale models ranging from the highest fidelity (without surrogate) to the lowest(linear regression surrogate) and has been demonstrated on the simulation of inflation of a disk-gap-band parachute in supersonicMartian atmospheric entry conditions. In this demonstration, the utilization of a neural network surrogate was able to achievespeedups of approximately relative to the highest fidelity, at an overall cost within one order of magnitude of a conventional(i.e. linear elastic model) material. The proposed discovery of a surrogate constitutive model by means of numerical coupontesting is analogous to the experimental testing procedure used to identify the parameters (e.g. the YoungâĂŹs modulus) ofconventional material models. A highlight of the proposed framework is that while experimental data is typically limited touniaxial tension (occasionally biaxial and/or shear data may also be available), numerical data suffers from no such limitationand we can readily explore the entire parameter space (i.e. all physically admissible combinations of normal and shear strains)during the discovery in order to characterize complex and unconventional materials. However, the generation of training datacan be a substantial cost. We propose an application of model-order reduction (see Appendix A) using a novel in-situ trainingapproach to accelerate this task; this will be the topic of a companion paper . ACKNOWLEDGMENTS
Philip Avery, Daniel Z. Huang, Johanna Ehlers and Charbel Farhat acknowledge partial support by the Jet Propulsion Laboratory(JPL) under Contract JPL-RSA No. 1590208, and partial support by the National Aeronautics and Space Administration (NASA)under Early Stage Innovations (ESI) Grant NASA-NNX17AD02G. Parts of this work were completed at the JPL, CaliforniaInstitute of Technology, under a contract with NASA. Optical microscope photographs were taken by Cheyenne Hua at theJPL. Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors and do notnecessarily reflect the views of JPL or NASA.
APPENDIXA PROJECTION-BASED REDUCED ORDER SURROGATE MICROSCALE MODEL
Here, a reduction/hyperreduction framework is presented for dramatically accelerating the solution of nonlinear dynamic multi-scale problems modeled using the multiscale formulation based on the concept of a locally attached microstructure overviewedabove. This framework constitutes a generalization of Zahr to include • a treatment of contact based on the method originally proposed in and featuring the application of a non-negative matrixfactorization scheme to the construction of a positive reduced-order basis for the contact forces, and • a novel training strategy based on the concept of a coupon test analogy.Unlike in where reduction was considered at all scales, we will restrict the exposition here to reduction of the microscaleonly. Specifically, the Proper Orthogonal Decomposition (POD) method is used to construct a PROM at the microscale, and a AVERY
ET AL computational approach based on the Energy Conserving Sampling and Weighting (ECSW) method is used to hyperreduceit. Training is performed offline (i.e. à-priori ) using a small multiscale coupon model.
A.1 Reduction of the primal unknowns
At the microscale (scale ), the number of primal dofs 𝑛 of the computational model is reduced by searching for the primalsolution 𝒖 of the problem of interest in a carefully constructed low-dimensional subspace, i.e., 𝒖 ≈ 𝑽 𝒚 (A1)where 𝑽 ∈ ℝ 𝑛 × 𝑟 is a Reduced Order Basis (ROB) representing a low-dimensional subspace, 𝒚 ∈ ℝ 𝑟 is the vector ofgeneralized coordinates of 𝒖 in this basis, and 𝑟 ≪ 𝑛 . The ROB is chosen to be orthonormal with respect to the identitymatrix, i.e., 𝑽 𝑇 𝑽 = 𝑰 . In this work, the ROB is constructed using POD and the method of snapshots . To this effect, let { 𝒖 (1)1 , … , 𝒖 ( 𝑚 )1 } be 𝑚 state snapshots at scale , i.e., solutions of (19) for different prescribed boundary displacements. We define the primal snapshotmatrix 𝒀 𝑢 as the 𝑛 × 𝑚 matrix whose columns are comprised of the snapshots 𝒀 𝑢 = [ 𝒖 (1)1 ⋯ 𝒖 ( 𝑚 )1 ] . The ROB 𝑽 is composed of the first 𝑟 principle components, or left singular vectors of 𝒀 𝑢 .Furthermore, it follows from (2) that the microscale constrained dofs’ displacement ̊ 𝒖 also lies in a low dimensional subspaceassociated with a vector of generalized coordinates identified as the column-wise vectorization of the right stretch strain tensor 𝑼 − 𝑰 , i.e. ̊ 𝒖 = 𝚷 ⎡⎢⎢⎣ ̊ 𝑿 ̊ 𝑿
00 0 ̊ 𝑿 ⎤⎥⎥⎦ vec ( 𝑼 − 𝑰 ) = ̊ 𝑽 ̊ 𝒚 (A2)where 𝚷 is a permutation matrix and ̊ 𝑿 is a matrix whose three columns represent the 𝑥 , 𝑦 , and 𝑧 nodal coordinates, respec-tively, of the constrained nodes on the microscale boundary. The definition of ̄ 𝒚 follows from to the notational convention (14);similarly a basis encompassing both unconstrained and constrained dofs can be represented, up to a permutation, as ̄ 𝑽 = [ 𝑽 ̊ 𝑽 ] . The dimensionality of the discrete governing equations (19) is reduced at scale by performing a Galerkin projection, i.e.,substituting (A1) in these equations and projecting the first of them onto the column space of 𝑽 . This leads to the PROM 𝑽 𝑇 𝒇 𝑖𝑛𝑡 ( ̄ 𝑽 ̄ 𝒚 ) + 𝑽 𝑇 𝑮 ( ̄ 𝑽 ̄ 𝒚 ) 𝝀 = 0 (A3a) 𝒈 ( ̄ 𝑽 ̄ 𝒚 ) ≥ (A3b) 𝝀 ≤ (A3c) 𝝀 𝑇 𝒈 ( ̄ 𝑽 ̄ 𝒚 ) = 0 (A3d)Despite the fact that the equations (A3) are characterized by a reduced dimensionality, their solution remains computationallyintensive due to the presence of the nonlinear term 𝒇 𝑖𝑛𝑡 . Indeed, the projection of this term implies that every evaluation of 𝑽 𝑇 𝒇 𝑖𝑛𝑡 requires the reconstruction of the full state using the approximation ̄ 𝑽 ̄ 𝒚 , the integration and assembly of the internalforce vector over the entire computational mesh, and its projection onto the subspace represented by the ROB 𝑽 . Because suchcomputations scale with the size 𝑛 of the high-dimensional model at level , they cannot be performed using limited resourcesor at low computational cost, and much less in real time. Hence, they constitute a substantial bottleneck in the solution of (A3).For this reason, a number of hyperreduction methods have been proposed to overcome this bottleneck introduced by nonlinearterms. For solid mechanics and structural dynamics problems, the ECSW method is our preferred hyperreduction method due toits desirable structure-preserving and numerical stability properties . However, in principle, any other hyperreduction methodcan be used to overcome the aforementioned computational bottleneck. VERY
ET AL As introduced in , the ECSW method amounts to a “mesh reduction” algorithm which samples a set of elements 𝒱 ′1 ⊂ 𝒱 and attributes to each sampled element 𝑒 a positive weight 𝛼 𝑒 > such that ̄ 𝑽 𝑇 ̄ 𝒇 𝑖𝑛𝑡 ( ̄ 𝑽 ̄ 𝒚 ) = ∑ 𝑒 ∈ 𝒱 ( ̄ 𝑽 𝑒 ) 𝑇 ̄ 𝒇 𝑖𝑛𝑡 𝑒 ( ̄ 𝑽 𝑒 ̄ 𝒚 ) ≈ ∑ 𝑒 ∈ 𝒱 ′1 𝛼 𝑒 ( ̄ 𝑽 𝑒 ) 𝑇 ̄ 𝒇 𝑖𝑛𝑡 𝑒 ( ̄ 𝑽 𝑒 ̄ 𝒚 ) = ̄ 𝒇 𝑖𝑛𝑡 𝑟 ( ̄ 𝒚 ) (A4)In the above expressions, the superscript 𝑒 designates the restriction of a global vector or matrix to element 𝑒 , and the so-calledreduced mesh | 𝒱 ′1 | ≪ | 𝒱 | can be computed using Lawson and Hanson’s Non-Negative Least Squares (NNLS) algorithm oralternative L1 minimization algorithms in a training step that seeks to minimize the size of the reduced mesh while maintainingan acceptable approximation error for the ensemble of the training data.In addition to achieving scalability with respect to the size 𝑛 of the PROM only in the computation of the components ofthe internal force vector corresponding to unconstrained dofs, the reduced mesh ensures that scale transmission is performedefficiently, i.e., without any operation whose computational complexity scales with | 𝒱 | . This is evident in the transmissionto finer scales where ̊ 𝒖 𝑒 is required for each 𝑒 ∈ 𝒮 ′1 while with regards to transmission to coarser scales, the homogenizedmacroscopic unsymmetric Biot stress tensor is approximated as vec( 𝑩 ) ≈ 1 | ℬ | ̊ 𝒇 𝑟 ( ̄ 𝒚 ) where 𝒮 ′1 ⊂ 𝒮 denotes the subset of surface elements contained in the reduced mesh 𝒱 ′1 , and ̊ 𝒇 𝑟 in general is the restriction ofthe total vector of reduced forces – both internal and contact – to the constrained generalized coordinates, i.e. ̊ 𝒇 𝑟 = ̊ 𝒇 𝑖𝑛𝑡 𝑟 ( ̄ 𝒚 ) + ̊ 𝑽 𝑇 ̊ 𝑮 ( ̄ 𝑽 ̄ 𝒚 ) 𝝀 although the microscale mesh can in some cases be constructed in such a way that the contact forces will contribute nothing tothis quantity. This requires a separation of at least one element between the contact surface and the boundary to be maintained. A.2 Reduction of the dual unknowns
At the microscale (scale ), the number of dual dofs 𝑛 𝜆 of the computational model can also be reduced by searching for thedual solution 𝝀 of the problem of interest in another carefully constructed low-dimensional subspace, i.e., 𝝀 ≈ 𝑾 𝒛 (A5)where 𝑾 ∈ ℝ 𝑛 𝜆 × 𝑟 𝜆 is a dual reduced order basis representing a low-dimensional subspace, 𝒛 ∈ ℝ 𝑟 𝜆 is the vector of generalizedcoordinates of 𝝀 in this basis, and 𝑟 𝜆 ≪ 𝑛 𝜆 . The dual ROB is chosen such that it has no negative elements.In this work, the dual ROB is constructed using Non-negative Matrix Factorization (NMF). To this effect, let { 𝝀 (1)1 , … , 𝝀 ( 𝑚 )1 } be 𝑚 state snapshots at scale , i.e., solutions of (19) for different prescribed boundary displacements. We define the dual snapshot matrix 𝒀 𝜆 as the 𝑛 𝜆 × 𝑚 matrix whose columns are comprised of the snapshots 𝒀 𝜆 = [ 𝝀 (1)1 ⋯ 𝝀 ( 𝑚 𝑘 )1 ] The dual ROB 𝑾 is composed of the columns of the left factor of the NMF of 𝒀 𝜆𝑘 .The dimensionality of the reduced governing equations (A3) is further reduced at scale by substituting (A5) and projectingthe gap function 𝒈 onto the column space of 𝑾 . Hyperreduction of the internal force is also accounted for. This leads to thePROM 𝒇 𝑖𝑛𝑡 𝑟 ( ̄ 𝒚 ) + 𝑽 𝑇 𝑮 ( ̄ 𝑽 ̄ 𝒚 ) 𝑾 𝒛 = 0 (A6a) 𝑾 𝑇 𝒈 ( ̄ 𝑽 ̄ 𝒚 ) ≥ (A6b) 𝒛 ≤ (A6c) 𝒛 𝑇 𝑾 𝑇 𝒈 ( ̄ 𝑽 ̄ 𝒚 ) = 0 (A6d)Typically, the evaluation of the gap function and its Jacobian does not require a reconstruction of the full state but only itsrestriction to the contact surface. Furthermore, Galerkin projection of the contact force term can be optimized by accounting forthe sparsity of the Jacobian. Specifically, only the row-wise restriction of 𝑮 to the contact surface is non-zero. Nevertheless, AVERY
ET AL these evaluations may still incur a substantial computational cost. In principal, hyperreduction can be applied to further acceleratethe evaluation of the reduced gap function and its Jacobian. This is an active topic of research but is not employed in thepresent work. However, we note that in the case of linear constraints the proposed reduction of the dual variables leads to termsinvolving reduced-order matrices that are precomputable , and as such does not generate any bottleneck in the online solutionof the reduced-order discrete microscale equations. Consequently, just like in the case of any other linear terms, their efficientprocessing does not require any hyperreduction.
References
1. Klinkel S, Govindjee S. Using finite strain 3D-material models in beam and shell elements.
Engineering Computations
International Journal for Numerical Methods in Engineering
Journal of Computational Physics
International journal for numerical methods in engineering
InternationalJournal for Numerical Methods in Engineering
ComputationalMechanics
International Journalfor Numerical Methods in Engineering
Computational Materials Science
Computer Methods in AppliedMechanics and Engineering
Computer methods in applied mechanics and engineering
Engineering Computations
Computer methods in applied mechanics and engineering
Computer methods in applied mechanics and engineering
Computa-tional mechanics
Computer methods in applied mechanics and engineering
Mathematical Elasticity: Volume I: three-dimensional elasticity . North-Holland . 1988.
VERY
ET AL
17. Yvonnet J, Monteiro E, He QC. Computational homogenization method and reduced database model for hyperelasticheterogeneous structures.
International Journal for Multiscale Computational Engineering
Rept. UCRL-CR-125780,Lawrence Livermore National Laboratory, Livermore, CA
SIAM Journal onOptimization
Computer Methodsin Applied Mechanics and Engineering
SIAM Journal on Optimization
SIAM Journal onscientific computing
ACM Transactions on MathematicalSoftware (TOMS) arXivpreprint arXiv:1905.12530 arXiv preprintarXiv:2004.00265
Journal of Spacecraft and Rockets arXiv preprint arXiv:1908.08382
Mechanical property determination for flexible material systems . PhD thesis. Georgia Institute of Technology, NorthAve NW, Atlanta, GA 30332; 2016.32. Vreman A. An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications.
Physics offluids
International Journal for Numerical Methods in Fluids
Journal of Computational Physics
International Journal for Numerical Methods in Fluids AVERY
ET AL
36. Main A, Zeng X, Avery P, Farhat C. An enhanced FIVER method for multi-material flow problems with second-orderconvergence rate.
Journal of Computational Physics
International Journal for Numerical Methods in Engineering
Unified multilevel adaptive finite element methods for elliptic problems . PhD thesis. University of Illinois atUrbana-Champaign, ; 1988.39. Borker R, Huang D, Grimberg S, Farhat C, Avery P, Rabinovitch J. Mesh adaptation framework for embedded boundarymethods for computational fluid dynamics and fluid-structure interaction.
International Journal for Numerical Methods inFluids arXiv preprint arXiv:1907.09632
International Journal forNumerical Methods in Engineering arXiv preprintarXiv:2004.00153
International Journal forNumerical Methods in Engineering
International Journal for NumericalMethods in Engineering
International Journal forNumerical Methods in Engineering
Quarterly of applied mathematics
Solving least squares problems . 15. Siam . 1995.49. Chapman T, Avery P, Collins P, Farhat C. Accelerated mesh sampling for the hyper reduction of nonlinear computationalmodels.
International Journal for Numerical Methods in Engineering
VERY
ET AL FIGURE 6
Reference and predicted strain-stress pairs for the training data set. AVERY
ET AL
FIGURE 7
Reference and predicted strain-stress pairs for the testing data set.
VERY
ET AL inlet outlet outlet outlet outlet outlet suspensionlinesreentry vehiclecanopy Ω ( FIGURE 8
Dynamic supersonic parachute inflation problem: system configuration (left); and embedding computational fluiddomain as well as embedded initial folded configuration (right). AVERY
ET AL
FIGURE 9
Time-evolutions of the deployment of parachute DGB system and the associated flow Mach number field.
VERY
ET AL FIGURE 10
Time-histories of the total drag generated during the dynamic, supersonic parachute inflation process: NASA’sCuriosity rover data (blue); FSI simulation with the NN-ReLU model-based (orange), NN-ReLU (weighted 20) model-based (green), and St. Venant-Kirchhoff (red) constitutive relations. FIGURE 11
Time-histories of the maximum von Mises stresses generated during the dynamic, supersonic parachute infla-tion process: FSI simulation with the NN-ReLU model-based (orange), NN-ReLU (weighted 20) model-based (green), and St.Venant-Kirchhoff (red) constitutive relations. AVERY
ET AL