Inverse analysis of metamaterials and parameter determination by means of an automatized optimization problem

In this paper, a novel parameter determination technique is developed for material models in continuum mechanics aimed at describing metamaterials. Owing to their peculiar mechanical properties and behaviors, such as extreme elasticity or high strength‐to‐weight ratio, metamaterials are of interest to be simulated by reduced‐order modeling by means of the generalized mechanics. Such models incorporate constitutive parameters to be determined; we develop an automatized optimization process specifically for obtaining metamaterials parameters. The process aims at minimizing a mechanically meaningful error function measuring the deviation of the continuum from a detailed description by using the Trust Region Reflective optimization method. The parameter identification procedure is tested for an exemplary extension experiment of a metamaterial, proving to be robust and reliable.

In this paper, a novel parameter determination technique is developed for material models in continuum mechanics aimed at describing metamaterials. Owing to their peculiar mechanical properties and behaviors, such as extreme elasticity or high strength-to-weight ratio, metamaterials are of interest to be simulated by reduced-order modeling by means of the generalized mechanics. Such models incorporate constitutive parameters to be determined; we develop an automatized optimization process specifically for obtaining metamaterials parameters.
The process aims at minimizing a mechanically meaningful error function measuring the deviation of the continuum from a detailed description by using the Trust Region Reflective optimization method. The parameter identification procedure is tested for an exemplary extension experiment of a metamaterial, proving to be robust and reliable.

Finite Element Method (FEM), inverse analysis, metamaterials, optimization INTRODUCTION
The design and construction of structures with peculiar mechanical properties are of paramount interest. The design and topology optimization of structures possessing high toughness such as composite materials [1][2][3][4] has been widely studied. Metamaterials are the family of materials that are architectured with tailored material properties by using an optimized topology with an inner substructure. Pantographic lattice, as a macro-scale example of metamaterials, is a structure consisting of two groups of parallel equispaced beams (fibers) which are connected to each other by cylindrical pivots at their intersections. The two groups of fibers are orthogonal in the initial frame [5] (see Figure 1). These days, samples of this structure are produced by virtue of 3D printing technology [6]. F I G U R E 1 A CAD model of a pantographic structure created in SALOME For describing the pantographic structures, it is possible to employ the Cauchy continuum theory. For a short history of the mechanics theory, see [7]. In this paper, "micro-scale" does not correspond to any specific length-scale but rather it means that at one (or more) smaller length-scale(s) than "macro-scale" (at which the phenomena are perceived by the naked eye), the material is comprised of complex substructures.
In the Cauchy continuum theory, the immediate neighborhood of a particle is described by the first derivative in space. This neighborhood is an infinitesimally small spherical region. In numerical mechanics, the radius of this region should be set to a value smaller than the length-scale of the structure. This makes the analysis of the pantographic lattice computationally costly since it has a micro-scale underlying structure, i.e., the pantographic structure has a micro-scale lengthscale. An efficient solution method is using generalized continuum mechanics [8], i.e., including second or higher gradients of displacement in the model, and thus, having additional material parameters than the classical theories [9,10]. In this approach, we do not need to consider the detailed geometry of the structure. Instead, a homogenized model is developed [11]. In this way, the length-scale will become relatively large. The effect of the microstructure is captured through the new coefficients added to the model [12,13]. In the end, we have a reduced-order homogenized model that is solved fast, numerically as in [14].
In [72], a predictive macro model is presented for characterizing large displacements and large deformations in planar pantographic structures assuming extensible beams. This continuum model is based on a meso-scale to macro-scale homogenization technique for formulating a fully nonlinear beam model that is developed based on a discrete system of extensional and rotational springs and masses. The considered energy expression depends on second-gradient displacement terms. In an improved model [73], the pantographic structure is described in a two-dimensional formulation embedded in three-dimensional space. The model takes account of large out-of-plane motions. There are second space derivatives of the displacement in the deformation energy formulation for considering the in-plane and out-of-plane bending and twisting of fibers. In [74], a numerical identification is performed to fit the parameters of the macro model of planar pantographic structures, where the total stored strain energy and two angles in the structure are considered. It is F I G U R E 2 Left: micro-scale model with the detailed information of the substructure leading to a high computational cost. Right: macro-scale model, i.e., metamaterial representing the same mechanical response of the micro-scale one by using the generalized mechanics leading to a reduced-order modeling with a low computational cost shown that the reduced-order method is able to model the structure accurately with a significantly lower computational cost. In [75], through an inverse analysis, the material parameters in the elastic energy description of the aforementioned improved model are determined. The in-plane and out-of-plane stiffnesses are determined by simulating bias extension and shear tests and employing the Levenberg-Marquardt optimization algorithm. In [76], the numerical identification of the constitutive parameters of the aforementioned predictive macro model is carried out through genetic algorithm [77] and Nondominated Sorting Genetic Algorithm II (NSGA-II) [78] optimization methods. Considering the above-mentioned literature, developing an automatized optimization process for identifying the constitutive parameters of metamaterials is the main objective of the underlying work.
In this research, we aim at implementing the conventional elasticity theory on a pantographic structure, and solve it numerically by following the computational methods by means of the finite element method (FEM) as in [79]. Then a macro-scale homogenized model is employed for modeling the structure at a reduced computational cost. The parameters of the macro-scale model are identified through an automatized optimization process, and the reliability of the model is assessed.
The paper is organized in the following way: in Section 2, we briefly explain the micro-scale nonlinear elastostatic model in continuum mechanics and a second-gradient macro-scale model. In Section 3, we identify the parameters of the macroscale model through an optimization problem. To this aim, a fully automatized optimization problem is developed and implemented in the Python language. In Section 4, the micro-scale and macro-scale models of a pantographic structure are implemented for obtaining the deformation in the structure. In Section 5, we present the results of the simulation of a tensile test on the pantographic structure obtained from micro-and macro-scale models. Also, the results of the parameter determination process are given. In Section 6, the conclusions and perspectives are discussed.

CONTINUUM MECHANICS MODELS
In this section, the micro-scale nonlinear elastostatic model in continuum mechanics and a second-gradient macro-scale model are briefly explained. Consider the micro-scale model with the substructure and its corresponding homogenized model as demonstrated in Figure 2.

Micro-scale model
In the Lagrangean frame, we assume that a continuum body initially occupying  0 ⊂ 3 deforms to the current frame occupying  ⊂ 3 , expressed in Cartesian coordinates. The position of the particles at the initial time, 0 , is identified by and they move to at the current time, , by undergoing the displacement ( , ). Therefore, we have By defining the deformation gradient as = , the Green-Lagrange strain tensor is defined as where we use the Einstein summation convention over repeated indices, comma denotes a partial derivative with respect to in the Lagrangean frame and is the Kronecker delta, which denotes an identity matrix. As the strain measure is nonlinear, the model is valid for large deformations such that it is capable of capturing geometrical nonlinearities.
The density of stored energy, m , is defined as a function of the Green-Lagrange strain tensor and the Lamé parameters of the material, , , as Here, for obtaining the governing equations of the structure, we exploit the formulation based on action principles [80]. An action like energy is used as where the primitive variable is the displacement . In Equation (4), the first term accounts for the kinetic energy where 0 anḋare the mass density at the initial frame and the time derivative of displacement, respectively. The third term, 0 , is the potential energy density (per volume) and is the specific (per mass) gravitational force. The second integral is for the work done on the surfaces wherêis the traction vector applied on the Neumann boundary  0 .
In quasi-static conditions, the continuum body deforms slowly such that the inertial terms are negligible. If the deformation caused by the mass of the body compared to the deformation due to applied force is negligible, the force term can also be ignored. According to the principle of least action, for arbitrary values of the test function , the variation of the action like the functional in Equation (4) vanishes, as follows: Employing the variational formulation leads to the weak form of the micro-scale model as below: where we apply the chain rule for computation of the first term as m , = m , .
By solving the weak form of Equation (6) numerically, the displacement is obtained for the whole body.

Macro-scale model
For characterizing the behavior of planar pantographic structures, the reduced-order model of [72] is employed here. In this model, a homogenization from meso-scale to macro-scale is carried out. In the formulation of this homogenized model, the stretching and bending stiffnesses of the fibers, as well as the shear stiffness of the cylindrical connecting pivots, are considered. The model is briefly explained in the following. The placement maps a material point from in the initial frame to in the current frame through the displacement . In the map ∶ Ω → 2 , Ω is chosen as a rectangular plane region. A Cartesian orthogonal coordinate system, with the unit base vectors of ( , ), is assumed to be aligned in the direction of the two families of fibers (see Figure 1).

=
( 1 2 ) = ( 1 , 2 ) = ( 1 + 1 ( )) + ( 2 + 2 ( )) The deformation gradient is ∇ and the vectors and show the orientation of the fibers after deformation (in the current frame), where || ⋅ || means the 2 norm. In this model, the "stretch of fibers" and the "fiber geodesic curvature" are defined by and , respectively, where, The vectors and contain second-gradient terms of displacement, which are incorporated in the stored energy definition by means of the generalized continuum. This second-order theory is a subset of the strain gradient theory, since contains the second-gradient of which is the gradient of the strain . The shear distortion, , is the change of the angle between two specific fibers from two families of initially orthogonal fibers, defined as = arcsin( ⋅ ).
The strain energy density of the macro-scale model, M , is defined as The strain energy density consists of the energy terms of elongation and geodesic bending of the fibers and the torsion of the connecting pivots. The e , g , s are the stretching, geodesic bending, and the shear stiffnesses, respectively, which are positive and constant constitutive parameters of the macro-scale model. The geodesic bending is the bending that occurs in a fixed plane and it has strain-gradient effects in the continuum theory [81].
Here, analogous to the previous section, we utilize the action principles by utilizing the energy density of the macroscale model in Eq. (15). Based on the principle of least action in Eq. (5), the weak form for macro-scale model is obtained as

OPTIMIZATION
The goal of the optimization problem is to identify the parameters of the macro-scale model by fitting its results with the micro-scale model. In the following, the details of the numerical identification process and the optimization algorithm are discussed.

Numerical identification
The identification of the parameters of the macro-scale model ( e , g , and s in Equation (15)) is done numerically through an optimization problem. The optimization procedure minimizes an error measure yet to be defined. The least restrictive hypothesis is that the micro-scale and macro-scale models have the same deformation energy. However, we emphasize that the displacements in both scales are different, as well as the local deformation energy density, although the total deformation energy is the same, applied by the same work done on boundaries. In a well-posed problem, the boundaries are given such that we expect that the same energy over the boundaries is transferred to the macro-scale model as well as to the micro-scale model. By the total deformation (or stored) energy, we mean the deformation energy stored in the whole structure. We use the energy difference between the micro-and macro-scale models as the error measure for optimization; any other additional measure like matching boundaries or the same curvature could be introduced by Lagrange multipliers [82].

Optimization algorithm
We automatize an optimization problem for determining metamaerials parameters, and implement it in the Python language. We have made use of the optimization submodule of the Scipy package [83]. A flowchart of the optimization process and the numerical identification is shown in Figure 3. In this process, at first, the micro-scale energy, bounds of the parameters, and the convergence tolerances are set as input parameters. Then, for the first iteration of the optimization, we need an initial guess of the parameters for calculating the macro-scale energy. In every iteration, the parameters are updated, until at least one of the convergence criteria stops the iteration, and the solution is the output.
TA B L E 1 Geometrical parameters used for constructing the geometry as suggested in [46], the formula for the initial guess of the parameters as suggested in [75], and the parameters determined the optimization problem The optimization problem is to minimize the sum of squared value of the error that is the difference between the energies at micro-and macro-scale models at every step of loading. The mechanical problem is as follows: Consider a rectangular plate as the macro-scale model representing the pantographic structure at the micro-scale under uni-axial tensile testing. One side is clamped, and the other end is loaded in equal quasi time steps up to the maximum loading. The stored energy of the structure is calculated after every time step of loading. The objective (error) function is defined as where is the number of steps. Then, the following optimization problem gives the best parameters ( e , g , s ) being the unknowns for the minimization problem, The initial guess for the parameters ( e , g , s ) are calculated as a function of the geometry and material of the pantographic structure, as suggested in [75]. The formulas for the initial guess of the parameters are presented in Table 1, where is the area moment of inertia with regard to -axis, and is the shear modulus of the material at the micro-scale. The unknowns are normalized by the values of the initial guess before being set in the optimization code. Therefore, the solution obtained from the optimization will be a factor of the initial guess values, and the parameters are multiplied by the initial guess values before being used in the macro-scale model. The unknown sensitivity of parameters to the error measure is normalized by dividing each parameter to their corresponding initial guess. Hence, the optimization procedure alters them with the same accuracy set by the convergence tolerances.
The employed optimization algorithm, responsible for updating the parameters and distinguishing the path toward the optimum point, is the Trust Region Reflective (TRR) algorithm. The TRR algorithm is an improved version of typical trust region algorithms. Trust region methods are well-known as powerful approaches for solving unconstrained nonlinear minimization problems with strong convergence properties. The idea of a trust region method for minimizing a function is as follows. It approximates the objective function, ( ), with a quadratic function, ( ), which simplifies the objective function in a neighborhood, , around the current point, . The ( ) is usually formed by the first two terms of the Taylor expansion of ( ) around . Then, a sub-problem is defined as computing a trial step and minimizing ( ) within . After every step, if the cost (objective) function at the point + has decreased, the current point is updated to + , and it is called a successful step. Otherwise, the current point is not changed, but the radius of the region is reduced for the next step [84].
In TRR, the trust region idea is generalized from unconstrained to bound-constrained nonlinear minimization problem. Using a novel reflection technique, in TRR, the convergence is accelerated, and it features high computational efficiency and robustness. The TRR is based on the interior reflective Newton algorithm, as proposed in [85]. The interior reflective method does the iterations within a confined interior of the feasible region, which is defined by the upper and lower bound constraints of the variables [86].

F I G U R E 4
The geometrical dimensions depicted on the pantographic structure After every iteration of the optimization, the numerical solution of macro-scale model is repeated, and the convergence of the optimization problem is assessed through three criteria. For this purpose, I) the reduction (change) of the cost function, II) the norm of the step , and III) the infinity norm of the scaled gradient of the objective function are calculated and checked by the criteria using the convergence tolerances. The scaled gradient is calculated by accounting for the presence of the parameter bounds. As soon as any of the above-mentioned criteria is satisfied, the optimization loop terminates, and the best constitutive parameters for the macro-scale model are obtained as the output of the problem.

GEOMETRY AND MODELING
In this section, a pantographic structure is considered, and the micro-scale and macro-scale models, i.e., the weak forms in Equation (6) and Equation (16), respectively, are implemented for obtaining the deformation in the structure. First, the geometrical and material properties of the structure are presented; then, the details of implementation of the models are discussed.

Geometry
The geometry of the pantographic structure considered here is from [46], where the structure was manufactured and inspected through experimental tests. The material of the structure is Polyamide PA 2200, a polymer modeled as an isotropic, linear elastic material with Young's modulus = 1600 MPa and Poisson's ratio = 0.3. The geometrical dimensions of the structure are presented in Figure 4 and Table 1.

Modeling and boundary conditions
For building the micro-scale CAD model of the structure as depicted in Figure 1, the open-source SALOME platform [87] is used. The model is meshed with tetrahedral continuum elements using NETGEN algorithm [88], as shown in Figure 5. two-dimensional (2D) rectangle is created and meshed with triangular elements with around 1300 degrees of freedom. We remark that the computational cost increases exponentially regarding the degrees of freedom. For a uniaxial tensile test, left and right ends are given by the Dirichlet boundary conditions. For both the 3D and the 2D models, the left side of the structure is clamped, and the right side is moved for 37 mm, applied in eight quasi time steps so that we will obtain the solution at every step of displacement.

Numerical implementation
The micro-and macro-scale continuum models discussed in Section 2 are implemented on the respective CAD models. For solving the weak forms of the two models numerically and applying the boundary conditions, the models are implemented into the finite-element code in the FEniCS platform [89], which is an open-source computing platform for automated solution of partial differential equations. The Newton-Raphson method is utilized for linearizing the nonlinear differential equations. For the convergence study of the macro-scale model, see the appendix.

RESULTS AND DISCUSSION
The results of the simulation of a tensile test on the pantographic structure obtained from micro-and macro-scale modelings are presented. The visualizations are done in ParaView [90]. By solving the weak form of the micro-scale model, the deformation of the structure, as the primitive variable of the problem, is obtained. Figure 6 shows the displacement in the pantographic structure for a displacement of 37 mm (17.6% normal strain) in the tensile test. The same boundary conditions for the tensile test are applied on the homogenized model of the macro-scale model. Through the optimization process, the parameters, e , g , s in Eq. (15) are determined and compiled in Table 1 these parameters in the macro-scale model, the deformation and energy density shown in Figure 7 are obtained. In order to compare the consistency of the models, the deformations calculated by the micro-and macro-scale models are plotted on each other in Figure 8.
As shown in Figure 8, the two approaches show very good consistency, and deformation along the boundary is adequate. The values of total stored energy are plotted and compared in Figure 9. As the loading is applied in eight time steps, the energy is calculated at each step. The graph denotes that the optimization algorithm has minimized the error between the energies of macro-and micro-scale models.

F I G U R E 9
The values of total stored energy from the macro-and micro-scale models In order to check the sensitivity of the model, the constitutive parameters are determined after every loading step, as shown in Table 2. The values in Table 2 are given as a factor of the initial guess values. By setting each set of the parameters of Table 2 in the macro-scale model, the energy is calculated and compared in Figure 10. The consistency of the plots in Figure 10 shows the robustness and reliability of the macro-scale model.
The results demonstrate that a reduced-order homogenized model simulates the behavior of a metamaterial with a complex structure which has a detailed substructure, and produces comparable results. The computational cost of the reduced-order model is significantly lower than the detailed three-dimensional model.

CONCLUSION
In this paper, a 3D micro-scale continuum model is implemented on a pantographic structure, and it is solved numerically in the FEniCS platform. On the other hand, a macro-scale homogenized model is employed for modeling the same structure, and the constitutive parameters of this model are identified by developing an automatized process by means of FEM computations and an optimization problem. In this optimization problem, the total stored energies obtained from the two models are compared, and the error is minimized. The optimization algorithm utilized here is the Trust Region Reflective method. The results show that the proposed procedure for parameter determination is robust and proves to be consistent with the micro-scale model. Being the error function defined only in terms of kinematic quantities, the application of this method to experimental model identifications based on full-field kinematic measurements is foreseen.

A C K N O W L E D G M E N T S
The work illustrated in this paper is part of N. Shekarchizadeh's thesis for the Ph.D. course in "Mathematical Models for Engineering, Electromagnetism and Nanosciences" at Sapienza University of Rome. Open access funding enabled and organized by Projekt DEAL.