We show how a theoretical framework developed for modelling nonuniform growth can model the shot peen forming process. Shot peen forming consists in bombarding a metal panel with multiple millimeter-sized shot, that induce local bending of the panel. When ap- plied to different areas of the panel, peen forming generates compound curvature proﬁles starting from a ﬂat state. We present a theoretical approach and its practical realization for simulating peen forming numerically. To achieve this, we represent the panel undergoing peen forming as a bilayer plate, and we apply a geometry-based theory of non-Euclidean plates to describe its reconﬁguration. Our programming code based on this approach solves two types of problems: it simulates the effect of a predeﬁned treatment (the forward problem) and it ﬁnds the optimal treatment to achieve a predeﬁned target shape (the in- verse problem). Both problems admit using multiple peening regimes simultaneously. The algorithm was tested numerically on 200 randomly generated test cases. NOMENCLATURE Ac Current local area Atar Target local area Atotal Total area of the plate in its initial conﬁguration a 22matrix deﬁning the local ﬁrst fundamental form (ar;t;ar;b) Two ﬁrst fundamental forms deﬁning the bilayer rest conﬁguration (aortho r;t;aortho r;b)Two ﬁrst fundamental forms deﬁning the bilayer rest conﬁguration for the case of local in-plane orthotropic growth (anew r;t;anew r;b) Two ﬁrst fundamental forms deﬁning the adjusted bilayer rest conﬁguration (arf;t;arf;b)Two ﬁrst fundamental forms deﬁning the bilayer rest conﬁguration obtained on the ﬁnal iteration b 22matrix deﬁning the local second fundamental form (ac;bc) The ﬁrst and the second fundamental forms deﬁning the monolayer current con- 2 Journal of Manufacturing Science and Engineering ﬁguration (anew c;bnew c)The ﬁrst and the second fundamental forms deﬁning the monolayer current conﬁg- uration, which was obtained by assignment of the adjusted rest conﬁguration and subsequent deformation (af;bf) The ﬁrst and the second fundamental forms deﬁning the monolayer ﬁnal conﬁgu- ration (ainit;binit) The ﬁrst and the second fundamental forms deﬁning the monolayer initial conﬁgu- ration (ar;br) The ﬁrst and the second fundamental forms deﬁning the monolayer rest conﬁgu- ration (anew r;bnew r)The ﬁrst and the second fundamental forms deﬁning the adjusted monolayer rest conﬁguration (aortho r;bortho r)The ﬁrst and the second fundamental forms deﬁning the monolayer rest conﬁgu- ration for the case of local in-plane orthotropic growth (atar;btar) The ﬁrst and the second fundamental forms deﬁning the monolayer target conﬁg- uration Cc Point cloud describing the current conﬁguration Cnew c Point cloud describing the current conﬁguration obtained by assignment of the adjusted rest conﬁguration and subsequent deformation Cf Point cloud describing the ﬁnal conﬁguration Ctar Point cloud describing the target conﬁguration dH The Hausdorff distance EBL Elastic energy of a bilayer EML Elastic energy of a monolayer ~ e0,~ e1;~ e2 Edge vectors deﬁning three edges of a triangular element Hc Current local mean curvature Htar Target local mean curvature h The plate thickness 3 Journal of Manufacturing Science and Engineering I The identity matrix K Number of triangular elements in the model kA Ratio between the target local area and the current local area kH Ratio between the target local mean curvature and the current local mean curva- ture M Maximal number of iterations N Number of predeﬁned peening regimes ~ m Mapping deﬁning position of each point of the mid-surface in R3based on its co- ordinates in U ~ n Unit normal vector ~ n0,~ n1;~ n2 Edge-directors deﬁning normals to the three edges of a triangular element ~ r Radius-vector deﬁning position of each point of the plate in space S The shape operator Sc Current shape operator Sr Rest shape operator Snew r Adjusted rest shape operator Star Target shape operator U Domain of R2englobing the plate mid-surface in coordinates (x;y) U1 Sub-domain of U ~ v0,~ v1;~ v2 Position vectors deﬁning three vertices of a triangular element (x;y) Lagrangian curvilinear coordinates parametrizing the plate mid-surface Y Y oung modulus z The through-thickness coordinate  Total eigenstrain 1 The ﬁrst eigenstrain moment  The upper threshold for jkHj " Eigenstrain tensor induced by the peening treatment "elStrain tensor reﬂecting the elastic material response 4 Journal of Manufacturing Science and Engineering "resResidual strain tensor ("t;"b) Eigenstrain induced in the top and bottom layers of a bilayer respectively ("t avg;"b avg) The average of the orthotropic eigenstrain in two principal directions on the top and bottom layers respectively ("cen;t;"cen;b)A group centroid on the plane ("t;"b) ("rf;t;"rf;b) Eigenstrain prescribed to the top and bottom layers by (arf;t;arf;b) (t;b) Angles between the local principal strain direction and the xaxis on the top and bottom layers respectively  Poisson’s ratio Threshold for the dimensionless error between the current and the target conﬁgu- rations used as a stop condition '0,'1;'2 Angles of inclination of the edge-director vectors ~ n0,~ n1; ~ n2with respect to the average of the adjacent face normals Dimensionless error between the ﬁnal and the target conﬁgurations INTRODUCTION Shot peen forming is a cost-effective technology for shaping large metal plates, such as air- plane wing skins, without dies. It consists in bombarding the surface of a component with a large number of millimeter sized shot made of steel, glass or ceramic. The velocity of a shot is sufﬁ- ciently high to plastically deform the upper layer of the plate upon impact and to stretch the plate locally. This effect causes local bending of thin components and leads to a convex curvature on the peened side [1]. Repeated impacts also induce a ﬁeld of compressive residual stress, that can improve fatigue life [2]. When developping a shot peen forming process, one is faced with two types of problems as schematized in Fig. 1: the forward problem andthe inverse problem [3]. The forward problem is formulated with the following question: which shape will the component adopt if it is peened ac- cording to a given pattern? The inverse problem denotes the following: given an initial shape of the component and the target shape, how should one peen the component to make it deform into the 5 Journal of Manufacturing Science and Engineering Fig. 1. The two problems of the shot peen forming simulation. a) The forward problem consists in the determination of the ﬁnal shape of the component given its initial shape and the peening pattern. b) The inverse problem consists in determination of the optimal peening pattern given the initial shape of the component and its target shape. target shape? A numerical solver for both problems is necessary to optimize the forming process. Thus, without numerical resolution of the inverse problem, the design of a peen forming procedure for each new component is a craft trial-and-error process plagued with risk and uncertainty. It lasts up to several months and implies many scrapped parts. On the other hand, numerical resolution of the forward problem is necessary to check the quality of the inverse problem resolution and to simulate the effect of additional treatments. A straightforward simulation of individual peening impact, such as the one conducted in [4], is precise but computationally expensive. For this reason, simpliﬁed multiscale simulation ap- proaches, such as the eigenstrain approach, were developed. It implies formulating the applied peening loads in terms non-elastic strains imposed on the component [5] [6]. The eigenstrain 6 Journal of Manufacturing Science and Engineering approach represents the treated plate as a thin bilayer where each layer undergoes a nonuniform plastic in-plane swelling or shrinking. The forward problem in this formulation can be numerically solved using shell ﬁnite element models [7] [8]. In addition, such models allow to implement an inverse problem resolution algorithm based on the topology optimization methods [9] [10]. A sim- ilar inverse problem resolution approach was also applied in [11] in the context of the laser peen forming process. Given that this method is based on numerical optimization, its speed decreases with the number of elements in the model. On the other hand, the inverse problem can be solved using an artiﬁcial neural network [8]. The neural network provides near-perfect accuracy and fast calculation on-line. However, it requires the generation of a large ﬁnite-element forward problem solution database and a long training phase for each new plate geometry. To overcome these issues during the inverse problem resolution, we turned to the theory of non-Euclidean plates [12] [13]. This theoretical framework lies at the intersection of mechanics and differential geometry. It precisely describes distortion of multilayer plates induced by the pre- scription of a nonuniform non-elastic strain. Prescription of such strain makes the surface metric non-Euclidean so that it does not satisfy the compatibility conditions of the Euclidean space, which gives name to the theory [13]. Experiments conducted in the ﬁeld of 4D printing with elastic poly- mer sheets have proven the accuracy of the theory of non-Euclidean plates in numerical forward problem resolution. Thus, this theory precisely predicted curvature of spherical, cylindrical and saddle shapes grown out of a ﬂat state by induction of a nonuniform plastic strain [14] [15]. Nu- merical simulation of growth for these three cases also showed good accordance with analytical shapes [16]. Moreover, this theory predicted the form of more complex shapes, such as heli- coid, catenoid or an orchid ﬂower grown out of a ﬂat state using polymer bilayers with oriented ﬁlaments [17]. The theory of non-Euclidean plates also provides instruments for an efﬁcient in- verse problem resolution. Such algorithms for the case of polymers with oriented local growth are reported in [17] and [18]. In this paper we propose to use the eigenstrain approach to represent the treated plate as a bilayer, and we resort to the theory of non-Euclidean plates to calculate the resulting distortion of the bilayer. To numerically solve the inverse problem, we created an iterative algorithm that implies 7 Journal of Manufacturing Science and Engineering resolving the forward problem on each iteration. The algorithm adjusts the peening pattern based on the discrepancy between the current shape, i.e., the shape obtained with the current pattern, and the target shape. The adjustment is done on a local scale using simple arithmetic operations, so the computation speed on this stage stays constant with an increasing number of variables. This stage does not require any preliminary training phase either. The algorithm constitutes a general approach for the inverse problem resolution in case of bilayers subjected to isotropic, i.e. non-oriented, local growth, such as those examined in [15]. When the iterative adjustment is ﬁnished, we group the peening pattern to make it practically applicable. In other words, we divide the pattern into zones treated with constant peening regimes. The number of available peening regimes and their intensities are pre-determined based on the practical constraints. We start this paper with the theoretical background section. First, we examine the eigenstrain approach that relates peen formed plates and swelling non-Euclidean bilayers. Next, we move to the theory of non-Euclidean plates, namely to the geometrical shape description and the forward problem resolution method that it implies. In this section we also formulate the inverse problem in terms of the theory of non-Euclidean plates. We then pass to the methodology section by present- ing our inverse problem resolution algorithm and an approach for its numerical implementation. The grouping and validation strategies are presented in the same section. The results of the val- idation campaign are presented subsequently, and ﬁnally the advantages and limitations of our approach are discussed. THEORETICAL BACKGROUND The eigenstrain approach and strain decomposition The term eigenstrains denotes all non-elastic strains arising in the material, such as plastic, thermal or piezoelectric strains [5]. The only type of eigenstrain generated by peen forming is the plastic strain. Indeed, numerous overlapping impact indentations plastically stretch the outer layer of material, and the rest of the material responds to this newly-introduced eigenstrain with the emergence of stress. In order to conserve its integrity and to balance the stress, the plate deforms elastically. In case of small strains the residual strain tensor "resis additively decomposed into the 8 Journal of Manufacturing Science and Engineering eigenstrain part "and the elastic part "el( [19], Appendix): "res="+"el: (1) This relation holds for peen forming because the process deals with thin plates and the peening- induced strains are small [4]. The elastic strain may affect the whole plate, while the eigenstrain is present only in the stretched outer layer. The thickness of the plastically deformed outer layer varies depending on the peening parameters and the treated material. With the eigenstrain ap- proach, the resolution of the forward problem for shot peen forming consists in introducing the eigenstrain over the whole shot peened area and determining the elastic springback. The through-thickness eigenstrain proﬁle We endow the mid-surface of the plate with two Lagrangian coordinates (x;y), and we assign a Lagrangian coordinate zin the through-thickness direction. The Lagrangian coordinates follow the plate as it deforms. We assume that the material is plastically isotropic, so the eigenstrains are the same in all in-plane directions: "xx(x;y;z ) ="yy(x;y;z ). Also,"zz=2"xxdue to plastic in- compressibility. For a small area around a point (x0;y0)on the mid-surface, the through-thickness eigenstrain proﬁle "xx(x0;y0;z) ="yy(x0;y0;z)can be measured directly using the X-ray diffrac- tion method [20]. Otherwise, it can be deduced from the residual stress proﬁle, which is deter- mined with such methods as hole drilling [21], layer-removal [21] [22] or the two cut compliance method [23]. In this case, the through-thickness eigenstrain proﬁle is reconstructed numerically from the measured residual stress proﬁle [7] [24]. Mechanically, the introduction of the eigenstrain can be modeled as slicing the plate into thin layers, stretching the outer layers separately following the eigenstrain proﬁle, and then gluing everything back. To numerically simplify the problem, we virtually idealize the eigenstrain proﬁle by assuming that the plate consists of two layers of equal thickness that can separately undergo nonuniform in-plane swelling or shrinking. The idealized eigenstrain proﬁle leads to the same 9 Journal of Manufacturing Science and Engineering deformed shape as the real one. Figure 2 illustrates both proﬁles induced on a uniformly treated plate. Mathematically, the idealization procedure consists in ﬁnding the local eigenstrain "t(x;y) to be introduced in the top layer and the local eigenstrain "b(x;y)to be introduced in the bottom layer. It is done by equating the total eigenstrain (x;y)and the ﬁrst eigenstrain moment 1(x;y) induced locally by the real eigenstrain proﬁle to those induced locally by "t(x;y)and"b(x;y)[7]. In the general case, and1are deﬁned as:  (x;y) =Zh=2 h=2"xx(x;y;z )dz; (2) 1(x;y) =Zh=2 h=2"xx(x;y;z )zdz; (3) wherehstands for the plate thickness. For the idealized bilayer proﬁle and1are expressed as:  (x;y) =h 2 "t(x;y) +"b(x;y) ; (4) 1(x;y) =h2 8 "t(x;y)"b(x;y) : (5) The idealized eigenstrain is positive on the side that undergoes the peening treatment and is negative on the other side. It should be noted that although "tand"bgive rise to the same 10 Journal of Manufacturing Science and Engineering in-plane extension and curvature as the real eigenstrain proﬁle, the idealized one generates a different residual stress proﬁle. The theory of non-Euclidean plates applied to the modeling of shot peen forming The theory of non-Euclidean plates allows to numerically determine the elastic response of thin bodies to an applied nonuniform non-elastic strain, e.g., biological growth or eigenstrain. The introduction of such strain leads to the emergence of stresses, which are entirely eliminated only if the plate adopts a so-called rest conﬁguration. The rest conﬁguration implies "res=", so it is generally not realizable without the loss of integrity. Instead of the rest conﬁguration the plate adopts an integral ﬁnal conﬁguration (ﬁnal shape) , which is still residually stressed. The theory of non-Euclidean plates relates the rest and the ﬁnal conﬁgurations through elastic energy and uses tools from differential geometry to describe the shape of plates. We adopted the approach described in [19] for its numerical implementation. Geometrical shape description In the framework of the theory of non-Euclidean plates, the plate shape is associated with the shape of its mid-surface [19]. We denote by Uthe domain of the plane containing the coordinates (x;y)that parametrize the mid-surface: (x;y)2UR2. The position of each point of the mid- surface in a 3D space is deﬁned by the mapping ~ m:U!R3. We adopt the Kirchhoff-Love assumptions, so the position ~ rof a point belonging to the plate is expressed as: ~ r(x;y;z ) =~ m(x;y) +z~ n(x;y); (6) where~ nis the unit normal vector. The mid-surface shape is described by the ﬁrst and the second fundamental forms , that are binary quadratic forms associated with a symmetric 22matrix. Both fundamental forms are local quantities varying smoothly along the surface. The ﬁrst fundamental form describes changes in the length of curves and areas of regions on the surface. In other 11 Journal of Manufacturing Science and Engineering      Fig. 2. The real and idealized eigenstrain proﬁles induced by uniform shot peening of a plate. We denote the plate thickness as h. The coordinate zgoes along the thickness and measures form the midplane. a) The plate undergoes uniform shot peening (side view). b) The through-thickness eigenstrain proﬁle is nonuniform along z, and its peak value "maxis close to the surface [7]. Its effect is accurately simulated by virtually dividing the plate into thin layers and imposing different eigenstrain to each of the layers thus reproducing the proﬁle shape. c) We idealize the eigenstrain proﬁle and represent the plate as a bilayer consisting of two layers of the thicknessh=2. The eigenstrain "tand"bassigned to each of the layers is derived from the real eigenstrain proﬁle by equating the total eigenstrain and the ﬁrst eigenstrain moment. d) Both real and idealized eigenstrain proﬁles lead to the same deformed shape, which is bent and stretched with respect to the initial state. 12 Journal of Manufacturing Science and Engineering words, it describes the local stretching of the surface. The 22matrix containing coefﬁcients of the ﬁrst fundamental form a(x;y)is computed as: a(x;y) =2 64@x~ m@x~ m @ x~ m@y~ m @y~ m@x~ m @ y~ m@y~ m3 75; (7) where@x~ m=@~ m=@x and@y~ m=@~ m=@y are two vectors tangent to the mid-surface at the point ~ m(x;y). If a certain area of the surface does not undergo any stretching, the ﬁrst fundamental form in this area is represented by the identity matrix I. Together with the ﬁrst fundamental form, the second fundamental form determines local cur- vatures on a surface. The matrix containing its coefﬁcients b(x;y)is computed as: b(x;y) =2 64@xx~ m~ n @ xy~ m~ n @yx~ m~ n @ yy~ m~ n3 75=2 64@x~ m@x~ n @ x~ m@y~ n @y~ m@x~ n @ y~ m@y~ n3 75; (8) where@xx~ m,@xy~ mand@yy~ mdenote the second derivatives of ~ m(x;y). If a surface is locally ﬂat, its second fundamental form at this area is described by the zero matrix, because all second derivatives of the mapping ~ mare orthogonal to the normal vector ~ n. The two fundamental forms deﬁne a unique surface up to solid body motions. The surface is integral if its fundamental forms are compatible , i.e., if they satisfy three partial differential equa- tions called the Gauss–Peterson–Mainardi–Codazzi (GPMC) equations [25]. Hence, the two ﬁnal fundamental forms of a non-Euclidean plate are compatible, while its restfundamental forms are not. This phenomenon illustrated in Fig. 3 is also called geometric incompatibility. 13 Journal of Manufacturing Science and Engineering Fig. 3. An example of geometric incompatibility inspired by [14]. The initial conﬁguration described by the fundamental forms (ainit;binit)is a ﬂat unstretched disc, so that ainit is the identity matrix and binit is the zero matrix along the whole sur- face. The rest conﬁguration (ar;br)prescribes isotropic in-plane stretching of the central part (green) while conserving the ﬂat shape, so that ar6=ainit andbr=binit. In the general case, brcan be different from binit as well. The rest fundamental forms are incompatible, because the adoption of the rest conﬁguration means superposition of the inner (green) part and the outer (blue) part of the disc and thus loss of integrity. In order to conserve its integrity, the disc adopts a curved ﬁnal conﬁguration (af;bf) described by compatible fundamental forms. However, the disc stays residually stressed in its ﬁnal conﬁguration. Elastic energy and the forward problem resolution In the terms of the theory of non-Euclidean plates, the forward problem consists in determining the ﬁnal conﬁguration as a function of the rest conﬁguration. First, we consider the simple case of an initially ﬂat monolayer plate, which is subjected to a rest conﬁguration described by the incompatible fundamental forms arandbr. The rest and ﬁnal conﬁgurations of such a plate are related through the elastic energy functional. We denote the ﬁnal fundamental forms as afandbf and we express the elastic energy EMLof an integral monolayer plate as [19]: EML=1 2Z Uh 4 a1 rafI 2 e+h3 12 a1 r(bfbr) 2 ep detardxdy: (9) In this expression we have introduced the elastic norm kAk2 e= Tr2(A) + 2 Tr(A2)with coefﬁ- cients =Y=(12)and =Y=(2 + 2). Here,Yis the Y oung modulus and is the Poisson’s 14 Journal of Manufacturing Science and Engineering ratio. A zero elastic energy means that the ﬁnal conﬁguration perfectly coincides with the rest one. If the rest fundamental forms are incompatible, an integral plate adopts an equilibrated ﬁnal con- ﬁguration that minimizes the elastic energy under constraints that afandbfbe compatible. The ﬁrst term of the elastic energy functional represents the stretching energy, and the second term deﬁnes the bending energy. If the plate is thin, the bending term is small in comparison to the stretching one, so the plate stretches as prescribed by the rest conﬁguration but adopts a different curvature. In other words, afin this case is close to ar, butbfcan be largely different from br. On the contrary, a thick plate adopts the curvature prescribed by the rest conﬁguration but stretches in a different way, so that bfis close to br[12]. Now let us consider a bilayer plate, where each layer is of thickness h=2. We suppose that the plate is initially ﬂat, and that its layers exhibit nonuniform in-plane swelling or shrinking. The eigenstrain introduced in each layer may vary along the surface, but it is constant along the layer thickness. Locally, the eigenstrain is different for each layer, so that the rest ﬁrst fundamental forms of each layer ( ar;tandar;b) are different. Matrices ar;tandar;bcontain information on the principal eigenstrain direction and magnitude on the top and bottom layers respectively. Essentially, the ar;tandar;bdescribe stretching that the layers would adopt if they were not attached together. By assuming that each layer expands uniformly across its thickness, all terms of the rest second fundamental form of each layer are zero along the whole surface: br;t=br;b=2 640 0 0 03 75: (10) Accordingly, ar;tandar;bfully describe the rest conﬁguration. The forward problem consists in ﬁnding afandbfthat describe the shape of the integral plate mid-surface after reconﬁguration. The reconﬁguration process for the bilayer case is presented in Fig. 4. Following [19], we express the elastic energy of the bilayer plate as the sum of the elastic energies of two monolayers of 15 Journal of Manufacturing Science and Engineering Inial con gura on Rest con gura on Final con gura on In-plane stretching of the upper layerAdoption of a spherical shape ( b t) ( ) () Fig. 4. An example of geometric incompatibility in a bilayer plate inspired by [19]. The initial conﬁguration described by the funda- mental forms (ainit;binit)is a ﬂat unstretched rectangular plate, so that ainit is the identity matrix and binit is the zero matrix along the whole surface. The rest conﬁguration (ar;t;ar;b)prescribes isotropic in-plane stretching of the top layer and shrinking of the bottom layer while conserving the ﬂat shape, so that ar;t6=ainit,ar;b6=ainit,br;t=binit andbr;b=binit. Since peen forming causes in-plane eigenstrain we do not consider cases when br;t6=0andbr;b6=0for our simulations, so that the rest conﬁguration is entirely described by (ar;t;ar;b). Adoption of the rest conﬁguration means dissection of the plate in two layers. Instead, the plate adopts an integral but residually stressed ﬁnal conﬁguration (af;bf)described by compatible fundamental forms. thicknessh=2. After integration over the total plate thickness we obtain: EBL=1 2Z Uh 8 a1 r;bafI 2 e+h3 24 a1 r;bbf 2 e+h2 8D a1 r;bafI ;a1 r;bbfE ep detar;bdxdy +1 2Z Uh 8 a1 r;tafI 2 e+h3 24 a1 r;tbf 2 eh2 8  a1 r;tafI ;a1 r;tbf ep detar;tdxdy:(11) The elastic energy inner product h;ieintroduced in this context deﬁnes the following operation: hA;Bie= Tr(A)Tr(B)+2 Tr(AB). Similarly to the monolayer case, the plate adopts a curved ﬁnal conﬁguration (af;bf)that minimizes the elastic energy. Generally, the ﬁnal conﬁguration is not unique, and moreover, the plate can get stuck in a conﬁguration corresponding to a local energetic minimum on its way to the global minimum. For a bilayer subjected to any rest conﬁguration (ar;b;ar;t), there exists an equivalent mono- layer that morphs into the same ﬁnal conﬁguration (af;bf)after being subjected to a rest conﬁg- uration (ar;br). The relation between (ar;b;ar;t)and(ar;br)is derived by equating the monolayer energy (Eqn. 9) and the bilayer energy (Eqn. 11) and is expressed as [19]: 16 Journal of Manufacturing Science and Engineering 8 >< >:ar=1 2(ar;b+ar;t) ; br=3 4h(ar;bar;t):(12) Inversely, 8 >< >:ar;t=ar2h 3br; ar;b=ar+2h 3br:(13) Here, the monolayer and bilayer plates are supposed to have the same initial geometry and the same total thickness h. The inverse problem resolution The inverse problem consists in determining the rest conﬁguration that leads to a target con- ﬁguration due to the elastic material response. In the bilayer case, it means ﬁnding the rest ﬁrst fundamental forms of the bottom and top layers ar;bandar;t, respectively as a function of the target shape described by atarandbtar. As each fundamental form is represented by a symmet- ric22matrix, the target conﬁguration is locally deﬁned by six scalar variables. At the same time, the rest conﬁguration also has six degrees of freedom, so the problem admits an analyti- cal solution provided by Eqn. (13). Moreover, the ﬁnal conﬁguration obtained with this solution is stress-free [19]. A practical application of the analytically calculated solution means inducing local orthotropic eigenstrain, so that "11(x;y;z )6="22(x;y;z ), where"11and"22are the local principal eigenstrains. In other words, one has to locally control the principal strain directions and the strain magnitude along both directions, and this – at each of the two layers separately. Such local control over the direction of expansions is generally not possible with shot peen forming, because here we assume that shot peening induces local isotropic in-plane strain. Thus, at each point we only 17 Journal of Manufacturing Science and Engineering control two degrees of freedom: "t(x;y)and"b(x;y). Hence, we have control over fewer degrees of freedom than input variables. In this case, a solution leading exactly to the target shape may not exist, so we can only numerically optimize ar;tandar;b. The uniqueness of solution is not guaranteed either. METHODOLOGY Our method for the inverse problem resolution consists in the iterative correction of the peening pattern on a local scale until convergence is reached within a tolerance and subsequent grouping of the pattern. The grouping algorithm divides the pattern in zones treated with predeﬁned peening regimes. Iterative adjustment of the peening pattern The iterative method idea is to adjust the rest conﬁguration by comparing the local stretching and curvature of the current shape with the stretching and curvature of the target shape. The current shape is computed at each iteration through a numerical resolution of the forward problem. For the initial guess, we use analytical expressions to approximately deﬁne an appropriate rest conﬁguration. The initial guess We characterize the target conﬁguration in terms of fundamental forms (atar;btar). Let us consider the monolayer rest conﬁguration aortho r;bortho r described by aortho r =atarandbortho r = btar, that in a general case prescribes local orthotropic in-plane strain. According to the expression for the elastic energy (Eqn. 9), imposition of this rest conﬁguration makes the plate adopt exactly the target conﬁguration, because the elastic energy in this case equals zero. Following [19], we express the equivalent bilayer rest conﬁguration aortho r;t;aortho r;b using Eqn. (13): 18 Journal of Manufacturing Science and Engineering 8 >< >:aortho r;t =aortho r2h 3bortho r; aortho r;b =aortho r +2h 3bortho r:(14) Equivalently, using the deﬁnition of aortho r;bortho r we rewrite: 8 >< >:aortho r;t =atar2h 3btar; aortho r;b =atar+2h 3btar:(15) Thus, application of the rest conﬁguration aortho r;t;aortho r;b leads to the target shape (atar;btar). However, this conﬁguration implies local orthotropic eigenstrain, which is not feasible with shot peen forming. We comply with this constraint and ﬁnd a suitable local isotropic eigenstrain based on this prediction. To that end, we ﬁrst compute the local eigenstrains in the principal directions in the top ("t 11,"t 22)and bottom ("b 11,"b 22)layers of the bilayer prescribed by aortho r;t andaortho r;b, respectively. Next, we take their averages "t avgand"b avgand impose them locally in all in-plane directions thus making the initial guess. To ﬁnd ("t 11,"t 22)and("b 11,"b 22)we perform a spectral decomposition of aortho r;t andaortho r;b, re- spectively [19]. At the top layer, the distortions prescribed by aortho r;t imply stretching by a factor of ("t 11+1)in the ﬁrst principal direction and by a factor of ("t 22+1)in the orthogonal second principal direction. The ﬁrst principal direction is rotated by an angle of twith respect to the x-axis. At the bottom layer, the stretch factors are ("b 11+ 1) and("b 22+ 1), and the ﬁrst principal direction is rotated by an angle of bwith respect to the x-axis. The initial conﬁguration is unstretched, so its ﬁrst fundamental form ainitis represented by the identity matrix: 19 Journal of Manufacturing Science and Engineering ainit=I: (16) Consequently, the eigenvalues of aortho t equal ("t 11+ 1)2and("t 22+ 1)2, and the eigenvalues of aortho bequal ("b 11+ 1)2and("b 22+ 1)2, so that: aortho r;j =2 64cos(j)sin(j) sin(j) cos(j)3 75T2 64("j 11+ 1)20 0 ("j 22+ 1)23 752 64cos(j)sin(j) sin(j) cos(j)3 75;forj=t;b: (17) We deduce the local eigenstrain in the principal directions ("t 11,"t 22)and("b 11,"b 22)from the eigen- values and calculate the average local eigenstrain "t avgand"b avgfor both layers: "j avg="j 11+"j 22 2;forj=t;b: (18) We impose the local isotropic eigenstrain "t avgand"b avgon the top and bottom layers, respectively, thus making the initial guess. The bilayer rest fundamental forms ar;tandar;bcorresponding to this strain are expressed as: ar;j=2 64("j avg+ 1)20 0 ("j avg+ 1)23 75;forj=t;b: (19) We substitute (ar;t;ar;b)to the bilayer elastic energy functional (Eqn. 11) and ﬁnd the current shape 20 Journal of Manufacturing Science and Engineering (ac;bc)that minimizes the elastic energy. Following Eqn. (12), it is possible to ﬁnd equivalent monolayer rest fundamental forms (ar;br). Prescription of the rest fundamental forms (ar;br)to a monolayer plate leads to the same cur- rent shape (ac;bc)as the prescription of (ar;t;ar;b)to a bilayer plate. Essentially, (ar;br)can be viewed as the rest conﬁguration imposed on the bilayer mid-surface. The forms arandbrmay be incompatible as they were deﬁned analytically, but acandbcare always compatible because they describe a surface in Euclidean space. Consequently, (ac;bc)are typically different from (ar;br). Adjustment of the local stretching In the previous subsection, we deﬁned a procedure to obtain a quick estimate of the inverse problem by solving for orthotropic expansions and averaging them to estimate isotropic expan- sions. In this subsection and the next, we seek to iteratively correct these expansions by using only local information. As the rest conﬁguration (ar;t;ar;b)prescribes local isotropic in-plane strain, we measure stretching in terms of local areas. The area Aof each region of the surface constrained by (x;y)2U1R2is expressed in terms of the ﬁrst fundamental form aasA=s U1p det(a)dxdy . The ﬁrst fundamental form is considered constant inside small regions, so we conclude that the current area of each small region Acand its target area Atarare related as: Atar Acp det (atar)p det (ac)=kA: (20) This means that if we locally multiply acby the coefﬁcient kA, then the current area will equal that of the target. However, we are only able to inﬂuence acindirectly through the adjustment of ar. Consequently, as a part of the iterative procedure, we multiply arbykAand thus obtain the equivalent monolayer rest fundamental form anew rto be imposed during the subsequent iteration: 21 Journal of Manufacturing Science and Engineering anew r=kAar: (21) Theanew rmay be different from the current ﬁrst fundamental form on the subsequent iteration anew c, and thus this correction of the rest fundamental form does not lead to an exact solution but allows to approach it. In other words, multiplication of arbykAdoes not correct the local area exactly by the coefﬁcient kAbut reduces the difference between the current local stretching and the target one. Adjustment of the local curvature We characterize the surface curvature in terms of the local mean curvature H. By deﬁnition, His the average of two local principal curvatures 1and2, that are computed as eigenvalues of the shape operator S=a1b[26]. We compute the local ratios kHbetween the current mean curvatures Hcand the target mean curvatures Htarand assign an upper threshold forjkHj: 8 >>< >>:kH=Htar Hcfor Htar Hc <; kH=sgnHtar Hc for Htar Hc :(22) Here,Htaris the average of the two eigenvalues of Star=a1 tarbtar, andHcis the average of the two eigenvalues of Sc=a1 cbc. The threshold is assigned in order to deal with special cases when jHcjis small. Provided that multiplication of a matrix by a constant multiplies its eigenvalues by the same constant, the multiplication of ScbykHwould make the current local mean curvatures equal to those of the target. As we are unable to adjust any of the current fundamental forms (ac;bc) directly, we inﬂuence them through adjustment of the rest fundamental forms (ar;br)in order to get: 22 Journal of Manufacturing Science and Engineering Snew r=kHSr=kHa1 rbr=kHkAa1 r kAbr=kHkA(anew r)1br: (23) Thus, we deﬁne: bnew r=kHkAbr: (24) Once the anew randbnew rare found, we compute the bilayer rest fundamental forms (anew r;t;anew r;b) as: 8 >< >:anew r;t=anew r2h 3bnew r; anew r;b=anew r+2h 3bnew r:(25) Next, we substitute (anew r;t;anew r;b)to the bilayer elastic energy functional (Eqn. 11) and ﬁnd the current shape (anew c;bnew c)that minimizes the elastic energy. Subsequent iterations and stop criterion We compare the current shape with the target shape and recalculate the bilayer rest funda- mental forms until a convergence criterion is satisﬁed. The convergence criterion is based on the calculation of the Hausdorff distance dHbetween the new current conﬁguration Cnew cdeﬁned by (anew c;bnew c)and the current conﬁguration from the previous iteration Ccdeﬁned by (ac;bc). We nondimensionalize dHby the square root of the total area of the plate in its initial conﬁguration Atotal. We stop iterating either when a predeﬁned maximal number of iterations Mis reached, or 23 Journal of Manufacturing Science and Engineering when the nondimensionalized Hausdorff distance becomes inferior to a chosen threshold : dH(Cc;Cnew c)pAtotal< : (26) Calculation of the adjusted eigenstrain We denote the bilayer ﬁnal rest fundamental forms, i.e., the ones obtained on the last iteration, as(arf;t;arf;b). To relate them with peen forming parameters, we determine the recalculated eigenstrain ("rf;t;"rf;b). Provided that the iterative adjustment implies only addition, subtraction and multiplication by a constant of diagonal matrices, the arf;tandarf;bare diagonal. Moreover, the imposed local eigenstrain ("rf;t;"rf;b)is isotropic, so arf;tandarf;bhave the following form: arf;j=2 64 "rf;j+ 120 0 "rf;j+ 123 75;forj=t;b: (27) Consequently, "rf;j=q a11 rf;j1;forj=t;b: (28) Numerical implementation We mesh the plate mid-surface with triangular elements and follow the energy calculation strat- egy presented in [17] [27]. The ﬁrst and the second fundamental forms are estimated separately for each element and are constant inside the element. The global elastic energy is calculated as 24 Journal of Manufacturing Science and Engineering a sum of local energetical contributions from all the elements. The ﬁrst fundamental form on a triangular element such as that schematized in Fig. 5 depends entirely on the coordinates of the vertices. The three vertices are deﬁned by position vectors ~ v0,~ v1 and~ v2, and the edge vectors constituting the triangle are expressed as: ~ e0=~ v1~ v0,~ e1=~ v2~ v1, ~ e2=~ v0~ v2. These vectors are tangent to the plane containing the triangle, so following Eqn. (7), the ﬁrst fundamental form on a triangular element is computed as: a=2 64~ e1~ e1~ e1~ e2 ~ e2~ e1~ e2~ e23 75: (29) In a general case of isotropic eigenstrain "rimposed on a triangular element, its rest ﬁrst funda- mental form aris expressed as: ar=2 64~ einit 1~ einit 1~ einit 1~ einit 2 ~ einit 2~ einit 1~ einit 2~ einit 23 752 64("r+ 1)20 0 ("r+ 1)23 75=ainit2 64("r+ 1)20 0 ("r+ 1)23 75; (30) where~ einit 1,~ einit 2,~ einit 3are the edge vectors in the initial conﬁguration and ainitis the element ﬁrst fundamental form in its initial conﬁguration. Otherwise, if an orthotropic in-plane strain is imposed, then the rest ﬁrst fundamental form aortho r becomes: aortho r =ainit2 64cos()sin() sin() cos()3 75T2 64("r 11+ 1)20 0 ("r 22+ 1)23 752 64cos()sin() sin() cos()3 75; (31) where"r 11and"r 22deﬁne the eigenstrain imposed in the principal directions and stands for the 25 Journal of Manufacturing Science and Engineering Fig. 5. A triangular mesh element and vectors that determine local fundamental forms: the vertex position vectors ( ~ v0,~ v1and~ v2), the edge vectors ( ~ e0,~ e1and~ e2) and the mid-edge normals ( ~ n0,~ n1and~ n2). The mid-edge normals are perpendicular to the edge, and their direction is determined by the angle of inclination with respect to the average normal of the two adjacent faces. Thus, the direction of~ n2is determined by the angle '2. This angle is measured with respect to ~ navg, which is the average normal of the two faces that share the edge ~ e2. angle between the ﬁrst principal direction and the xaxis. The second fundamental form deﬁnes the surface curvature, so the information about surface normals is required. In this connection we introduce a unit normal vector ~ ni; i= 1;2;3at the center of each edge of the mesh (the edge-director ). This vector is normal to the edge, and its angle of inclination 'i; i= 1;2;3in the plane perpendicular to the edge provides a supplementary degree of freedom. This angle is measured with respect to the average of the adjacent face normals [27]. A ﬁnite-difference approximation of the derivatives appearing in Eqn. (8) yields the following expression for the second fundamental form of a triangular element: b=2 64~ e12(~ n0~ n2)~ e12(~ n1~ n0) ~ e12(~ n1~ n0)~ e22(~ n1~ n0)3 75=2 64~ e12(~ n0~ n2)~ e1~ n0 ~ e1~ n0~ e22(~ n1~ n0)3 75: (32) According to Eqn. (11), the global elastic energy for a plate composed of Ktriangular elements is expressed in terms of local fundamental forms as: 26 Journal of Manufacturing Science and Engineering EBL=1 2KX k=1hk 8 (a1 r;b)k(af)kI 2 e+(hk)3 24 (a1 r;b)k(bf)k 2 e +(hk)2 8D (a1 r;b)k(af)kI ;(a1 r;b)k(bf)kE eq det (ar;b)k +1 2KX k=1hk 8 (a1 r;t)k(af)kI 2 e+(hk)3 24 (a1 r;t)k(bf)k 2 e (hk)2 8D (a1 r;t)k(af)kI ;(a1 r;t)k(bf)kE eq det (ar;t)k:(33) The local plate thickness hkcan be different for each element. Similarly, the local Y oung’s modulus Ykand the Poisson’s ratio kcan vary along the plate. As follows from Eqn. (30), the element’s rest ﬁrst fundamental forms ar;tandar;bfor the isotropic growth case are deﬁned by the imposed local eigenstrain "tand"band the initial ver- tex positions ~ vinit i; i= 1;2;3. At the same time, according to Eqn. (29) and Eqn. (32), the two ﬁnal fundamental forms afandbfare deﬁned by the ﬁnal vertex positions ~ vf i; i= 1;2;3;and ﬁnal angles of inclination of the edge directors 'f i; i= 1;2;3. Numerically, the forward problem consists in minimizing the global elastic energy functional (Eqn. 33) with respect to ~ vf iand'f iprovided with ~ vinit i,"tand"b. We perform the minimization using a quasi-Newton L-BFGS algorithm [28]. The gradients of the elastic energy functional required by the minimization algorithm are computed an- alytically following [29]. The corresponding programming code for the forward problem resolution was developed by van Rees et al. and is publicly accessible [30]. In this implementation, the inverse problem resolution consists in ﬁnding the local eigenstrain "rf;tand"rf;bto be imposed on each triangular element. This means that the iterative correction of the rest fundamental forms ar;tandar;bis executed on a local scale for each element separately, whilst the forward problem is resolved on each iteration on a global scale, thus reﬂecting the mechanics of the plate. 27 Journal of Manufacturing Science and Engineering Algorithm 1 The inverse problem resolution The iterative loop 1:while condition (26) is not satisﬁed and the number of iterations is below maximum do 2: foreach triangular element do 3: Find principal curvatures (c 1;c 2)as eigenvalues of the shape operator Sc=a1 cbc 4: Compute the current mean curvature as Hc= 0:5 (c 1+c 2) 5: Compute the current area as Ac=p det (ac) 6: Compute the ratios kAandkHfollowing Eqn. (20) and Eqn. (22) respectively 7: Compute the monolayer rest fundamental forms arandbrfollowing Eqn. (12) using only local information 8: Compute the adjusted monolayer rest fundamental forms anew randbnew rfollowing Eqn. (21) and Eqn. (24) respectively 9: Compute the adjusted bilayer rest fundamental forms anew r;tandanew r;bfollowing Eqn. (25) 10: end for 11: Substitute anew r;tandanew r;bto the bilayer elastic energy functional (Eqn. 33) and minimize it to solve the forward problem and ﬁnd the current conﬁguration (anew c;bnew c) 12:end while The initial guess 1:foreach triangular element do 2: Compute the ﬁrst fundamental form of the initial shape ainitfollowing Eqn. (29) 3: Compute the monolayer fundamental forms of the target shape atarandbtarfollowing Eqn. (29) and Eqn. (32) respectively 4: Find principal curvatures tar 1andtar 2as eigenvalues of the target shape operator Star= a1 tarbtar 5: Compute the target mean curvature as Htar= 0:5(tar 1+tar 2) 6: Compute the target area as Atar=p det (atar) 7: Compute the orthotropic bilayer rest fundamental forms aortho r;t andaortho r;bfollowing Eqn. (15) 8: Find eigenvalues j 1andj 2of(a1 initaortho r;j)forj=t;b 9: Compute the orthotropic eigenstrain as "j ii=q j i1fori= 1;2andj=t;b 10: Compute the average eigenstrain "t avgand"b avgfollowing Eqn. (18) 11: Compute bilayer rest fundamental forms ar;tandar;bby substituting "t avgand"b avgrespec- tively for"rin Eqn. (30) 12:end for 13:Substitute ar;tandar;bto the bilayer elastic energy functional (Eqn. 33) and minimize it to solve the forward problem and ﬁnd the current conﬁguration (ac;bc) 28 Journal of Manufacturing Science and Engineering Algorithm 1 The inverse problem resolution Final step after exiting the iterative loop 1:foreach triangular element do 2: Compute the eigenstrain "rf;tand"rf;bprescribed by arf;t=anew r;tand by arf;b=anew r;b respectively following Eqn. (28) 3:end for Grouping of the peening pattern In the general case, the eigenstrains ("rf;t;"rf;b)provided by Alg. 1 are different for each element and can take any real values. Peen forming often deals with smoothly curved target shapes, so ("rf;t;"rf;b)may also vary smoothly along the surface given that these two values depend on the target shape curvature. We call the eigenstrain pattern provided by the Alg. 1 the free pattern. From a practical point of view, each pair ("rf;t;"rf;b)represents a peening regime. However, a limited number of regimes is available when peening a real part. Thus, we divide the pattern into zones with uniform prescribed eigenstrain and obtain a grouped pattern. We associate all triangular elements with points on a plane with cartesian coordinates "t;"b , and the coordinates of each point kare determined by the eigenstrain "rf;t k;"rf;b k assigned to the corresponding element k, as illustrated in Fig. 6. We divide the points in groups, and the group centroids are determined by the predeﬁned peening regimes. We denote the centroid of a group nby "cen;t n;"cen;b n . Each point "rf;t k;"rf;b k is attributed to the group with the closest centroid in terms of Euclidean distance. When all points are grouped, we homogenize the eigenstrain inside each group, i.e., we assign the eigenstrain "cen;t n;"cen;b n to all triangular elements that fall into the groupn. Consider N predeﬁned peening regimes such as the ones illustrated in Fig. 7. Regime i= 1;2;:::;N induces expansions "t iand"b ion the treated and opposite layers, respectively. In ad- dition, we consider the lack of treatment "t 0="b 0= 0. Since the top and bottom surfaces can be peened independently, there are (N+ 1)2possible treatment combinations. Each treatment combination gives rise to a group centroid. Combining regime i= 1;2;:::;N on the top surface with regime j= 1;2;:::;N on the bottom surface leads to the following expansions of the top and 29 Journal of Manufacturing Science and Engineering    b Fig. 6. Graphical representation of grouping of the eigenstrain pattern. The triangulated ﬂat initial conﬁguration (a) and the target conﬁguration - wavy shape (b) - are the input data for the inverse problem resolution. We divide the eigenstrain pattern (c) into three zones treated uniformly (d). A plane with cartesian coordinates ("t;"b)(e) illustrates the grouping from the numerical point of view. The colored points correspond to the eigenstrain assigned to each element of the triangular mesh. The four centroids are denoted by colored squares. They are generated by one peening regime and the lack of treatment as an additional regime. The grouping is based on calculation of the least Euclidean distance from the points to the group centroids. In the presented case there are no points close to the centroid denoting treatment from both sides, so the corresponding group is empty. Once the points are divided in groups, we homogenize the eigenstrain for all the elements attributed to the same group.             Fig. 7. Calculation of the group centroids. The group ijimplies treatment from the bottom side with regime iand treatment from the top side with regime j. Its centroid is determined by a pair of parameters "cen;t ijand"cen;b ij, and each of them is a superposition of eigenstrain generated by the two regimes that form this group: "cen;t ij="t i+"t j,"cen;b ij="b i+"b j. In this example, the regime iis more intense than the regime j, so"cen;t ij<"cen;b ij. bottom surface: ("t i;"b i)+("t j;"b j) = ("cen;t ij;"cen;b ij). Figure 7 illustrates this principle. Figure 6 e) also provides a cartesian representation of (N+ 1)2centroids for the case N= 1. Numerical validation of the inverse problem solver We generated target shapes numerically to test our algorithms for the iterative inverse problem resolution and grouping. To ensure that the target shapes were achievable with peen forming, we generated them by assigning a random peening pattern to the initial conﬁguration and then solved the forward problem. The random peening patterns were generated following Alg. 2. When the target shapes were generated, we solved the inverse problem for each of them 30 Journal of Manufacturing Science and Engineering Algorithm 2 Generation of random peening patterns 1:Mark 1 to 6 random points on the top and bottom surfaces of the plate 2:foreach point do 3: Draw a square of random size (but not bigger than the plate size) centred on the point 4: Assign randomly one of the available peening regimes to the square 5: ifthe square protrudes beyond the plate area then 6: Translate the part that protrudes symmetrically on the other side of the plate 7: end if 8: ifthe square superimposes with a previously drawn square on the same side then 9: Erase the previously assigned regime in the superimposing area and leave only the latest one 10: end if 11:end for following Alg. 1 and then grouped the peening pattern. The predeﬁned regimes were ﬁxed as those that were originally used to generate the target shapes. To quantify the error, we solved the forward problem for the free and the grouped patterns. We thus obtained two ﬁnal shapes for each test case and compared them with the target shape by calculating the nondimensionalized Hausdorff distance : =dH(Cf;Ctar)pAtotal; (34) where Cfstands for the ﬁnal conﬁguration and Ctarstands for the target conﬁguration. The overall process for the inverse problem validation is schematized in Fig. 8. RESULTS We generated 200 random patterns with Alg. 2 and applied them on a ﬂat square plate ( 11m). The plate thickness was arbitrarily assigned in each case and ranged from 2 mm to 15 mm. The Poisson’s ratio was also arbitrarily picked between 0.32 and 0.36. The plate was meshed with 1152 triangular elements. The forward problem resolution took 10-30 seconds for one shape, depending 31 Journal of Manufacturing Science and Engineering Initial shape Generation of a random peening pattern (Algorithm 2) Random pattern + Random thickness + Random Poisson's ratioThe forward problem resolution Target shape Free pattern Final shape #1Pattern discretization Grouped pattern Final shape #2 Calculation of the Hausdorff distanceThe inverse problem resolution (Algorithm 1) The forward problem resolutionThe forward problem resolution Fig. 8. A validation process to test the quality of the inverse problem resolution. The initial shape is ﬁxed as a ﬂat 1 1 m plate, its thickness and Poisson’s ratio are chosen randomly. These parameters are kept constant throughout the whole validation process. A random peening pattern is generated following Alg. 2 and assigned to the initial shape. Different colors on the peening pattern correspond to different peening regimes applied from the top, bottom or both sides. We solve the forward problem taking the random pattern as input, and the result is used as the target shape for the inverse problem validation. The eigenstrain pattern is determined following Alg. 1 and grouped using pre-determined regimes. Finally, the forward problem is resolved for the free and the grouped patterns. The difference between the target shape and the two ﬁnal shapes is quantiﬁed with the nondimensionalized Hausdorff distance . on the pattern and the plate thickness. Thus, the forward problem resolution took longer time for thinner plates with bigger treated areas due to larger deﬂection of these plates. For the ﬁrst 100 test cases ( series 1 ), we made only one peening regime available, so the entire treated area was peened with the same parameters. For the second 100 test cases ( series 2), we assigned randomly one of four available peening regimes to each square on both sides. We considered real peening regimes presented in [7]. The authors of this paper deduced idealized eigenstrain proﬁles from the residual stress measurements performed on the treated specimens. The idealized eigenstrain proﬁles were formulated as one uniformly expanding layer of a constant thickness. The layer thickness and the eigenstrain magnitude were different for each regime. We reformulated the idealized eigenstrain proﬁles in terms of "t;"b by equating the total eigenstrain and the ﬁrst eigenstrain moment 1induced by "t;"b and by the one expanding layer. Table 1 summarizes the eigenstrain "t;"b induced by each of the four regimes applied on a 5 mm thick plate from the top side. Figures 9 and 10 present the free and the grouped patterns along with the convergence curves 32 Journal of Manufacturing Science and Engineering Table 1. The in-plane eigenstrain induced by the regimes used to generate the random peening patterns. The in-plane eigenstrain is presented for the case of a 5 mm thick plate. The peening regimes represent four real treatments inducing different eigenstrain proﬁles, which are examined in [7]. Regime"t103"b103 1 2.5 -0.4 2 3.2 -0.6 3 1.7 -0.4 4 2.0 -0.4 for two particular test cases from series 2: a low-error case and a high-error case. Figures 9 b) and 10 b) show that the free pattern on the ﬁnal iteration is locally close to the originally generated random pattern. Due to that, most of the elements are grouped correctly, so that the eigenstrains prescribed by the grouped and the random patterns to these elements become equal (Fig. 9 c and Fig. 10 c). However, in each case there are elements that are attributed to a wrong group. In the low-error case this happens only for several elements. Consequently, the dimensionless error is lower for the grouped pattern than for the free pattern. In the high-error case the grouped pattern undergoes the checkerboard problem, meaning that the pattern locally alternates two peening regimes over a certain area (Fig. 10 c). A large area affected by the checkerboard problem increases the in comparison with the free pattern (Fig. 10 e). Nevertheless, the regimes that are mixed up in the checkerboard-affected zone have only a slight difference in terms of the induced eigenstrain, so the increases up to 0:23% at most. In terms of convergence, the most important correction is done on the ﬁrst iteration after the initial guess, as illustrated by Fig. 9 e) and Fig. 10 e). Thus, the ﬁrst iteration decreases the by 65%on average, and the becomes inferior to 0:1%for all the cases. All the subsequent iterations together decrease the obtained on the ﬁrst iteration by 60% on average. Because of the local nature of the eigenstrain adjustment, the solution does not converge to the exact target shape but to a shape which is close to the target. In other words, after several iterations the plateaus at a low but ﬁnite level (Fig. 9 e and Fig. 10 e). The may slightly grow during the subsequent iterations, but the stop criterion (26) terminates the iterative process as soon as this happens. 33 Journal of Manufacturing Science and Engineering Random pattern Target shapeFree pattern a) d)b) e)c)Low-error case 0 1 2 3 4 5 6 7 8 9 1000.20.40.60.81Grouped pattern Ω =0.08% Ω =0.02% Free pattern Dimensionless error (%) IterationsGrouped pattern Fig. 9. One of the test cases from series 2 with low dimensionless error. The plate is 7 mm thick, and the Poisson’s ratio equals 0.34. a) The random peening pattern generated for this test case (top layer). Different colors on the peening pattern denote different peening regimes. b) The free peening pattern on the ﬁnal iteration. c) The grouped peening pattern. A visual comparison of (c) and (a) shows that almost all elements were attributed to a correct group. Consequently, grouping decreased the dimensionless error . d) The target shape induced by the pattern (a). The deformations are at their original scale. e) The convergence curve showing the dimensionless error on each iteration. The optimization required 10 iterations, that was the maximum ﬁxed for this test. The results of the numerical validation are presented with histograms in Fig. 11. The inverse problem solver provided free peening patterns that led to the target shape with the inferior to 0:35% for both series of tests (Fig. 11, top). The pattern optimization needed between 2 to 10 iterations, depending on the target shape. The histograms in Fig. 11 show that the pattern grouping has decreased the in most of the cases. Thus, the original pattern was perfectly reproduced for all test cases in series 1, so the after grouping became less than 103%. The corresponding for series 2 was bigger because of the higher complexity of the grouping problem: four available regimes induced twenty-ﬁve group centroids for series 2, while there were only four centroids induced by one available regime for 34 Journal of Manufacturing Science and Engineering 0 1 2 300.511.522.53Ω=0.09% Ω =0.23% Random pattern Target shapeFree pattern a) d)b) e)c)High-error case Grouped pattern Free pattern Dimensionless error (%) IterationsGrouped pattern Fig. 10. One of the test cases from series 2 with high dimensionless error. The plate is 4 mm thick, and the Poisson’s ratio equals 0.34. a) The random peening pattern generated for this test case (top layer). Different colors on the peening pattern denote different peening regimes. b) The free peening pattern on the ﬁnal iteration. c) The grouped peening pattern. A visual comparison of (c) and (a) indicates elements that were not attributed to a correct group. The deformations are at their original scale. However, grouping decreased the dimensionless error in this case. d) The target shape induced by the pattern (a). e) The convergence curve showing the dimensionless error on each iteration. The optimization required 3 iterations to converge. series 1. The group centroids for series 2 were situated close to each other, so several test cases were signiﬁcantly affected by the checkerboard problem, as illustrates Fig. 10 c). This explains the increased error after grouping for 18 cases out of 100 from series 2. DISCUSSION The described inverse problem resolution algorithm relies on the assumption that the peening treatment induces isotropic expansions. In practice, however, peen forming sometimes induces different eigenstrain along the xandyaxes:"xx(x;y;z )6="yy(x;y;z ). Such anisotropic expan- sions are due to plastic anisotropy of the treated material, which is especially explicit for rolled 35 Journal of Manufacturing Science and Engineering Fig. 11. Histograms evaluating the dimensionless error between the target shapes and the ﬁnal shapes obtained during the numerical validation of the inverse problem solver. All in all, 200 test cases were considered. They were divided in two series of 100 cases each. The target shapes in series 1 were obtained with one peening regime, and the same peening regime was ﬁxed as the only available for the grouping stage. The target shapes in series 2 were obtained with four different peening regimes, and the four predeﬁned regimes were available for grouping. Figures (a) and (b) show the dimensionless error induced by application of the free pattern for series 1 and 2 correspondingly. Similarly, ﬁgures (c) and (d) represent the dimensionless error induced by the grouped pattern. aluminum sheets, and to prestressing the component in one direction before treatment, i.e., stress peen forming. This effect is examined in detail in [31]. For uniform plastic anisotropy which does not vary over the area of the plate, the inverse problem resolution algorithm can be easily adapted by introducing a ﬁxed eigenstrain anisotropy coefﬁcient in the model: "j xx=(1 + ) (1 )"j yy;forj=t;b: (35) 36 Journal of Manufacturing Science and Engineering This relation may be imposed after adjustment of the bilayer rest fundamental forms ar;tandar;bon each iteration. Thus, the forward problem will be solved taking into account the plastic anisotropy. Given that the rest fundamental forms are numerically adjusted for each triangular element separately, the efﬁciency of the inverse problem resolution depends on the consistency between the target mesh and the initial mesh. More precisely, it depends on the mapping ~ mbetween the initial 2D shape and the target 3D shape. The general requirement for the mapping is to preserve the shape of each triangular element as well as possible. This minimizes the local eigenstrain assigned by the algorithm and makes the computed free pattern smoother. This problem was not faced during the numerical validation because the target shapes were derived from the initial shapes through the forward problem resolution, so they were optimally meshed by default. The mesh consistency can be ensured by ﬁxing the target shape mesh and by its mapping onto the initial 2D geometry. The ﬁxed initial geometry is an important constraint for the mapping because it involves a ﬁxed 2D boundary. The mapping can be done using the methods oriented on maximal preservation of local angles, such as the Least Squares Conformal Mapping (LSCM) algorithm [32]. Next, local mesh distortions with respect to the target mesh can be minimized using a numerical optimization algorithm. For example, the L-BFGS algorithm that we use for the global elastic energy minimization can cope with this task. CONCLUSION The theory of non-Euclidean plates in combination with the eigenstrain approach provides an extensive theoretical framework for the modeling of shot peen forming. The eigenstrain approach represents the treated plate as a bilayer undergoing nonuniform eigenstrain, and the theory of non-Euclidean plates accurately solves the forward problem for this case. The deformed shape is calculated through minimization of the global elastic energy following analytical gradients. The iterative inverse problem resolution is based on the comparison of geometrical properties of the plate in its current and target conﬁgurations. The adjustment of the prescribed eigenstrain on each iteration is done on a local scale involving simple arithmetic operations and takes negligible amount of time. A low number of iterations (not more than 10) ensures fast resolution of the inverse 37 Journal of Manufacturing Science and Engineering problem. The eigenstrain formulation of the inverse problem makes the algorithm applicable for any type of processes that induce small isotropic eigenstrain. These include, among others, laser peen forming of metal plates or 4D printing of shape-shifting polymer structures. The precision of the inverse problem resolution is independent of plate thickness and its mechanical properties. It is, however, dependent on consistency between the initial and the target meshes. The pattern grouping makes the inverse problem solution practically applicable. It adds uncer- tainty to the solution, but in many cases it decreases the induced error. The grouping algorithm can locally mix up the peening regimes having slightly different intensities. A method for correction of the local grouping errors would enhance the solution quality. The future work implies experimental validation of the proposed inverse problem simulation technique. It will reveal practical constraints that may cause simulation error. Among others, we will examine inﬂuence of the peening parameters and of the material plastic anisotropy on the induced eigenstrain. ACKNOWLEDGEMENTS The authors gratefully acknowledge ﬁnancial support from The Fonds de recherche du Québec – Nature et technologies (FRQNT) and The Aluminium Research Center (REGAL).