We present a methodology for automated forming of metal plat es into freeform shapes using shot peening. The methodology is base d on a sim- ulation software that computes the peening pattern and simu lates the eﬀect of its application. The pattern generation requires pr eliminary experimental characterization of the treatment. The treatm ent is applied by a shot peening robot. The program for the robot is generate d auto- matically according to the peening pattern. We validate the methodology with a series of tests. Namely, we form nine aluminum plates i nto dou- bly curved shapes and we also shape model airplane wing skins. The article describes the complete workﬂow and the experimental r esults. Keywords: Shot peen forming, Optimization, Modeling, Automation, Non-Euclidean plates, Eigenstrains 1 Springer Nature 2021 L ATEX template 2Simulation and automation of aluminium panel shot peen form ing 1 Introduction Shot peen forming (peen forming) is an industrial process fo r forming metal plates into smoothly curved shapes. It consists in projecti ng a stream of rigid particles (shot) towards a metal plate surface at high veloc ity. Numerous over- lapping impacts induce plastic expansion of the outer layer of material, which causes bending of the plate in the treated segment. The set of treated segments over the ﬂat plate surface is called the peening pattern. Alt ering the pattern shape or the treatment eﬀectiveness alters the deformed sha pe as well. This versatility makes peen forming an eﬃcient instrument for ac hieving complex curvature proﬁles on large plates, such as airplane wing ski ns [1] or rocket tank bulkheads [ 2]. In contrast to stretch forming, which is often used for for ming of such components, shot peen forming requires no punch matr ices, dies or intermediate heat treatments [ 3,4]. In addition, shot peen forming allows to shape composite aluminum panels that constitute fuselage s kins [5]. Nowadays, the market proposes a wide range of programmable p een form- ing robots, such as those mentioned in Refs. [ 6] and [7]. Nevertheless, there exists no publicly accessible numerical tool that would aut omatically control a peening robot, bridging the gap between the numerical model (CAD) and the robot path. Such a tool must generate an optimal peening patt ern for a given target shape, i.e., solve the inverse problem , then translate the pattern into a program for the peening robot, and simulate the eﬀect of its a pplication, i.e., solve the forward problem . At least one commercial company claims to have developed an appropriate tool, but the adopted approach is c onﬁdential [ 6,8]. Smaller companies thus use manual peen forming, which is che aper but less eﬃcient. It is a trial-and-error process that provides lowe r repeatability than the automated one. The quality of forming in this case relies completely on the operator’s experience, which involves costly training. Mo reover, it can lead to chronic health problems for the operator [ 8]. In the matter of peen forming simulation, a straightforward modelling of every impact would be computationally unachievable for the case of large industrial parts [ 9]. However, the simulation can be accelerated by means of theeigenstrain representation of the peening-induced loads [ 10]. In the peen forming case, it implies formulating the loads as an anelast ic strain (eigen- strain) introduced in the treated segments analogous to the rmal expansion, so that the individual impacts are not simulated [ 11,12]. The deformation of the plate caused by the eigenstrain is described as the elastic m aterial response, i.e., the springback [ 13]. It can be computed using the shell ﬁnite element method [ 12,14]. This method can also be reformulated in terms of the the- ory of non-Euclidean plates [ 15], which is speciﬁcally designed to describe the reconﬁguration of plates subjected to eigenstrain. The eig enstrain approach together with the theory of non-Euclidean plates provide a f ramework for the inverse problem resolution [ 16]. In terms of the eigenstrain approach, the peening pattern is described by the eigenstrain magnitude and its distribution over the sur face of the plate. However, from a practical point of view, the pattern is descr ibed by the shot Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing3 stream parameters and by the shape of the treated segment. Th e eigenstrain magnitude can be related with the shot stream parameters thr ough direct impact simulation on a small representative volume of mater ial [9,17]. This requires building a speciﬁc simulation model that reﬂects p lastiﬁcation of material induced by the impacts. Otherwise, the eigenstrai n can be quantiﬁed experimentally using small specimens made of the same mater ial as the plates subjected to forming. The specimens must be uniformly treat ed with the same parameters that are applied for treatment of the plates. The peening-induced plastic strain, i.e., eigenstrain, is either measured dire ctly using X-ray diﬀrac- tion [18,19], or computed based on the residual stress measurements [ 12]. These strategies require costly measuring equipment and ca n be replaced by an eigenstrain-based simulation software. Indeed, in such case the eigenstrain is deduced from the deﬂection of a uniformly treated specime n [20,21]. More precisely, the eigenstrain is adjusted in the way that the si mulated deﬂection corresponds to the measured one. Although this method does n ot allow recon- structing the full residual stress proﬁle, it allows quanti fying the macroscopic bending and stretching induced by the treatment, and since i t does not require any speciﬁc equipment or numerical models, it can be eﬃcient ly applied in industry. The eigenstrain induced in an aluminum plate by a perfectly s ymmetrical impact can be anisotropic [ 22]. The main causes for this eﬀect are the material plastic anisotropy [ 23,24] and the initial stresses induced by the fabrication process [ 25,26]. Another potential cause is the prestress, i.e., pre-bend ing of the plate with clamping supports during peening, which is no t present if the plate is held ﬂat during treatment. The inﬂuence of the three causes can be quantiﬁed experimentally using an additional equipment [ 26]. Otherwise, the eigenstrain anisotropy can be characterized using an eigen strain-based sim- ulation software, so that the anisotropy is numerically adj usted to ﬁt the simulations to the experiments. In terms of this method, a pa rticular challenge lies in the fact that the inﬂuence of the eigenstrain anisotr opy on the ﬁnal shape depends on the geometry of the plate [ 26]. Consequently, the character- ization should ideally be done using exactly the same plates as those that are supposed to be shaped afterwards. Precise reproduction of the peening pattern shape is crucia l for the process eﬃciency. One possible approach implies masking of the untr eated segments with subsequent uniform treatment of the plate. The masks ar e glued to the surface of the plate to absorb the energy of the shot and to pro tect the mate- rial underneath. A thick adhesive tape or cardstock cope wel l with this task. The masks are shaped either manually or using a laser cutter a nd are manu- ally glued to the component [ 25,27]. This approach is cumbersome when the applied pattern changes from one component to another, beca use a diﬀerent mask needs to be fabricated for each pattern. For this reason , an automated peen forming system would require a diﬀerent approach to rem ove this manual and time consuming component of the workﬂow. Springer Nature 2021 L ATEX template 4Simulation and automation of aluminium panel shot peen form ing In the existing literature on shot peen forming, the experim ental validation of the numerical forward problem resolution is mostly limit ed to uniform treat- ment of small coupons [ 21,28]. Such experiments do not require an automated peenformingsystemduetothesimplicityoftheappliedpatt ern.Moreover,the existing literature on the numerical inverse problem resol ution in this domain is mostly limited to the description of the developed algori thms [29] or to their numerical validation [ 14,16]. Only Faucheux and Miao et al. have validated their inverse problem resolution strategy experimentally on complex shapes such as a cylinder, a saddle and a wave [ 25,27]. The peening patterns in this research were applied using masking, so the overall peen for ming workﬂow is not completely automated. Moreover, the developed deﬂec tion was of the order of millimeters for 1 ×1 m panels, which is too small to reveal nonlinear geometric eﬀects appearing with large deﬂections. Here, we present a fully automated peen forming workﬂow impl ying a pro- grammable peening robot. We describe the numerical tools, t he steps for their calibrationandthespeciﬁcallydesignedinstallationsfo rtherobot.Thenumer- ical tools include the simulation software and the translat ion software. The former solves both the forward and the inverse problems and a llows to quan- tify the eigenstrain anisotropy. The latter automatically converts the peening pattern into the nozzle paths. Hence, the shot stream “paint s” the peening pattern on the plate, so the necessity for masking is elimina ted. Together, the numerical tools form a “CAD-to-path” solution allowing the user to input a desired shape and to generate the robot program with a minima l intervention. The automated workﬂow was developed for experimental valid ation of the numerical forward and inverse problem solvers described in Ref. [16]. The forward problem solver validation involves randomly gener ated and strongly curved target shapes that have deﬂection of 10 −31 mm. The inverse problem solver validation consists in shaping parts of a model airpl ane wing skin of a predeﬁned geometry. 2 Methodology The methodology section begins with a description of the und erlying princi- ples of our simulation software. Next, we present the transl ation software and the peen forming robot. We then examine our strategy for the e stimation of the experimental error implying 3D scanning. Subsequently , we describe the experimental steps for the software calibration. They incl ude characterization of the shot stream, determination of the eigenstrain induce d by a given treat- ment and quantiﬁcation of the eigenstrain anisotropy on a gi ven component. All the elements mentioned above constitute together the au tomated peen forming workﬂow, which is schematized in Figure 1. The end of the methodol- ogy section is dedicated to the experimental validation of t he forward and the inverse problem solvers. Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing5 g) 3D Scanner h) Error Mapa) Target shapec) Peening patterne) Peening programb) Simulation softwared) Translation software f) Peening robot i) Shot stream j) Eigenstrain k) Anisotropy Fig. 1 Schematic overview of the automated shot peen forming workﬂow. The n umerical model of the target shape (a) is input into the simulation software (b). The latter computes an optimal peening pattern (c), that must be applied to form the plat e into the target shape. The peening pattern is input into the translation software (d), wh ich generates the peening program (e). The peening program includes the nozzle paths and a set of p eening parameters that allow the peening robot (f) to reproduce the peening pattern. Th e peening program is loaded into the peening robot in the form of a programming code. When the treatment is ﬁnished, we characterize the curved shape of the plate using a 3D scanner (g ) and trace an error map (h), which shows discrepancy between the target and the scanned sha pes. The simulation and the translation softwares require experimental calib ration, which consists in characterizations of the shot stream (i), of the eigenstrain (j) ind uced by the peening treatment, and of the eigenstrain anisotropy (k). The shot stream is ch aracterized with ﬂat dummy specimens, the eigenstrain is measured using 76 ×19 mm aluminum coupons, and the anisotropy is quantiﬁed with the help of square aluminum plates. Ar t by IMPAKT Scientiﬁk. 2.1 The simulation software 2.1.1 Representation of peening-induced loads with eigenstrain Following Ref. [ 16], we virtually divide the plate into two layers of equal thic k- ness to simulate the curving eﬀect of peen forming. We assign in-plane swelling to one layer and in-plane shrinking to the other in the spots w here the plate is treated. The swelling layer is located on the side that und ergoes treatment. Springer Nature 2021 L ATEX template 6Simulation and automation of aluminium panel shot peen form ing Theswellingandshrinkingareassignedintheformofeigens trainandareanal- ogous to thermal expansions [ 12]. The eigenstrain magnitude is constant in the through-thickness direction for each layer but varies alon g the plate surface according to the peening pattern. We denote the eigenstrain tensors imposed on the top and bottom layers as εtandεbrespectively. The incompatibility between εtandεbforces the plate to bend. We endow the component with Lagrangian coordinates x,y,z, where the x−andy−axes follow the plate mid-surface. Both eigenstrain tensor s in these coordinates have the following form [ 25]: εj= εj xx0 0 0εj yy 0 0 0−(εj xx+εj yy) forj=t,b. (1) For a plate with uniform properties and initial stresses, th e proportion betweenεj xxandεj yyis constant along the whole component. It is expressed with the anisotropy coeﬃcient χ[26]: /braceleftBigg εj xx= (1+χ)εj ∗, εj yy= (1−χ)εj ∗,forj=t,b, (2) whereεj ∗= (εj xx+εj yy)/2. If the material is plastically isotropic and free of stresses then εj xx=εj yy=εj ∗. 2.1.2 Geometrical shape description Provided that peen forming deals with thin plates, we associ ate the plate shape with the shape of its mid-surface. We numerically disc retize the mid- surface with triangular elements. Following Ref. [ 30], we describe the shape of each element with two matrices: aandb. The former represents the local ﬁrst fundamental form and the latter represents the local second fundamental form, that respectively characterize the local stretching and th e local curvature [ 31]. The set of fundamental forms computed for each element uniqu ely describes theshapeofatriangulatedsurface.Foratriangularelemen tdeﬁnedbyvertices /vector v0,/vector v1,/vector v2and edges/vector e0=/vector v1−/vector v0,/vector e1=/vector v2−/vector v1,/vector e2=/vector v0−/vector v2the two matrices are computes as: a=/bracketleftbigg /vector e1·/vector e1/vector e1·/vector e2 /vector e2·/vector e1/vector e2·/vector e2/bracketrightbigg ; (3) b=/bracketleftbigg /vector e1·2(/vector n0−/vector n2)−/vector e1·/vector n0 −/vector e1·/vector n0/vector e2·2(/vector n1−/vector n0)/bracketrightbigg . (4) Here,/vector nistands for a vector normal to the edge /vector ei,i= 1,2,3. The angle of inclination of /vector niwith respect to the plane enclosing the triangular element i s Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing7 variable and thus allows a more ﬂexible shape description. F igure2illustrates the components of a triangular element that participate in c omputation of the fundamental forms. The local principal curvatures can be computed for each elem ent as the eigenvalues of the shape operator S=a−1b[32,33]. The local mean curvature His the average of the two principal curvatures, and the curva ture along any direction may be deduced from the principal curvatures usin g Euler’s formula [32]. 2.1.3 The forward problem resolution In terms of peen forming, the forward problem consists in the determination of the deformed shape of the component provided with the eige nstrain pat- tern. For this, we rely on the approach presented in Ref. [ 30]. We assume that the plate is free to deform and that the treatment is appl ied instantly from both sides. We also assume that the initial conﬁguratio n is ﬂat, and we denote the initial lengths of the element edges as /vector e0,/vector e1,/vector e2. We consider the eigenstrain constant inside a triangular element. An aniso tropic eigenstrain (εxx,εyy) prescribes new restedges/vector er 0,/vector er 1,/vector er 2: /vector er i= (εxx+1)/vector eix+(εyy+1)/vector eiyfori= 1,2,3. (5) Here,/vector eixis the projection of /vector eion thex-axis and/vector eiyis the projection of /vector eion they−axis. New edges induce a new restﬁrst fundamental form ar: ar=/bracketleftbigg /vector er 1·/vector er 1/vector er 1·/vector er 2 /vector er 2·/vector er 1/vector er 2·/vector er 2/bracketrightbigg . (6) A diﬀerent eigenstrain ( εt xx,εt yy) and (εb xx,εb yy) imposed on the top and bottom layers respectively induce diﬀerent ﬁrst fundamental form sar,tandar,b. Con- currently, all components of the rest second fundamental fo rmsbr,tandbr,b equal zero because the edge normals stay perpendicular to th e plane englobing the ﬂat initial conﬁguration. Fig. 2 A triangular mesh element and the vectors serving for computation of the local fundamental forms. The vectors /vector v0,/vector v1,/vector v2deﬁne positions of the element vertices in the global coordinate system. The edge vectors /vector e0,/vector e1,/vector e2are computed as the diﬀerence between the corresponding vertex vectors. The vectors /vector n0,/vector n1,/vector n2are traced in the center of the correspondingedgesanddeﬁnethelocalcurvatureofthemesh.Theyareno rmaltotheedges, but their angle of inclination with respect to the element constitutes a deg ree of freedom. Springer Nature 2021 L ATEX template 8Simulation and automation of aluminium panel shot peen form ing Prescription of a diﬀerent eigenstrain to the top and bottom layers does not make the plate delaminate. Instead, the plate adopts a cu rved conﬁgura- tion described with the ﬁnalfundamental forms afandbf. Hence, from the mathematical point of view, the forward problem resolution consists in the determination of afandbftakingar,tandar,bas the input. This is done through numerical minimization of the global elastic energ y functional EBL [16]: EBL=1 2K/summationdisplay k=1/bracketleftbiggh 8/vextenddouble/vextenddouble(a−1 r,b)k(af)k−I/vextenddouble/vextenddouble2 e+h3 24/vextenddouble/vextenddouble(a−1 r,b)k(bf)k/vextenddouble/vextenddouble2 e +h2 8/angbracketleftbig/parenleftbig (a−1 r,b)k(af)k−I/parenrightbig ,(a−1 r,b)k(bf)k/angbracketrightbig e/bracketrightbigg/radicalBig det(ar,b)k +1 2K/summationdisplay k=1/bracketleftbiggh 8/vextenddouble/vextenddouble(a−1 r,t)k(af)k−I/vextenddouble/vextenddouble2 e+h3 24/vextenddouble/vextenddouble(a−1 r,t)k(bf)k/vextenddouble/vextenddouble2 e −h2 8/angbracketleftbig/parenleftbig (a−1 r,t)k(af)k−I/parenrightbig ,(a−1 r,t)k(bf)k/angbracketrightbig e/bracketrightbigg/radicalBig det(ar,t)k.(7) Here, the index kis attributed to the fundamental forms belonging to the elementk, andKdenotes the total number of elements in the model. The matrixIis the 2×2 identity matrix. The elastic norm/vextenddouble/vextenddouble·/vextenddouble/vextenddoubleand the elastic energy inner product/angbracketleftbig ·,·/angbracketrightbig denote the following operations: /vextenddouble/vextenddoubleA/vextenddouble/vextenddouble=αTr2(A)+2βTr(A2), (8) /angbracketleftbig A,B/angbracketrightbig =αTr(A)Tr(B)+2βTr(AB), (9) whereAandBare arbitrary matrices, α=Yν/(1−ν2) andβ=Y/(2+2ν). The values νandYstand respectively for the Poisson’s ratio and the Young modulus of the treated material. We minimize the elastic energy functional with the help of th e Limited- memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optim ization algo- rithm [34] using the analytically computed gradients [ 35]. 2.1.4 The inverse problem resolution In terms of peen forming, the inverse problem consists in the determination of the eigenstrain tensors εtandεbimposed on each element taking the target shape as input. Another necessary input is the initial shape geometry. Our algorithm for the inverse problem resolution relies on full consistency between the initial and the target meshes. For this reason, the initi al mesh is generated as a function of the target mesh, which is illustrated in Figu re3. It is done in two steps. Firstly, the target mesh is projected on the initi al geometry using the Least Squares Conformal Mapping (LSCM) algorithm [ 36]. Secondly, the initial vertex positions are numerically optimized using t he L-BFGS algorithm Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing9 Conformal mapping + Minimization of distortion b) Initial geometrya) Target mesh c) Initial mesh Fig. 3 Generation of the initial mesh. The mesh of the curved target shape (a) is projected on the ﬂat initial geometry (b) using the LSCM algorithm. The dist ortion of each element in the projected mesh is then minimized with respect to the target mesh using the L-BFGS algorithm. The resulting mesh (c) is used as the initial mesh during the inverse problem resolution. to minimize the element distortion with respect to the targe t mesh. The initial shape boundary is kept ﬁxed during the optimization. The target shape is described with two fundamental forms atarandbtar. Our algorithm for the inverse problem resolution is iterati ve. Each iteration implies an adjustment of the eigenstrain pattern and subseq uent computation ofthenewcurrentshapeprovidedbytheadjustedpattern.Th eeigenstrainten- sors are adjusted for each triangular element separately ba sed on comparison of the current shape with the target shape. More precisely, o n each iteration, we multiply εt ∗andεb ∗by the ratio between the local mean curvatures of the target shape Htarand of the current shape Hc: εj,new ∗=εj ∗Htar Hcforj=t,b. (10) Next, we solve the forward problem for ( εt,new ∗,εb,new ∗) and thus obtain a new current shape. These two steps are repeated until stabiliza tion of the cur- rent shape. The anisotropy coeﬃcient χis pre-deﬁned with a calibration step described in Section 2.5.3and kept constant throughout the iterations. The initial guess is made in terms of ar,tandar,bas presented in Ref. [ 16]. Thus, we ﬁrst compute the ar,tandar,bas: /braceleftBigg ar,t=atar−2h 3btar, ar,b=atar+2h 3btar.(11) Then, for each layer, we take the average of the eigenstrains prescribed by ar,t andar,bin each direction and use this eigenstrain as the initial gue ss. Springer Nature 2021 L ATEX template 10Simulation and automation of aluminium panel shot peen form ing This algorithm constitutes a simpliﬁed version of that desc ribed in Ref. [16]. Namely, the eigenstrain is adjusted in the way that it lead s to the target curvature , butthe stretching of thetarget shapeisnot taken intoaccount. This formulation was adopted because the target shapes implied l arge rotations but small stretch with respect to the initial state. The ﬁnal step of the inverse problem resolution is the patter n grouping [16], which is done after the last iteration. The aim of this proc ess is to split the peening pattern into segments treated with practically attainable peen- ing regimes. Each attainable peening regime is related to th e corresponding inducedeigenstrainduringthecalibrationphase.Thegrou pingalgorithmtakes the eigenstrain assigned to each element and reassigns it th e closest eigenstrain among the limited number of attainable ones. This removes gr adual variation of the eigenstrain along the plate surface and divides the pa ttern into strictly deﬁned peened segments. When the grouped pattern is compute d, we solve the forward problem taking this pattern as input. This step serv es for preliminary estimation of the error between the target shape and the ﬁnal shape developed with the ready-to-use pattern. 2.1.5 Inverse identiﬁcation of the anisotropy coeﬃcient This problem is similar to the inverse problem, except for th e condition that the local eigenstrains εt ∗andεb ∗are pre-deﬁned and the anisotropy coeﬃcient χmust be determined. We assume that χis a global property of the treated material, so it is uniform along the plate surface. The initi al mesh is generated based on the target mesh in the same way as it is done during the inverse problem resolution. The problem is solved through the iterative adjustment of χbased on the comparison of the current and target shapes. We compute loca l curvatures along thex-direction for both shapes: κc xfor the current shape and κtar xfor the target shape. We average their absolute values κc avgandκtar avgand then adjust χas follows: χnew= (1+χ)κtar avg κcavg−1. (12) We then adjust the local eigenstrain tensors according to Eq n. (2) and solve the forward problem. Starting with an initial guess of χ= 0, the pro- cess is repeated until the convergence of κc avgandκtar avg, providing us with an evaluation of the anisotropy coeﬃcient χ. 2.2 The translation software The translation software ﬁlls the gap between the numerical peening patterns and the actual peen forming process. The peening pattern tha t it takes as input is provided in the form of a mesh with the eigenstrain as signed to each of the elements. The pattern is divided into segments with di ﬀerent prescribed Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing11 eigenstrain that correspond to treatment with diﬀerent pee ning regimes. The relation between each peening regime, i.e., each set of the p eening parameters, and the corresponding eigenstrain is deﬁned during the cali bration phase and savedinthesoftware.Thesoftwareﬁrstdetectsbordersoft hesegmentstreated with the same peening regime. It then generates the nozzle pa ths that ﬁll each segment and assigns the peening parameters corresponding t o each segment. Finally, the software wraps up the paths and the peening para meters in a programming script written in Karel programming language [ 37], which is then loaded in the peening robot. The script controls all paramet ers of the robot operation and does not require any manual adjustment. The translation software is able to fulﬁll a peened segment w ith two types of nozzle paths, which are presented in Figure 4. The ﬁrst path type is called “Zigzag” and implies straight parallel nozzle paths. The se cond type is called “Circular”. The paths generated in this way constitute a set of concentric closed lines, and the longest line is the outline of the treat ed segment. The “Circular” paths better ﬁll complexly shaped segments. How ever, the neigh- boring paths of this type may not always have a constant oﬀset , that results in a less uniform coverage. The software allows adjusting several common parameters fo r both path types. Thus, the outline of the treated segment can be smooth ed with respect to its actual border prescribed by the peening pattern. This allows ﬁltering out featurestoosmalltopeenwithrespecttotheshotstreamwid th.Moreover,the outlinecanbeoﬀsetwithrespecttotheborder,whichallows scalingthepeened segment without resorting to the simulation software. The o ﬀset between the paths can also be adjusted depending on the shot stream width , which is measured during the calibration phase. The path oﬀset is adj usted in the way that the peening coverage in the treated segment is as unifor m as possible. Figure5illustrates the adjusted outline and the concept of the path oﬀset. The speed of movement of the nozzle along the peening paths ca n also be outline of the peened segment on the top side nozzle paths"Circular" paths "Zigzag" paths outline of the peened segment on the bottom side Fig. 4 Two types of nozzle paths generated by the translation software. The exa mple pattern implies treatment from the top and bottom sides with the same peening regime. Springer Nature 2021 L ATEX template 12Simulation and automation of aluminium panel shot peen form ing Path offset Border according to the peening pattern Smoothed border Outline Nozzle paths Fig. 5 Adjustment of the peened segment executed with the translation software. The outline of the peened segment is smoothed with respect to its initial border and oﬀset to increase the size of the treated segment. The path oﬀset is adjusted to ensure un iform coverage. adjusted depending on the peening regime and is usually set t o 50-100 mm/s. The peening robot is not able to block the shot stream during t he program execution, so the nozzle moves between diﬀerent peened segm ents and diﬀerent paths inside one segment at a high speed of 1 500 mm/s. 2.3 The peening robot The peening robot involved in the experiments is based on a Ca nablast shot peening machine. It possesses a FANUC M20iA programmable ro botic arm, which guides the peening nozzle. The robotic arm and the trea ted component are installed inside a closed peening cell shown in Figure 6. We developed a speciﬁc peen forming ﬁxture to treat square pa nels of 320 × 320 mm. The ﬁxture is designed in two main parts: the ﬁxed and t he removable frames. The ﬁrst ensures a precise location of the peened par t in terms of the robot’s reference frame and the latter clamps the plate alon g its perimeter. This system eases the installation process and allows remov ing the plate from the shot peening cell while keeping it under the clamping con straints. The ﬁxture is mounted on a rotary axis, so the plate can be trea ted from both sides with no need for manual reinstallation. The rotar y axis turns the plate by 180◦in 2 seconds, and during this process the programming script points the nozzle aside, that prevents undesired impacts. 2.4 3D scanning and error map Onceaplateistreated,wecompareitscurvedshapewithnume ricalprediction. WeremovetheplatefromtheframeandwescanitusingaHexago nStereoScan Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing13 Media         NozzleFANUC M20iA Removable frame Fixed frame Rotary axis Fig. 6The interior on the robotic peen forming cell. The programmable ro btic arm FANUC M20iA guides the peening nozzle, and the media is supplied to the nozzle th rough the elastic hose. The treated plate is clamped by the removable frame along its perimeter. The removable frame is mounted on the ﬁxed frame attached to the rotar y axis. neo 3D scanner. Before scanning, the plate is covered with Ma gnaﬂux SKD- S2 aerosol to make its surface non-reﬂective. The result of 3 D scanning is a meshed surface describing the plate shape. Next,wenumericallysuperimposeeachscannedshapewithth ecorrespond- ing predicted shape and we compute the Hausdorﬀ distance dH[38] between the two shapes. We nondimensionalize the Hausdorﬀ distance by dividing it by the deﬂection of the target shape ∆ z, which is computed as the smallest dimension of a box bounding the plate. Thus, we obtain a nondi mensionalized error Ω z: Ωz=dH ∆z. (13) Figure7illustrates the workﬂow for determining Ω z. The optimal alignment of the two shapes is done through least squares estimation of th e transformation parameters [ 39]. Springer Nature 2021 L ATEX template 14Simulation and automation of aluminium panel shot peen form ing 2.5 Calibration steps 2.5.1 The shot stream characterization We ﬁrst estimate the shot stream intensity for each peening r egime that we use. It is done following the standard procedure involving t ype A Almen strips and the Almen gage [ 40]. Each peening regime used for the treatment of the plates is applied until saturation, and the intensity measu rement is necessary to estimate the appropriate number of nozzle passages. In ot her words, we estimate how many times the nozzle must pass along the paths p rescribed by the translation software to achieve saturation in the whole treated segment. The nozzle paths are lines that direct the center of the nozzl e, but the actual shot stream has a ﬁnite width. Hence, we measure the di stribution of shot inside the stream to deﬁne the nozzle path oﬀset. It is ma inly inﬂuenced by the distance from the nozzle to the treated component. Wep erform this test using a 418 ×127 mm ﬂat dummy specimens. We make one peening passage in the center of the specimens along the longitudinal direct ion translating the nozzle at a high speed of 1 500 mm/s and then scan the specimen w ith a doc- ument scanner. Next, we post process the scanned image, i.e. , we highlight the white pixels corresponding to impacts and we darken the back ground corre- sponding to the untreated surface. The picture resolution i s 4830×1470 pixels, and we count the number of white pixels in each pixel column ou t of 1470. Figure8shows an example of the shot distribution measured in this wa y. We set the paths oﬀset to the largest value that ensures uniform coverage of the treated segment with a maximal variation in coverage of 1%. T his value equals approximately half of the shot stream width but is speciﬁc fo r each particular distribution. z 0(mm) a  b c  d  Fig. 7 Determination of the dimensionless error Ω z=dH/∆zbetween the simulated shape and the shape of the treated plate. a) The simulated shape and its boun ding box. The deﬂection ∆ zis computed as the smallest dimension of the bounding box. b) A photo o f the plate treated following the computed pattern that leads to the si mulated shape. The non-reﬂective surface is the treated part. c) The meshed shape scanned with a 3 D scanner. d) Map of the local distances dbetween the simulated and the scanned shapes that were optimally aligned. The map is traced on the simulated shape. The Hausd orﬀ distance dHis the maximal value of local distances d. Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing15 Fig. 8Spatial distribution of shot on a 418 ×127 mm ﬂat dummy specimen that underwent one passage of the peening nozzle in the longitudinal direction. The n ozzle travelled at a speed of 1500 mm/s, and the nozzle distance in this case equaled 152.4 mm. Th e graph indicates the percentage pof the peening-aﬀected surface as a function of the transverse coordinate d. More precisely, we virtually divided the specimen into 1470 columns tr aced in the longitudinal direction, and the percentage pdenotes the ratio of the peening-aﬀected area in each column with respect to its total area. 2.5.2 Determination of the eigenstrain Every combination of peening regime and treated material in duces diﬀerent εt ∗andεb ∗. We determine εt ∗andεb ∗experimentally using rectangular coupons of size 76 ×19 mm. The coupons are cut with waterjet cutting from the same aluminumplatesasthosethatareshapedafterwards.Thecou ponsareclamped in the Almen holder and treated until saturation. The bending behavior of the coupon is governed by the dimensi onless load ΓB[26]: ΓB=3 2(1+ν)∆ε, (14) where ∆ε=εt ∗−εb ∗. For small values of Γ B, the curvature κmeasured at the center of the treated coupon is the same in any direction. Mor eover, in this case the load Γ Bis proportional to the dimensionless curvature κh[26]. The treated coupon has a convex shape ( κis negative), and the relation takes the form: κh=−µΓB. (15) The proportionality constant µcan be determined using numerical simulations (see Appendix A). Therefore, ∆ εcan be deduced from κ. According to the simulations, the linear relation holds at least up to Γ B= 8×103. The curvature κ, in turn, is deduced from the coupon deﬂection u, which is measured using the Almen gage [ 41]. We can approximate the shape of a curved coupon to that of an elliptic paraboloid [ 42], which allows to express u as: Springer Nature 2021 L ATEX template 16Simulation and automation of aluminium panel shot peen form ing u=−κl2 x+l2 y 2. (16) Here,lx=15.87mmand ly=7.94mmarethedimensionsoftherectangleformed by the supports of the Almen gage [ 41]. The dimensionless curvature κhdepends on ∆ εand does not depend on the sum (εt ∗+εb ∗), which only inﬂuences stretching of the plates [ 16]. However, shot peening induces small eigenstrain [ 9], so we conclude that the peening- induced stretching is also small and assume that εt ∗=−εb ∗. Combining Eqns. (14–16) yields:   εt ∗=2hu 3µ(1+ν)(l2x+l2y); εb ∗=−2hu 3µ(1+ν)(l2x+l2y).(17) For better stability of the results, we treat four coupons wi th each regime and compute their average deﬂection uavg. We then substitute uavgforuin Eqn. (17) and thus ﬁnd εt ∗andεb ∗. 2.5.3 Determination of the anisotropy coeﬃcient Aluminum plates are produced by rolling, and this process is one of the key factorsinducingtheeigenstrainanisotropy.Namely,thep eeninginducedeigen- strain is diﬀerent in the rolling direction (L) and in the in- plane transverse direction (T). This eﬀect alters the local curvatures induc ed by peening in the L- and T-directions. Thus, if peen forming is simulated a ssuming that the material is isotropic, then in practice the local curvature in one direction is bigger than that simulated, and the local curvature in the ot her direction is smaller [ 26]. A similar eﬀect is also observed during peen forming of met al laminates containing oriented ﬁlaments [ 43]. We quantify the eigenstrain anisotropy using the “wave” exp eriment. For this, we use the plates that have the same size and are produce d in the same way as those used for the forming experiments. We treat two pl ates following the “wave” pattern illustrated in Figure 9. This pattern consists in treating half of the plate from one side and the other half from the othe r side. The orientation of the “wave” pattern controls the bending dire ction of the plate. Hence, for the ﬁrst plate, the halves are separated along the T-direction thus forcing the plate to bend in the L-direction. On the contrary , the same pattern rotated by 90◦makes the second plate bend along the T-direction. We scan both treated plates using the 3D scanner and obtain th eir meshed numerical models. These models serve as target shapes for nu merical determi- nationoftheeigenstrainanisotropycoeﬃcientasdescribe dinSection 2.1.5.We orient the coordinate system in the way that the x-axis points in the “wave” direction, i.e., in the L-direction for the ﬁrst plate and in the T-direction for the second plate. To compute the κtar avg, we take into account only the elements Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing17       200 mm200 mm Fig. 9The “wave” experiment. The “wave” peening pattern is applied on two 320 ×320 mm sized aluminum plates. The pattern forces the plates #1 and #2 to bend alo ng the L- direction (a) and the T-direction (b) respectively. The plates are scanned in 3D after peening. Figures c) and d) show the deformed state of the plates #1 and #2 respectivel y. For each plate, we compute the average absolute value of κxfor all the elements contained in the 200×200 mm square traced in the center of the plate (blue), and thus obtain κtar avg. that fall into a 200 ×200 mm square traced in the center of the plate (see Figure9). The curvature of the other elements is not taken into accou nt to reduce the inﬂuence of the edge eﬀects. In this way we determine two anisotropy coeﬃcients: χLfor the ﬁrst plate andχTfor the second plate. One of them is positive and the other is n egative. We then compute the average of their absolute values χ∗= (|χL|+|χT|)/2 and thus ﬁnd the eigenstrain anisotropy coeﬃcient for the gi ven material. For the subsequent simulations we make a convention that the x-axis follows the principal eigenstrain direction, so that χ∗is positive. Theoretically, χLmust equal −χT, so the “wave” treatment of one plate should be suﬃcient to determine χ∗. However, exact equality is rarely observed because of practical constraints. For example, the simulat ions do not take into account the initial stresses, which are diﬀerent for each pa rticular component, or the inﬂuence of clamping. For this reason, treatment of tw o plates stabilizes the measured result and allows to check if εt ∗andεb ∗were correctly determined in advance. Moreover, the usage of more than two plates could further improve the precision of this experiment. Springer Nature 2021 L ATEX template 18Simulation and automation of aluminium panel shot peen form ing 2.6 Experimental validation of the simulation software The experiments consisted in shaping 320 ×320 mm aluminum plates. For this, we ﬁxed six peening regimes described in Table 1. The peening equipment was charged with the cut wire CW28 media having shot diameter of 0 .71 mm. 2.6.1 Validation of the forward problem solver We shaped nine plates following randomly generated peening patterns and compared the developed shapes with the numerical predictio ns. The plates were2.06mmthickandweremadeofaluminum6061-T6.Thepatt ernsimplied one peening regime listed in Table 1as #1. To ensure that each pattern induces a diﬀerent shape, we gene rated the patterns based on random target shapes. At ﬁrst, we generate d nine diﬀerent shapes using the peaksfunction embedded in Matlab programming language. We then solved the inverse problem for each shape and thus obt ained the peening patterns. We took into account the measured coeﬃcie nt of eigenstrain anisotropy, i.e., we imposed the ratio (1 −χ∗)/(1+χ∗) between the induced eigenstrain in the L- and T-directions. The target shapes were generated without any constraints on their deﬂec- tion. The inverse problem solution showed that their develo pment implied much larger eigenstrain magnitude than those achievable wi th regime #1. Thus, according to the forward problem resolution, the simu lated shapes induced by the grouped pattern were less curved than the init ially generated target shapes. Consequently, we used these less curved simu lated shapes as reference for the experiments. 2.6.2 Validation of the inverse problem solver The inverse problem resolution was validated using 320 ×320×1.6 mm plates madeofaluminum 2024-T3. Theobjectiveofthisvalidation w astoshapeparts of the upper and lower skins of a model airplane wing. The wing had constant cross-section in the form of NACA 5430 airfoil [ 44] with the total chord length RegimeMass ﬂow rate (kg/min)Nozzle distance (mm)Pressure (kPa)Intensity (A)Saturation (Number of passes) #1 2.27152.4 248.2 23.5 5 #2 152.4 172.4 19.2 4 #3 152.4 137.9 17 4 #4 152.4 103.4 13.6 4 #5 152.4 68.9 10.9 8 #6 304.8 48.2 8 9 Table 1 Parameters of the peening regimes used in the experimental campaign. The n ozzle was oriented perpendicularly to the treated surface. The mass ﬂow rate det ermines the mass of media that passes through the nozzle. The nozzle distance is the distance from the nozzle to the treated surface. The pressure is created by compressed air. The peening int ensity and the number of passes to saturation were measured following the stand ard procedure using type A Almen strips [ 40]. The speed of the nozzle travelling was ﬁxed at 100 mm/s. Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing19 of 500 mm (see Figure 10). The two target shapes had the same area as the ﬂat plates: 320 ×320 mm. We ﬁrst solved the inverse problem for both target shapes and obtained the peening patterns. We then solved the forward problem to pred ict the dimen- sionless error Ω zwith respect to the target shapes. Next, we translated the peening patterns into programs for the peening robot and tre ated the alu- minum plates. Finally, we compared the experimentally deve loped shapes with the target shapes constituting the NACA airfoil in terms of Ω z. The peening patterns consisted of segments treated with the regimes #2- 5 from Table 1. The plates were also uniformly treated from both sides with regime #6 before the application of the peening pattern. Suc h low intensity pre-treatment is commonly used in industry to enhance the fa tigue life of the components. It also makes the surface ﬁnishing more uniform . 3 Results 3.1 Calibration steps The width of the peening-aﬀected band was approximately 4 cm for the nozzle distance of 152.4 mm and approximately 8 cm for the noz zle dis- tance of 304.8 mm. The largest path oﬀsets that provided unif orm coverage corresponded to 19 mm and 38 mm respectively. We determined the eigenstrain using 76 ×19 mm aluminum coupons (see Section2.5.2). The 2.06 mm thick coupons made of aluminum 6061-T6 were treated to saturation from one side with regime #1. The 1.6 mm thick coupons made of aluminum 2024-T3 were ﬁrst treated to saturation fro m both sides with regime #6 and then treated to saturation from one side wi th the regimes #2-5. The average deﬂection of the four coupons treated with each regime is presented in Table 2. Table2also shows the eigenstrain εt ∗induced by each regime. It was computed as a function of deﬂection uavgusing the Eqn. ( 17). We observed a slight diﬀerence in deﬂection of less than 1 % of average in each bunch of the four coupons. It is explained by the presence of i nitial stresses coming from the manufacturing process. To characterize the eigenstrain anisotropy (see Section 2.5.3), we treated the 6061-T6 plates following the “wave” pattern with regime #1. The 2024- T3 plates were ﬁrst treated uniformly with regime #6 and then underwent the “wave” treatment with regime #2. The anisotropy coeﬃcients revealed by the Fig. 10 The two target shapes for the 320 ×320 mm sized panels made of aluminum 2024- T3. The target shapes constituted parts of the upper (blue) and lower (o range) airplane wing skins. The wing was shaped as the NACA 5430 airfoil with the chord l ength of 500 mm. Springer Nature 2021 L ATEX template 20Simulation and automation of aluminium panel shot peen form ing “wave” experiment for the two materials are presented in Tab le3. The 6061- T6 plates demonstrated a 4.6 times stronger eigenstrain ani sotropy than the 2024-T3 plates. 3.2 The forward problem solver validation The nine random peening patterns used for validation of the f orward prob- lem solver are presented in Figure 11. The patterns and the corresponding simulated shapes were numbered from #1 to #9. Figure 12(a-f) shows three representative cases demonstrating the smallest, the bigg est and an average dimensionless error Ω z. This ﬁgure also traces the local dimensionless mean curvatureHhof each shape, showing that the curvature was qualitatively well simulated. The simulations were performed on the same i rregular mesh as that generated by the 3D scanner to allow element-wise cur vature compar- ison. Thus, the irregular mesh was projected onto the initia l geometry using the LSCM algorithm and optimized using the L-BFGS algorithm in the same way as during the inverse problem resolution (Figure 3). The dimensionless error lied in a compact range between 20.5 % and 28.5% for seven shapes out of nine. This means that the Hausdorﬀ dis tancedHgen- erally increases with ∆ z, which is illustrated in Figure 12(g). The shapes #4 and #7 demonstrated signiﬁcantly smaller Ω zthan the others (8% and 12% respectively), which is explained by the following reason. The developed jig did not completely prevent deformation of plates during treatment. It clamped the perimeter but let the interior par t deform. At the same time, the peening pattern was ﬁrst fully reproduced fro m the top side RegimePre- treatment (regime)Material h(mm)uavg (mm)σu (mm)εt ∗×103 #1 - 6061-T6 2.05 0.341 0.003 1.19 #2 #6 2024-T3 1.60.543 0.005 1.47 #3 0.302 0.015 1.2 #4 0.443 0.004 0.82 #5 0.174 0.007 0.47 Table 2 Results of the eigenstrain measurements performed using the bunches of fou r 76×19 mm sized aluminum coupons (see Section 2.5.2). The coupons were ﬁxed during peening with an Almen holder, and their deﬂection was measured with an Alm en gage. The thickness of the coupons is denoted by h. The average deﬂection in each bunch of four coupons is denoted by uavg, and the parameter σudescribes the standard deviation in deﬂection. The eigenstrain induced in the top layer is denoted by εt ∗. In this study we suppose that the eigenstrain εt ∗equals−εb ∗, which stands for the eigenstrain induced in the bottom layer. Material χLχTχ∗ 6061-T6 -0.143 0.193 0.168 2024-T3 0.013 -0.06 0.0365 Table 3 The anisotropy coeﬃcients determined using the “wave” experiment for the t wo materials. The coeﬃcient χ∗is the average of the absolute values of χLandχT. The coeﬃcientχ∗was used as the eigenstrain anisotropy coeﬃcient during the simulations. Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing21 # # # # # # # #  # -     de trea  t   -      e t   -     m side tr  e t    Fig. 11 The nine peening patterns that were applied during the experimental f orward prob- lem validation. Each pattern implied treatment from both sides. The p lates were clamped along their perimeter, so the clamped part did not undergo any treatment . For this reason, each pattern has margins along the plate edges. The patterns were app lied using regime #1 without any pretreatment. and then from the bottom side. Consequently, when the treatm ent was applied from the bottom side, the plates were already bent towards th e top side. This created an eﬀect of prestress that decreased the eigenstrai n in the segments treated from the bottom side [ 7]. Accordingly, the curvature in these segments was smaller than the simulated one. The induced prestress varies locally along the plate and is a lso dependent on the peening pattern. We quantiﬁed this eﬀect with the para meterψthat indicated the ratio of areas treated on the top and bottom sid es: ψ=Atop−Abot Atotal, (18) whereAtopandAbotare the areas of the segments treated on the top and bot- tom sides respectively, and Atotalis the total plate area. Figure 12(h) plots Ωzas a function of ψ. This graph shows that the plates #4 and #7 had the smallestψamong the others, which means that these plates were prestre ssed less than the others when peened from the bottom side. Accord ingly, these Springer Nature 2021 L ATEX template 22Simulation and automation of aluminium panel shot peen form ing51  1 5 20 2 5 3  3 52 3 456781234 56 789 (mm) ∆ (mm) (%)  5 0  5 5 1 1 520 2 5 3 12 34 56 78 9  Best case (#4)Simulated shapes Average case (#8) Worst case (#3) Scanned shapes ! "$ "g %h %=8% =25% =28.5%Best case (#4) Average case (#8) Worst case (#3) d) e) f) 0.0035Hh 0 -0.004 0.0035Hh 0 -0.004 Ω Ω Ω Ω Fig. 12 Results of the forward problem validation. We present three representat ive simu- lated shapes, which showed the smallest (a), an average (b) and the biggest (c) discrepancy with the experimentally developed shapes. The 3D scanned shapes developed in pr actice are marked as d), e) and f) respectively. The colormap indicates the local dim ensionless mean curvatureHh. The discrepancy was quantiﬁed in terms of the dimensionless error Ω z(Eqn. 13). The graph g) traces the Hausdorﬀ distance dHfor all nine cases depending on the plate deﬂection. The graph h) traces Ω zas a function of the parameter ψ, that indicates the dif- ference between the treated areas on the top and bottom sides (Eqn. 18). Each point on the graphs is marked with the number of plate that it corresponds to. shapes showed the smallest Ω z. Essentially, the problem of an undesirable pre- stress may be mitigated in practice by constant alternation of the peened sides during treatment. Moreover, this problem becomes less impo rtant for stiﬀened industrial parts, because the stiﬀeners block bending of th e component in the longitudinal direction unless they are not treated themsel ves [45]. Another issue that increased Ω zin all the cases consisted in the shot stream width. Thus,theshotstreamwasnotsuﬃcientlynarrowtopre ciselyreproduce the borders of the treated segments, which is illustrated in Figure13. This eﬀect had stronger negative inﬂuence on patterns prescribi ng highly irregular Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing23 borders to the peened segments. Hence, the plates #5 and #9 de monstrated a highΩ zwhilehavingarelativelysmall ψ(seeFigure 12h),becausetheborders of the corresponding peened segments were strongly curved ( see Figure 11). On the contrary, the plates #1 and #6 demonstrated a moderate Ωzwhile having a large ψ, because the corresponding patterns prescribed smoothly shaped peening segments. The described eﬀect may be mitigat ed technically with a nozzle that provides a narrower shot stream. The third source of error consisted in the forces and the mome nts induced by clamping, which inﬂuenced local curvatures inside the tr eated areas. Thus, the simulations prescribed uniform curvature inside each t reated segment, but the experiments revealed that the curvature along the perim eter of the treated segmentwashigherthaninthemiddleofthesegment.Thiseﬀe ctwasobserved onallplatesandisvisibleinFigure 12(a-f).Theinﬂuenceofthiseﬀectdepends on the shape of the peening pattern, and its quantiﬁcation is left for further studies. Nevertheless, we expect the inﬂuence of clamping t o be weaker for real size industrial panels that have much bigger radii of curvat ure. 3.3 The inverse problem solver validation The results of the inverse problem solver validation using t he 2024-T3 plates are presented in Figure 14. The fact that these plates were peened from only onesidereducedtheeﬀectoftheundesiredprestress.Never theless,theoutlines of the treated segments were blurred by the wide shot stream, and clamping of the plates altered their curvature, which caused discrep ancy between the simulations and the scanned shapes. The discrepancy between the simulated shapes and the target shapes showed that the ﬁxed peening regimes did not suit perfectly f or the chosen target shapes. Thus, the simulated Ω zequaled 1.2 % for the upper part and 1.4 % for the lower part with the free pattern, that was adjust ed iteratively and was free of constraints, and grew up to 5.5 % and 12.9 % resp ectively with the grouped pattern, that implied only the ﬁxed peening regi mes. The main reason for this is that the precise reproduction of the targe t curvature required &'( ') * + , . /de treatment)7 + *8eatment 7 +7 )8 n : ; n < * /= n .> 8: ? < n ) top side treatment8 n : ; n < * /= n .> 8: ? < n )7 + *8eatment Fig. 13 Diﬀerence between the prescribed peening pattern and the pattern that was act u- ally applied during treatment. a) The treated area prescribed by the pat tern #3 on the top side. b) A photo of the top side of the plate that was treated followi ng pattern #3. The contour of the treated segment is smoother than the one prescribed by the pa ttern, because the shot stream was too wide to reproduce the contour precisely. Springer Nature 2021 L ATEX template 24Simulation and automation of aluminium panel shot peen form ing U @ @ A B C D E GL H I A B C D E G TJ K M N O P Q J R N P a)b) S T V W XJ O N YP Q J R N P c)d) =5.5% Z [\ ] ^^Z [ _ _ ^^[ ` f i j kSl J m m N YP Q J R N P f)[ ` ` k[ ` o k T p q r s u v q v v w s w x p w y z { | }~ p   p  r s u v q v v w s w x p w y zPeening regimes{ | { | { | { w p q v v w s w x e)Ω Ω Ω Ω Fig. 14 Results of the inverse problem solver validation. The target shapes const itute upper (a) and lower (b) parts of the airplane wing skin reproducing the NACA 5430 airfoil. The inverse problem was solved for both shapes. We simulated deformation induced by the prescribed patterns for the upper (c) and lower (d) skins and compared the simulated shapes with the target shapes. The aluminum plates were treated according to th e patterns, and the 3D scans of the upper (e) and lower (f) skins were compared to the targ et shapes. The deﬂection ∆ zis computed as a function of the target shape. The parameter Ω zdescribes the dimensionless error with respect to the target shape (Eqn. 13). The plates were made of aluminum 2024-T3 and had the size of 320 ×320×1.6 mm. They underwent uniform treatment from both sides with regime #6 before application of the peen ing patterns. gradual transition of intensities from the front line of the wing towards the back line. However, having four predeﬁned peening regimes w as not suﬃcient to ensure smoothness of this transition. Moreover, the high est required eigen- strain for the upper (0.0023) and for the lower (0.0016) part s was higher than that provided by the most intense regime #2 (0.0015). Conseq uently, the cur- vature induced by the grouped pattern along the front line wa s inferior to the target one. Figure 15provides details on the free and the grouped patterns for both target shapes. Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing25 The upper plate showed a bistable behaviour after treatment . Thus, the computed ﬁnal shape appeared to be an unstable state of equil ibrium, and the plate tended to curve along the diagonal axis in an unconstra ined state. For this reason, we applied a weight on two corners of this plate d uring scanning, which is a common practice during the accuracy measurements in the peen forming industry. The weight the plate to keep the prescribe d longitudinal symmetry, but the two corners had to be discarded from the num erical model provided by the 3D scanner (see Figure 14). We used the two plates to assemble a physical model of the NACA 5430 airfoil, which is presented in Figure 16. The forward and the back parts of the model wing, which are colored in black in Figure 10, were printed in 3D using an FDM printer charged with the PLA plastic. These parts also served as supports for the aluminum plates. The plates were drilled an d assembled with the printed parts using binding barrels and screws. The cons traints induced by the binding have eliminated the shaping error, so the plat es perfectly ﬁt the NACA proﬁle. This situation is often the case in industry , where the parts are shaped up to a pre-deﬁned tolerance, and the error is enti rely eliminated only when the part is installed.                        a)b)            c)d) ∗ ¡ ¢ £Peening regimes¡ ¢¤¡ ¢¥¡ ¢ ¦¡ § ¨ © ª ª § « § ¬00.0023 Fig. 15 The eigenstrain patterns prescribed by the inverse problem to shape the ai rplane wing skin parts. Before grouping, the patterns assign gradual transi tion of the eigenstrain along the plate surface to shape the upper (a) and the lower (b) skins. The grouping process divides the patterns for the upper (c) and lower (d) parts into segment s with constant eigenstrain. The eigenstrain in the grouped segments corresponds to the on e induced by the ﬁxed peening regimes. Springer Nature 2021 L ATEX template 26Simulation and automation of aluminium panel shot peen form ing Fig. 16 The model of an airplane wing skin shaped as the NACA 5430 airfoil w ith a chord length of 500mm. The upper and lower skins were shaped using the shot peen f orming, and the forward and the back parts (black) were printed in 3D. The str ucture is assembled using binding barrels and screws. 4 Conclusion We developed a complete workﬂow for the shot peen forming aut omation. We generate the numerical peening pattern and we automaticall y translate it into a program for the shot peening robot. The program admits appl ication of multiple peening regimes following complex patterns and th us allows to shape complex doubly-curved shapes. The workﬂow implies prelimi nary characteri- zation of the shot stream eﬀect in terms of eigenstrain and ch aracterization of the eigenstrain anisotropy. These preliminary stages are n ecessary to generate the peening pattern and the nozzle paths. Simple yet eﬃcient procedures were designed for this purpose. The workﬂow was validated experimentally. The validation i mplied treat- ment following randomly generated patterns and shaping ski ns of a model airplane wing. The treatment was performed at high Almen int ensities of 10.9- 23.5A to reveal all possible errors appearing with large def ormations. The validation showed promising results but revealed an error b etween the simu- lated shapes and the practically developed shapes from 8% to 28.5% of the total deﬂection. We expect this error to be smaller for indus trial parts that have bigger radii of curvature and include stiﬀeners, becau se these compo- nents are less subjected to the undesirable prestress and to the inﬂuence of clamping supports. The inverse problem solver validation a lso emphasized the importance of an optimal choice of the peening regimes for ea ch target shape. The vast experimental campaign indicated directions for th e workﬂow improvement. The simulation error was provoked by an incomp lete consistency between the simulated and the real peening conditions. The m ajor diﬀerence lies in the fact that the software simulates immediate appli cation of the peen- ing treatment from both sides, while in practice it was ﬁrst f ully applied from one side and then from the other. The future experiments will imply uniform application of the pattern from both sides, i.e., alternati on of the peened sides after each pass of the peening nozzle. In addition, the peeni ng nozzle will be modiﬁed in order to provide a narrower shot stream and thus to reproduce the pattern more accurately. The eﬀect of clamping on the eigens train distribution will also be quantiﬁed. Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing27 Acknowledgments. The authors gratefully acknowledge ﬁnancial support from The Natural Sciences and Engineering Research Council of Canada (NSERC), from The Fonds de Recherche du Qu´ ebec – Nature et Te chnologies (FRQNT)andTheMinist` eredel’ ´Economieetdel’InnovationduQu´ ebec,from The Aluminium Research Centre – REGAL, and from the industri al partner of this project – Aerosphere Inc. The authors also acknowled ge the consult- ing support provided by Aerosphere Inc. An important contri bution to this project was made by Rapha¨ el Paradis, who fabricated the pee n forming ﬁx- ture during his internship at the Aerospace Technology Cent er (CTA), and by Nada Abbassi, who performed the experimental manipulation s for validation of the inverse problem solver during her internship at Polyt echnique Montreal. Declarations Funding and Conﬂicts of interests/Competing interests This study was funded by: •The Natural Sciences and Engineering Research Council of Ca nada (NSERC) (funding reference numbers RGPIN-2019-07072 and R GPIN- 06412-2016); •The Fonds de Recherche du Qu´ ebec – Nature et Technologies (F RQNT) (funding reference number 0000295154); •The Fonds de Recherche du Qu´ ebec – Nature et Technologies (F RQNT) and The Minist` ere de l’ ´Economie et de l’Innovation du Qu´ ebec (funding reference number LU-210888); •The Aluminium Research Centre – REGAL; •Aerosphere Inc. The authors have no relevant ﬁnancial or non-ﬁnancial inter ests to disclose. The authors have no conﬂicts of interests/competing intere sts to declare that are relevant to the content of this article. Appendix A Proportionality between the dimensionless load and curvature We performed a speciﬁc series of simulations in order to ﬁnd t he proportion- ality constant µrelatingκhand ΓB(Eqn.15). Namely, we took a range of the dimensionless loads Γ Band simulated its application on aluminum coupons of size 76 ×19 mm and various thickness h. All in all, 100 testcases for Γ B ranging between 0 and 0.008 were simulated. For each testcas e,εt ∗was positive andεb ∗was negative, but the ratio εt ∗/εb ∗was chosen randomly between -1 and 0. The material was considered plastically isotropic. Addi tionally, the thick- nesshranged randomly from 2 to 4 mm for each testcase. The Poisson’ s ratio νwas ﬁxed equal to 0.33, that is a typical value for aluminum. W e meshed the coupons with 760 triangular elements and solved the forw ard problem as described in Section 2.1.3. For all the testcases, the curvature κmeasured at Springer Nature 2021 L ATEX template 28Simulation and automation of aluminium panel shot peen form ing the center of the coupon varied by less than 0.015% in diﬀeren t directions that proves the results presented in Ref. [ 26]. Appendix B Inﬂuence of the eigenstrain anisotropy Theeigenstrainanisotropyismajorlycausedbythemateria lproperties,butits inﬂuence is not always perceptible during uniform treatmen t of the 76 ×19 mm coupons. Thus, the 1.6 mm and 2.06 mm thick coupons develop th e same local curvature in both directions when treated with the reg imes described in Table1. This is explained by the fact that the curvature of uniforml y peened rectangular specimens depends on [ 25] [26]: ΓB= ΓBD2 h2, (B1) whereDis the characteristic length of the specimen. Namely, recta ngular specimens adopt spherical shape for low values of ΓBindependently of the eigenstrain anisotropy. Thus, the applied treatment was no t intense enough to make the 76 ×19 mm sized coupons develop diﬀerent curvatures in the L- and T-directions. The 320 ×320mm sized plates have 17.7 times bigger ΓBthan the coupons of the same thickness treated with the same regime. C onsequently, this type of plates was used for determination of the eigenst rain anisotropy. Appendix C Inﬂuence of pre-treatment on the induced eigenstrain We performed additional measurements to characterize the e igenstrains induced in the 2024-T3 panels that did not undergo the prelim inary double- side treatment with regime #6. First, we measured the eigens trains induced in the 76×19 mm coupons cut of the same plates (see Section 2.5.2). The results of this measurement are summarized in Table C1. The induced deﬂection of the coupons was bigger than in the case with pre-treatment, s o the induced eigenstrains were bigger as well. This phenomenon is explai ned by hardening of the surface generated by the pre-treatment, which is disc ussed for example in Ref. [46]. The inﬂuence of hardening was more signiﬁcant for less int ense peening regimes (Table 1). Thus, the absence of pre-treatment increased the eigenstrain induced by the regime #5 by 169%, while the eigen strain induced by the regime #2 was only increased by 30%. The graphical comp arison of the eigenstrain induced in both cases is presented in ﬁgure C1. Next, we characterized the eigenstrain anisotropy using th e ”wave” exper- iment (see Section 2.5.3). Thus, we treated two 2024-T3 panels following the ”wave” pattern with regime #2 without pre-treatment. Th e experiment revealed a mild eigenstrain anisotropy, which is character ized in terms of the coeﬃcientχ∗in Table C2. The principal eigenstrain direction has changed Springer Nature 2021 L ATEX template Simulation and automation of aluminium panel shot peen form ing29 with respect to the case with pre-treatment, but the value of χ∗is of the order of 10−3, which means that the eigenstrain is almost isotropic. Give n that the eigenstrain anisotropy demonstrated by the pre-treated pl ates was mild as well (Table3), we conclude that the pre-treatment did not have a tangible inﬂuence on this parameter. Regime uavg(mm) σu(mm) εt ∗×103Increase in εt ∗ #2 0.706 0.008 1.91 169% #3 0.629 0.015 1.7 92% #4 0.582 0.012 1.58 42% #5 0.467 0.001 1.26 30% Table C1 Results of the additional eigenstrain measurements performed using the bunches of four 76 ×19×1.6 mm sized aluminum coupons (see Section 2.5.2). The coupons were made of aluminum 2024-T3 and did not undergo any pre- treatment. They were ﬁxed during peening with an Almen holder, and the deﬂection was measu red with an Almen gage. The average deﬂection in each bunch of four coupons is denot ed byuavg, and the parameter σudescribes the standard deviation in deﬂection. The eigenstrain induced in the top layer is denoted by εt ∗. The right column shows the increase in εt ∗with respect to the pre-treated coupons peened with the same regimes (Table 2).­eening regime ® ¯ ® ° ® ± ® ²00.5 ³ ³ ´ ²2 ∗× 103 µ ¶ · ¸ ¹ º »e-tr¼ ½ ¸ ¾¼ ¿ ¸ µ ¶ · ¸ ¹ À Á ¸ º » ¼ µ ¸ » ¼ ½ ¸ ¾¼ ¿ ¸ Fig. C1 Graphical comparison of the eigenstrain εt ∗induced by the same peening regimes on pre-treated and not pre-treated coupons. The coupons had the size of 76 ×19×1.6 mm and were made of aluminum 2024-T3. The pre-treatment was done usin g the regime #6. The parameters describing each peening regime are presented in Table 1. The exact values ofεt ∗for each case are presented in Tables 2andC1. χLχTχ∗ -0.013 0.003 0.008 Table C2 The anisotropy coeﬃcients determined using the “wave” experiment for the 2024-T3 plates that were not pre-treated. The “wave” pattern was app lied using the regime #2. The eigenstrain anisotropy coeﬃcient χ∗is the average of the absolute values ofχLandχT. Springer Nature 2021 L ATEX template 30Simulation and automation of aluminium panel shot peen form ing