Electro-aeromechanical modelling of actuated membrane wings

This paper presents a numerical investigation on the aeromechanical performance of dynamically actuated membrane wings made of dielectric elastomers. They combine the advantages of membrane shape adaptability, which produces increased lift and delayed stall, with the benefits of simple, lightweight but high-authority control mechanism offered by integral actuation. High-fidelity numerical models have been developed to predict their performance and include a fluid solver based on the direct numerical integration of the unsteady Navier–Stokes equations, an electromechanical constitutive material model and a non-linear membrane structural model. Numerical results show that harmonic actuation can either increase or reduce the overall aerodynamic efficiency of the wing, measured as the mean lift-to-drag ratio, depending on the ratio between the actuation frequency and the natural frequency of the membrane. In addition, the definition of a reduced-order model based on POD modes of the complete high-fidelity system provides an insight of the main characteristics of the dynamics of the coupled system. & 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Wind tunnel testing of membrane wings (Winter and Von Helversen, 1998;Newman, 1987;Rojratsirikul et al., 2009) has shown that they can offer a superior aerodynamic performance in low Reynolds number flows. This includes delayed stall, higher lift, and enhanced manoeuvrability and agility as compared with their rigid counterparts. These aerodynamic advantages are mostly related to the membrane compliance that allows the passive shape adaptation to the pressure gradients and energy-efficient coupling with the unsteady flow phenomena (Rojratsirikul et al., 2010). However, membrane wings also require a supporting frame, which in general will degrade their performance. The shape and the stiffness of the supports strongly influence the dynamics of the wing (Arbos-Torrent et al., 2013). With that as a starting point, further gains would be possible through wing actuation. While this has mostly been attempted using mechanical actuation (Stanford et al., 2008), integral electromechanical actuation could provide a much higher control flexibility. In particular, the use of dielectric elastomers (DEs) as material for the wing would allow the coupling of the advantages of the passive adaptation offered by the wing compliance, with the benefits of a simple, lightweight control mechanism from embedded actuation. DEs are composed of a thin polymeric layer sandwiched between two compliant electrodes. The application of a voltage through the thickness causes the compression of the material in that direction and its in-plane extension (Suo, 2010). This defines a variation of the tension in the wing that can be used to tune its dynamic behaviour.
( -), separation and recirculation are the main sources of unsteadiness in the interaction with the membrane, and cause the strong coupling between the fluid and the membrane vibrations (Rojratsirikul et al., 2009(Rojratsirikul et al., , 2010(Rojratsirikul et al., , 2011. On the structural side, Arbos-Torrent et al. (2013) have experimentally shown that mean camber and modes of membrane aerofoils are strongly influenced by the size and shape of the leading-and trailing-edge supports, which in practice, can never be neglected. The membrane compliance determines large structural displacements and the soft material used is usually accompanied by a rate-dependent constitutive behaviour. Rojratsirikul et al. (2011) Bleischwitz et al. (2015) investigated the effects of the wing finite aspect ratio on the aerodynamic performance, identifying dominant spanwise vibrations generated by the interaction of the tip vortices with the membrane structure.
To the authors' knowledge only Hays et al. (2013) and Curet et al. (2014) have to date experimentally tested the use of DEs for integral actuation of membrane wings. From the computational point of view, the only high-fidelity models in the literature for actuated membrane wings consider only two of the three fields involved and the problem (see Fig. 1, right). The goal of this paper is to investigate a comprehensive simulation approach including the aerodynamic, structural and electro-mechanical aspects for the design and performance evaluation of actuated DE membrane wings.
Earlier numerical models of passive flexible membrane wings are based on linear elastic structural models coupled with an inviscid (Newman, 1987) or viscous (Smith and Shyy, 1995) flow description. Tiomkin et al. (2011) andGordnier (2009) investigated the influence of membrane rigidity, pretension, angle of attack and Reynolds number for a 2D membrane in laminar flow ( Re ¼2500-10 000) for static and dynamic cases. When considering DEs, a comparison with experimental results indicates that a geometric non-linear structural model and hyper-elastic constitutive behaviour need to be considered instead of linear models (Wissler and Mazza, 2007;Mokarram et al., 2012). Lian et al. (2003) and Lian and Shyy (2007) considered a hyper-elastic nonlinear constitutive model, to deal with the large structural displacements, with a e N transition model for the fluid, to predict the separation point. From the fluid point of view, separation, recirculation and vortices interactions require the use of transition models, if RANS models are adopted, opportunely tuned LES models, or DNS simulations.
All previous high-fidelity models are required for simulations, but a simplified description could be beneficial, and after even mandatory, in the context of parametric investigations and control system design. Several model reduction techniques have been proposed in the literature for aeroelastic systems, including Volterra series (Raveh, 2001), proper orthogonal decomposition (POD) (Feng and Soulaïmani, 2007), balanced truncation and co-prime factorisation (Hesse and Palacios, 2014). POD has also been widely used in the postprocessing of results for unsteady aerodynamics (Rowley, 2005), structural dynamics (Kerschen et al., 2005) and, in a more limited way, fluid-structure interaction problems (Attar and Dowell, 2005;Lieu et al., 2006;Liberge and Hamdouni, 2010). Some efforts have already been done in the use of POD to analyse the flow field around a membrane wing . Linearity is often assumed (Khalil et al., 2007) in this case but some works have proposed parameter dependent functions to take into account the variation of the system properties with its states (Attar and Dowell, 2005). This work will investigate the potential gain in performance of dynamically actuated membrane wings in low Reynolds number flows and the challenges associated to full-system simulation. Section 2.1 presents a non-linear elastic constitutive model for the DE, with a linear electrostatic model for the Maxwell stresses. This has shown to be a good approximation to the material behaviour in practical situations, but the model is also capable of dealing with electrostrictive effects if necessary. This is the first predictive model of the static performance of inflated VHB4905 membranes. In Section 2.2 the structural model is coupled using an implicit time marching algorithm with non-matching spatial discretisations, with a flow solver based on the direct numerical integration of the Navier-Stokes equations. In Section 3.2 the proposed dynamic aeroelastic model will be compared with numerical results from the literature to verify the implemented time-marching simulation method. Section 3.3 explores the variation of the aerodynamic performance obtained with the actuation of the membrane wing. In Section 3.4, POD of the full aeroelastic system is used for model reduction and system identification to estimate the membrane shape as a function of the actuation.

Actuated membrane model
It is usual practice in the literature on membrane wing modelling (Tiomkin et al., 2011;Gordnier, 2009) to assume ideal membrane behaviour for the structural model, neglecting the thickness effects because only passive membranes were considered. This is however not sufficient for integrally actuated membranes where the electrostatic stress tensor needs to be defined in all three dimensions. Commonly the effect of the Maxwell stress is modelled as an equivalent through-thethickness compressive stress or an in-plane relaxation stress (Suo, 2010), using the incompressibility assumption. For high voltages this might not be a valid assumption, requiring the computation of the real stress tensor. This model will be solved using the solid eight-nodes three-dimensional hybrid elements with translational degrees of freedom in Abaqus (Dassault System, 2013). Their formulation reduces the volumetric locking that can be encountered with incompressible or nearly incompressible materials, and it is hence suitable for hyperelastic materials (Dassault System, 2013). Comparisons with membrane elements showed the absence of volumetric locking for the range of deformations considered. Large structural displacements are required to consider a non-linear geometric structural description of the membrane as well as a nonlinear constitutive model. Those are described next in some detail to justify the modelling assumption choices used in the definition of the material response. To our knowledge there are no predictive models in the open literature of the real performance of inflated DE membranes capable of considering boundary conditions found in wing applications.

Hyperelastic material model for DEs
The constitutive model for DE is developed in the finite deformation framework (Bonet, 2001;Bonet and Wood, 2008;Simo and Hughes, 1998). Let 0 be the reference configuration of a body at time t 0 and ϕ the mapping function from 0 to the current configuration at time t (Fig. 2). Let F be the deformation gradient of ϕ expressed in the reference configuration and consider its multiplicative decomposition (Flory, 1961) into its volumetric and isochoric components, that is, where F measures the isochoric deformation of the body and J F det = ( ) measures changes in volume. F and J are used to define the elastic deviatoric and hydrostatic stresses in the deformed configuration, respectively.
In this work, an isotropic constitutive model is assumed. This has been experimentally verified on DEs up to a stretch of around 1.5 (Schmidt et al., 2011), but would need further investigation for higher stretches. The mechanical free-energy function of an isotropic material can be expressed in terms of the invariants of the right Cauchy-Green deformation tensors C F F T =¯¯ (Bonet and Wood, 2008): Thus the free-energy function is of the form where U J ( ) is the volumetric energy function accounting for the variation in volume defined by J and I I , is the deviatoric free energy function.
In this work, the elastic constitutive law for DE will be the nearly incompressible model for the deviatoric part due to (Gent, 1996) where μ and J m are the elastic shear modulus and the limiting stretch of the material, respectively. This model is used because it consists of two parameters of immediate physical significance and it predicts well the stiffening behaviour characteristic of membranes when the deformation approaches the value of the limiting stretch. Both constants are determined from experiments. The volumetric energy function is expressed by where K is the compressibility modulus of the material. The value of the material constant, as in the deviatoric part, has to be determined experimentally. The stresses are the conjugates with respect to the free energy function of the stretches in the directions considered (Bonet and Wood, 2008). The second Piola-Kirchhoff stress tensor S represents the stresses expressed in the reference configuration and is obtained as where S J are the stresses due to the contribution of the volumetric energy function, S ∞ is the contribution of the elastic deviatoric stresses and the projection tensor is   with  being the fourth-order identity tensor. The mechanical stress tensor in the deformed configuration m σ is finally obtained through a push-forward operation on S: The constitutive model is completed by adding the effect of the electrostatic stresses using Maxwell's stress tensor. Given an electric field vector E, the corresponding stress tensor e σ is defined as (Suo, 2010) where ϵ is the material dielectric constant, E is the electric field vector in the deformed coordinate system and I is the identity second order tensor. The material dielectric constant ϵ is usually expressed as r 0 ϵ = ϵ ϵ where 0 ϵ is the vacuum dielectric constant and r ϵ is the material relative dielectric constant. The total stress is finally obtained as σ (Suo, 2010): . 9 m e σ σ σ = + ( ) These constitutive equations have been implemented in the finite element (FE) package Abaqus (Dassault System, 2013). An implicit iterative scheme is used to obtain the stresses in terms of C at the end of each time increment. The equilibrium system is solved using an iterative Newton-Raphson scheme, with the tangent stiffness matrix  in the reference configuration obtained as (Suchocki, 2011) (6). The exact tensor is computed analytically and finally written in the deformed coordinate system by means of a push-forward operation The tangent stiffness matrix is recomputed at every iteration considering the actual geometry and stresses in the material.

Material coefficients
As it was already discussed, the material coefficients for the free-energy function (4) and the Maxwell stress tensor (8) need to be experimentally determined. Fox andGoulbourne (2008, 2009) reported the performance of an actuated inflated membrane made of acrylic dielectric elastomer (VHB4905 from 3M), that can be representative of the behaviour of membrane wings. The static data from their work will be used to identify the material coefficients of the constitutive model described above. An analytical membrane model will be used in the fitting process, which will be based on the least-square approach. The model uses the spherical cap assumption of the deformed shape, which is valid when the maximum height of the central membrane point is lower than its in-plane radius (Shian et al., 2010). With this assumption we have that where z is the displacement of the central point of the membrane, R is the radius of the correspondent spherical cap and a is the initial membrane radius. Combining (11) with the Young-Laplace law allows us to relate the height of the spherical cap z with the in-plane tension T and the pressure difference across the membrane P Δ : The tension needs to be recomputed iteratively considering the material model selected and the updated value of the stretch due to the out-of-plane displacement that is obtained from where λ and p λ are the final and the initial membrane stretch, respectively. When the membrane is actuated, the total inplane stress is the sum of the mechanical and Maxwell stresses.
The geometrical data used in the model to reproduce the experimental data on VHB4905 from Fox and Goulbourne (2009) are initial membrane thickness h 0.5 mm = , membrane radius a 88.9 mm = and p λ ¼3.5. The coefficients obtained were found to be close to J 20 kPa, m μ = ¼100, K 3.5 10 6 = × Pa and r ϵ ¼ 2.7. Those are the values that will be used in this work to represent the DE material behaviour. The comparison of the present FE model with experiments is showed in Fig. 3, demonstrating a very good agreement. In addition, a comparison of the results of the FE model with the non-linear analytical model used to identify the coefficients further confirmed the absence of volumetric locking with hybrid solid elements used to model the membrane. The maximum voltage considered is 2.5 kV in the experiments, corresponding to a nominal applied electric field of 61 MV/m, which set the upper limit of validity of the model that will be considered in this work. Similar values of the shear modulus have been proposed in the literature (Fox and Goulbourne, 2009;Schmidt et al., 2011) using an Ogden material model. (Lu et al., 2012) proposed a higher range of values of the shear modulus for the Gent model and a similar J m when fitting equibiaxial stretch experiments. The discrepancy is determined by the different boundary conditions since no out-of-plane loadings were considered in their work. In addition, the static configuration considered there is not fully relaxed, and includes some effects from the viscoelastic stresses. The same methodology for coefficient identification is applied to the experimental data available the VHB4910 membrane (Fox, 2007) and the material coefficients obtained are J K 24 kPa, 90, 3.8 10 m 6 μ = = = × Pa and r ϵ ¼4.8.

Aeroelastic model
A viscous flow solver is used to compute the low-Reynolds flow around the membrane. As opposed to the structural model, that required a number of modifications to the conventional theory, the direct solution of Navier-Stokes implemented in STAR-CCMþ (CD-Adapco, 2013) will be used here. The solution is based on a finite-volume discretisation of the compressible equations. A coupled approach is used to solve the flow and energy equations, which requires more memory but offers better stability and velocity in the convergence. Second-order accurate integration in space and time are used. The maximum number of sub-iteration for each time-step is set to 10 and the limit relative tolerance for the monitored flow properties is set to 10 4 − . The fluid domain considered is shown in Fig. 4 together with a zoom on the wing mounted on rigid supports. The leading and the trailing edge are modelled as rigid, aerodynamically shaped supports. Both supports and the membrane are characterised by a no-slip wall boundary condition. The outer boundaries of the domain extend for 100 chords length in every direction. At the domain inlet the velocity magnitude and direction are set as boundary conditions, while at the outlet a constant pressure is specified. The other two boundaries are modelled as free-stream boundaries, where the direction of the velocity, the Mach number and the pressure are specified. In case a 2D problem is considered, periodic boundary conditions are imposed in the lateral surfaces. The domain is discretised using a structured mesh whose size ranges from 0.25 10 6 × to 2 10 6 × elements. For each case considered, a mesh-refinement and time-step sensitivity study has been conducted. For the dynamic cases, the relative errors on the mean value of C l and C d as well as the amplitudes and frequencies of their oscillations were considered to establish convergence. The maximum relative error of those quantities was set to be below 2%. In addition to the sensitivity to the grid-resolution and time step, for each case considered the element size and time step have been compared with the characteristic viscous length and time scales, respectively (Valen-Sendstad et al., 2011): where ν is the kinematic viscosity of the fluid and u ν is the friction velocity defined as (Pope, 2000) u SS , 15 ij ij average of the elements size, the characteristic length of the mesh elements is l l l l mesh x y z 3 = with l l , x y and l z being the dimensions of the hexahedral cells. The correct spatial and temporal resolution of the smallest turbulence scales requires (Pope, 2000) l l l 1 , 16a where t Δ is the time step of the implicit solvers.

Fluid-structure coupling
The coupling of viscous flow and structural solvers used in this work requires the definition of a common interface. It has to occupy the same spatial positions in the models, allowing the mapping of the nodes of the elements of the two meshes considered. Once neighbouring nodes and elements are identified, the fluid solver uses the shape functions of the structural solver for the interpolation of the pressure and viscous forces and the definition of the equivalent set of nodal loads to be used in the finite element solver. The new nodal displacements and velocities from the FE solver are then passed to the fluid solver for the deformation of the mesh and the computation of the new flow field. The mesh deformation process moves the mesh nodes according to a linear interpolation law that is depending on their relative distance from the moving and fixed boundaries. An implicit coupling was found necessary, considering the strong interaction of the fluid and the membrane structure in this problem (Fig. 1, right). When the relative difference in the norm of two consecutive field exchanges is below a defined tolerance, 10 4 − for the results in this work, convergence is established. The exchange of data between the solvers uses a bridge in the RAM memory (CD-Adapco, 2013).

Model reduction
Once the simulation environment has been outlined, it remains to describe the process to generate reduced-order models used to identify the key dynamic features of the response to actuation (and subsequently, if needed, although this is not done in this work, to design a control scheme). Proper orthogonal decomposition (POD), also known also as discrete Karhunen-Loéve decomposition (Karhunen, 1946;Loeve, 1955), is used here to reduce the size of the aeroelastic system and facilitate control design. This model reduction is very appealing in the context of system identification, because it significantly reduces the number of variables that have to be determined, but it may have very slow convergence and only local validity. The time history of selected state variables of the aeroelastic system is obtained from the high-fidelity simulation. The snapshot matrix, S, whose columns contain the values of the system variables at each time step, is assembled and used for the calculation of the correlation matrix, C S S T = (Buljak, 2012). The eigenvectors of C define the basis of the projection matrix ϕ that is used to project the time history of the state variables and obtain the time evolution of the amplitude of each base considered. The amplitudes are the new variables of the reduced-order system. Their number depends on the number of eigenvectors of C retained. The choice is done on the basis of the convergence to the high-fidelity solution and the percentage of the total energy retained (Buljak, 2012). In this work it is set to be greater than 99.5%.

Numerical results
In order to explore the key dynamic characteristics of systems showing full fluid-electrical-mechanical coupling, the numerical results in this paper have been limited to 2D membranes. First of all, the constitutive model of Section 2.1.2 will be used to reproduce the static experimental results on a dielectric tunable lens from Shian et al. (2010). This is done to prove that the material model is predictive of the behaviour of the DE considered. Results with a linear elastic constitutive model will be then compared with those of Gordnier (2009), to verify the correct coupling of the solvers. Then, the constitutive model of the DE described in Section 2.1.1 will be used to investigate the effects of the integral actuation on the aerodynamic performance of the actuated wing with the same geometry. In the final part, the proper orthogonal decomposition will be used to define a reduced-order model of the full aeroelastic system and identify the fundamental dynamic characteristics of the coupled system.

Electromechanical model verification: tunable lens model
The material models for DE described in Section 2.1, with the material constant for VHB4905 and VHB4910 obtained in Section 2.1.2, are here used to reproduce the experimental results of (Shian et al., 2010) on a tunable lens and demonstrate the predictability of the constitutive models. The lens is composed of two membranes made of transparent DE membranes and mounted on a rigid frame (Fig. 5). A transparent fluid is encapsulated between the membranes. The passive membrane is made of VHB4905 while the active one is made of VHB4910. When the active membrane is actuated, the in-plane tension is reduced causing a redistribution of the enclosed fluid and a variation of the focal length. For a detailed description of the experiment the reader is referred to the cited paper.
A comparison of the static model with experimental data from Shian et al. (2010) is presented in Fig. 6. A very good agreement is obtained up to 4.5 kV. The discrepancy at the higher voltages can be attributed to the assumption of a constant dielectric value for the material that neglects electro-strictive effects.

Aeroelastic model verification: passive membrane wings
Next, we will verify the ability of the current approach to capture the unsteady aerodynamics on passive membranes mounted on rigid supports. The support geometry is from the experimental setup by Rojratsirikul et al. (2009). The chord is c 0.1366 m = , the thickness is h 0.2 mm = , and the membrane is mounted on the supports with zero prestrain 0 δ ¼0. The material is modelled with a linear elastic constitutive law with Young's modulus E 3346 Pa = , Poisson ratio ν ¼ 0 and density 473 kg m     (2009), who assumed a membrane structural model with constant and uniform thickness and only out-of-plane displacements. The flow field is solved there using a sixth-order Navier-Stokes solver, allowing for a coarser mesh than the one used in this work. Fig. 8 shows a very good agreement of the mean deformed shapes and mean C p distributions on the membrane for both cases considered.
Numerical results showed convergence to very similar values to those obtained by Gordnier (2009) and the small differences can be attributed to thickness effects (which are only included in the current model) that also led to some differences in the shape of the supports. In fact, from Fig. 8 the main differences in the results are in proximity to the leadingedge support (x ¼0). The scattering of the points for C p at 8 α =°at x c / 0 = is due to the sharp edge that is generated by the membrane deformation at the attachment with the leading-edge support. Increasing the incidence angle, the effect of wake instabilities increases in importance causing the membrane to oscillate over the mean configuration with increasing amplitudes. When the angle of attack is 8°, the system exhibits a self-excited oscillatory behaviour, with a dominant third membrane mode with a nondimensional frequency of St¼ 1.45, as reported by Gordnier (2009). This comparison demonstrated the suitability of the current solution method to represent the dynamics of the membrane wings at low Re and moderate angles of attack.

Membrane wing actuation
The electromecanical material model described in Section 2.1.1 is here coupled with the aeroelastic model from Section 2 to investigate the expected performance of an integrally actuated 2D membrane wing made with VHB4905 elastomer. Curet et al. (2014) experimentally investigated actuated DE membrane wings made with a similar material, the acrylic VHB4910 elastomer, and found a strong dependency of the aerodynamic performance with the actuation frequency. With excitations near natural frequencies of the system, the wing exhibits the maximum increase in the averaged values of the lift and drag coefficients, due to the enhanced camber effects. The maximum increase in efficiency, though, is not obtained for that value of voltage frequency, but for a lower one because the increase of C l close to the resonance is gradual, while the C d manifests a sharper peak.
In the proposed numerical model, it is considered a Re number lower than the ones from Curet et al. (2014), 5.7 10 4 × and 8 10 4 × . This is due to two reasons: Firstly, the passive aeroelastic model has been verified for Re¼2500 in the previous section, and those results will serve as starting point for the investigation of the actuated membrane. In addition, with Re ¼2500, the flow is laminar and can be modelled using a 2D-assumption. This reduces the size of the model that has to be solved and allows for a faster parametric investigation of the actuated aeroelastic system. For these reasons, the same flow conditions and geometric parameters of Section 3.2 are used with Re¼2500 and 4 α =°. The Mach number is 10 4 − . The initial thickness is set to h 0.05 mm 0 = and the membrane prestretch is set to 1.02, to allow for a high compliance of the wing and benefit from the passive wing aerodynamic performance. The value of prestretch is smaller than typical prestretch levels in DEs and it has been chosen to allow for large structural deformations and test both the structural and fluid solvers. It has to be considered that the resulting thickness of the prestretched wing is compatible with the usual thickness of DE actuators, μ as well as the through-thickness elements, and the bigger elements, at the domain boundaries, are around twice the size of the membrane chord. In the passive case, the equilibrium configuration presents a steady camber with maximum amplitude of 2% at x c / 0.44 = , and a lift coefficient C 0.47 l = . The flow induced deformation determines an increase of the membrane stress of 3.4%. Considering as reference the deformed configuration, a linear analysis predicts a natural frequency in vacuum at St¼0.55. A spectral analysis of the system shows that, because of the effect of the surrounding flow, the first natural frequency of the aerolastic system is reduced to St¼0.37.
The mean deformed shapes and aerodynamic coefficients of the passive case are considered as reference values for the dynamic actuation of the membrane using a through-the-thickness voltage. The amplitude of the harmonic voltage is fixed to V 500 V 0 = , generating a nominal electric field similar in magnitude to the one used by Curet et al. (2014) and well within the limit of 61 MV/m from Section 2.1.2. The actuation frequency, St v , is varied to investigate its effect on the aerodynamic performance of the wing. The resulting forcing of the system is the second power of the voltage, t 1 cos 4 St defining the Maxwell stress. The forcing input in the system evolves with a frequency that is twice the voltage frequency, St v . If the actuation frequency is set to zero, the input of the system has a constant value of . This is equivalent to an actuation with a constant voltage of V 2 0 , which is set to be the magnitude of the voltage when the static actuation case is considered. To reduce the high-frequency noise in the membrane response, the voltage is initially ramped with a 5th order polynomial. In addition, the first two actuation cycles are not considered in the results to remove the effects of the transient and focus on the periodic response of the wing. In both static and dynamic cases, the actuation determines a reduction of the mean membrane tension of approximately V h 2 0 2 ϵ , causing an increase in camber and hence in the mean lift and drag coefficients. The mean aerodynamic coefficients with actuation are presented in Fig. 9. The first observation is about the high increase in performance obtained with a static actuation that, because of the low increase in drag, appears to be beneficial for the overall performance of the wing. This is equivalent to considering a passive wing with a pretension that is reduced by the value of the resulting Maxwell stress. The increased compliance allows for higher lift, but does not allow for a dynamic tuning of the dynamics of the system, that is what is required to operate this kind of wings.
When the wing is dynamically actuated, there is an initial drop of the lift coefficient, which then increases as the frequency approaches the first resonance of the aeroelastic system. The drag coefficient, instead, keeps growing with the increase in the actuation frequency. The rate of variation of the lift in Fig. 9 is almost steady in the frequency range before the first natural frequency of the system. On the other hand, the C d increase is steeper and has a peak that is slightly delayed as compared with the C l one. This determines a peak in the wing efficiency (Fig. 10) at around St ¼0.12 and a significant drop Fig. 9. Variation of the aerodynamic coefficients with frequency, Re¼ 2500 and Fig. 10. Variation of the mean aerodynamic efficiency with frequency, Re¼ 2500 and for higher frequencies due to the increase in drag and drop in lift. It is interesting to understand what causes the significant drop in performance after the first resonance of the system. For this reason two cases will be taken in consideration, the actuation with a frequency of St 0.08 v = , presenting a positive, net increase in the wing efficiency and the St 0.16 v = , where the aerodynamic performance is deteriorated. Fig. 11 reports the time evolution of the aeroelastic system for a cycle of actuation with St 0.08 v = (two cycles in the response) starting at non-dimensional time t 2 = ⁎ . Fig. 11(a)-(f) shows the positive and negative vorticity contour plots at different times during the actuation cycle and Fig. 11(g) shows the evolution of the normalised aerodynamic coefficients C l and C d , membrane maximum amplitude, y a , and Maxwell stress, V 2 . The values are normalised with respect to the maximum value of the cycle. The time history is a function of the adimensional time t ⁎ . At the beginning of the cycle, Fig. 11(a), the membrane is near its minimum deformation where the displacement of the central point as well as the lift coefficient assume their lower values ( Fig. 11(g)). As the voltage increases, the inertia effects and tension relaxation due to the actuation determine the increase in camber, augmenting lift and drag. The maximum amplitude point moves downstream chordwise following the low pressure region determined by the leading-edge vortex travelling towards the trailing edge. At t 2.1 = ⁎ , Fig. 11(b), the vortex detaches at approximatively 60% of the chord, causing a sudden drop in the aerodynamic coefficient Fig. 11(c). The camber continues to increase, the maximum amplitude point moves upstream and the membrane finally assumes the convex shape typical of the first natural mode and the static equilibrium condition (Fig. 11(d)). Lift and drag coefficients increase in parallel to the deformation. The maximum amplitude point continues to move upstream and reaches its maximum, Fig 11(e). The maximum camber amplitude is delayed from the maximum value of the voltage due to the phase delay introduced by the damping of the system. The deformation starts to reduce back to the initial position in the cycle, Fig. 11(f). When the voltage gets to a value close to zero, the effect of the tension relaxation on the membrane is negligible, the rate of variation of the amplitude is reduced and due only to inertial effects. This explains the variation of the slope of the amplitude curve, generating a v-shaped path around t 2.5 = ⁎ in Fig. 11(g). A spectral analysis shows that the two dominant frequencies in the response are St¼0.16, twice the actuation frequency, St v , and the first natural frequency of the system, St¼0.31. A POD decomposition of the membrane deformation shows that the main modes excited are the first and second one. , the variation of the aerodynamic coefficients is mainly driven by the frequency of actuation, but a higher number of modes than the previous case are excited. A POD decomposition of the displacement history shows that amplitude of oscillations of modes from 3 to 10 are of the similar order of magnitude of the first two. The presence of these modes causes the reduction of the lift due to the local negative camber. In addition, higher inertia forces define a negative membrane amplitude that defines a broad valley in the C l history, reducing the time averaged value. These seems to be the main differences that cause such a different aerodynamic response. The numerical results are consistent with the experimental observations from Curet et al. (2014), despite the different flow regimes considered. In fact, in both cases, excitation at the resonance determines a degradation of the aerodynamic efficiency of the membrane. Actuation at a slightly lower frequency determines an overall positive aerodynamic enhancement.

System identification
A reduced-order model for the wing response presented in Section 3.3 is now obtained to identify the key dynamics of the system. The passive membrane in equilibrium is considered as reference and a uniform voltage V 500 V = is applied. The frequency is varied linearly up to St 0.2 v = in a time window of 240 t ⁎ with a time step of t 0.02 s Δ = ⁎ . This range is chosen to capture the two first resonant frequencies of the coupled system. The displacement history of the points in the suction surface of the membrane is stored and the vertical position of the points defines the system variables. The chord-wise displacement of the points is negligible, so their positions are assumed to be constant. In this work, only the structural  degrees of freedom are considered in the POD, to limit the number of data that has to be post-processed for the model reduction and allow for a preliminary assessment of the potentials of the POD for the aeroelastic problem of actuated membrane wings. The vertical coordinates of all the points of the membrane at every time step are used to build the snapshot matrix as described in Section 2.4. The number of basis is selected on the basis of the convergence of the reconstructed solution to the high-fidelity one. Fig. 12 shows that three basis capture very well the higher frequency content of the displacement history of the full solution. The basis selected are showed in Fig. 13(a)-(c), they account for the 99.8% of the total energy of the system (Fig. 13(d)) The modes plotted in Fig. 13 are similar to the classic first three eigenmodes in vacuum of an elastic string clamped at the extremities. They show some differences from the sinusoidal shapes due to the effect of the surrounding fluid that is captured thanks to the POD reconstruction.
The displacement, velocity and acceleration histories of the points are projected onto the 3-POD basis functions. The aeroelastic system is then modelled as a mass-normalised second-order linear system: where C ⁎ and K ⁎ are the damping and stiffness matrices, respectively, and G ⁎ is the gain matrix of the input signal V 2 , which is proportional to the Maxwell stress applied as indicated in Eq. (8). The coefficients of the four matrices are obtained imposing that (17)  , and are consistent with the predictions from the linearised model and the variations of the lift coefficients in the high-fidelity simulations. The higher frequency is larger than the frequencies range in the dynamic excitation in the present model reduction and thus it is not considered representative of the real system dynamic at high frequencies, but is required in the numerical model for the modulation of the system output.
The C l of the membrane is now estimated as a linear function of the states of the system: where C ou [ ] is the output matrix determined solving the system in the least square sense. The linear model obtained is used to produce an estimate of the numerical results from the high-fidelity system under purely harmonic excitation. The comparison of the C l for the high-fidelity system and the reduced-order linear model is presented in Fig. 14. It can be seen that the key dynamics is captured by the reduced-order description of the system. Since the POD is based only on the structural information, there is a high frequency content that is filtered out in the linear system and appears only in the flow behaviour. This is not exciting the structural part, and hence it is not captured in the ROM. Fig. 15 compares the spectral content of the structural response, Fig. 15(a), and the lift coefficient, Fig. 15(b), from high- Fig. 14. Comparison of the C l history of high-fidelity and reduced order models for the membrane under harmonic actuation, V 500 V 0 = .
fidelity and linear models for the actuation case at St The dominant forcing frequency is at St¼0.30 (twice the voltage frequency f v ). It drives the evolution of the system and it is the maximum peak in Fig. 15(a) and (b), followed by a second peak at the second resonance frequency of the aeroelastic system. The higher frequencies modulating the response of the lift coefficient in Fig. 14 are clearly visible in Fig. 15(b), but are filtered by the structural system. For the linear case, the frequency content of both structural response and lift coefficient is fundamentally the same, with a peak at the forcing frequency and the absence of the high-frequency content of the C l evolution. This explains the difference in the comparison of the high-fidelity and linear model predictions of Fig. 14. We conclude that this very low order description still provides a good estimate of the response and could facilitate the design of the control system for the integrally actuated wing.

Conclusions
This work has introduced a high-fidelity computational framework for the performance evaluation of integrally actuated membrane wings. Experimental investigations on such configurations are available in the literature, but to the authors' knowledge, this is the first time a high-fidelity model has been reported. A hyperelastic constitutive material model has demonstrated to be predictive of the performance of actuated inflated membranes. The fluid model is based on the direct numerical integration of the Navier-Stokes equations and its coupling with the structural solver has been verified against results from the literature. The integral actuation of the 2D wing showed that for the cases considered an increase of 15% of the aerodynamic performance of the wing can be obtained, although further investigation of the total energy balance needs to be considered to provide a fairer estimate of the potential gains. Finally, a methodology was introduced for the reduced order modelling of the actuated membrane based on proper orthogonal decomposition. The high-fidelity system was excited with input signals to characterise the dynamics of the actuated wing in the frequency range of interest, and the POD projection allowed the definition of a very low order model that can capture the main features of the system. The lift predictions of the actuated wing from the identified linear system have been compared with the high-fidelity model results, showing a good agreement and demonstrating the suitability of a reduction approach based only on structural information for the reduced-order modelling of this class of aeroelastic systems.