A morphing wing with cellular structure of non-uniform density

This paper proposes an optimization design method for the modular cellular structure of non-uniform density, which is filled into the morphing wing to generate variable torsion angle. By actively adjusting the distribution of the span-wise torsion angle, the lift distribution on the wing surface can be properly adjusted to avoid the problem of aeroelastic divergence or reduce the bending moment at the wing root. This ability is validated using CFD simulation. In the optimization framework proposed, the adaptive gradient algorithm is used to suppress the divergence of iteration. A finite element model with geometrical nonlinear effects is then proposed to correct the errors of the linear analysis and verify the effectiveness of the optimization method. This design is shown to be able to reduce the overall weight of the structure and achieve control of the macro mechanical performance of the wing. The work provides a general optimization design method for similar modular structures, allowing independent programmable adjustment of the parameters of each single structural cell.


Introduction
The concept of morphing has attracted widespread attention from present aircraft designers [1]. It is found that morphing wings can improve the overall flight efficiency by adapting to different flight conditions by adjusting their shapes [2], and can improve their performance within a certain airspeed range, reduce vibration [3], increase maximum lift, reduce drag [4,5], or enhance aircraft maneuverability [4,6,7] and aeroelastic performance [8,9]. In addition, due to the reduction of discontinuous rudder surfaces, morphing wings can further improve the stealth performance [2]. * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. For a long time, continuous increase in flight speed and altitude has increased the demand for structural rigidity because of aeroelastic problems. This means deformation of wings has become more difficult. This contradiction results in almost complete stalling of development in morphing wings. With continuous improvement of awareness of aircraft designers, materials and structural design technologies, engineers are now examining the structural design of aircraft from a completely different perspective. More diverse materials, processes and design methods allow more precise and effective configuration, and make the research of morphing wings revitalize again. Some researchers tried to produce morphing wings, blades or some adaptive structures driven by actuators of piezoelectric materials, and have developed some ingenious and effective deformable devices [5,10,11]. However, the weight of the basic structure still remains a formidable challenge, resulting in the limited application of heavy traditional actuators [12].
As a novel structural concept, the cellular structure may be a feasible way to solve the problem of weight in morphing wings. Using the stacked structure of modular cells, the macroscopic stiffness similar to that of elastic materials can be obtained at a mass density typical of an aerogel [13]. To improve the flexibility of wings, cellular structures are often used to make the elastic interlayers of morphing skins over the past few decades [14]. Nowadays, cellular structures have begun to be used to manufacture the main structures of morphing wings [13,15,16], rather than just a part of the skin. The essence of this concept is a three-dimensional ordered porous structure composed of periodic lattices and trusses, which are also called lattice structures sometimes. Because of its light weight, high specific strength and physical properties such as heat insulation, energy absorption and noise reduction, the cellular structure is considered to be one of the most promising lightweight, strong and tough materials. In addition to these excellent overall properties, decoupling and adjusting mechanical properties in single-material systems are also long-standing research goals of this kind of metamaterials [17]. Some advances in cellular structures have demonstrated that the overall performance of the structure can be controlled by tuning the parameters and topologies of single cells, and Zheng et al [18] have extended this result from the nanometer to the centimeter scale. Lattices of cellular structures can be made of polymers, metals or ceramics, exhibiting ultra-stiff mechanical properties at various scales [19].
However, when applying this extremely novel structure concept to the design of wings, designers rely heavily on their personal engineering experience to complete the transition from cells to a whole structure. Little has been done to develop a general design and optimization method of morphing wings. How to analyze and optimize the distribution of various cells still remains unclear. This makes the requirement for experience in cellular structure design very demanding. Due to the lack of relevant studies, previous design schemes used uniform density cells to fill the interior structure, and only few design parameters can be adjusted. Furthermore, considering the significant increase in the number of design variables input for non-uniform density cells, it will be difficult for general optimization algorithms to obtain a convergent solution [20].
In this work, we design a morphing wing with variable torsion angle using non-uniform density cellular structures to fill the wing. The structure is optimized by independently considering the relative density of each cell, so as to improve structural efficiency and reduce structural weight while ensuring adequate flexibility.
The primary difficulties are to obtain a universal design method for non-uniform density cellular structures, and to address the problem of convergence in optimization iteration. We establish an analysis model of cellular structures with the Euler-Bernoulli beam theory, and then test a traditional structural optimization algorithm and an adaptive gradient algorithm to compare their performance in solving this optimization problem. The adaptive gradient algorithm, AdaDelta, is finally taken. It is found to have strong robustness in solving large-scale inputs optimization problems, and has been applied to the training of complex neural network models [21].
Using the AdaDelta method, the convergence speed can be significantly improved, while some iterative divergence problems can be suppressed. Based on this universal optimization design method, we design a morphing wing with variable torsion angle as a case study to ground these arguments. A finite element model with geometrical nonlinear effects and Timoshenko beams is proposed to correct the calculation errors arising from the linear analysis process.
The work presented here contributes to the field of deformable components with cellular structures by demonstrating a universal method for analyzing and optimizing actively deformable elastic wings and cellular structures of variable densities. Using modular and reconfigurable cells, the adjustable macro performance and reversible assembly ability can be obtained.

Design and optimization of cellular structure
Before we fill the inside of the wing with the cellular structure, the major structural options must be determined, such as the geometrical topology, the characteristic length of the cell and the cross-section information of the poles in each cell.
Research on the macroscopic mechanical properties shows that regular octahedral lattice with high coordination number mainly transmit loads through the axial direction [15,19]. It is shown in figure 1, where the poles in compression and tension are respectively shown in yellow and blue. Considering the engineering calculation of axial force is more accurate, we can achieve a satisfactory analysis accuracy in the design phase with regular octahedral lattices. Therefore, this article takes regular octahedral lattice as the design object.
As an example, a straight wing with variable torsion angle is discussed to show a general optimization design method of morphing wings with cellular structure. This wing has a 300 mm × 580 mm × 36 mm size and NACA0012 airfoil, which could be used in a low speed small aircraft.
The remaining two features, lattice size and relative density, will be determined by optimization. To fill the cellular structure into the wing, the middle part is separated as the designing area, as shown in figure 2. The leading and trailing edges are not filled due to geometric constraints. In addition, the designing area is simplified into a flat block for better modeling and calculation efficiency.

Analytical model of lattice
To analyze the cellular structures, it is necessary to model and study its mechanical behavior during the loading process. Because of the typical feature of rigid frame structures, it is appropriate to use the Euler-Bernoulli beams to model the cellular structures. Euler-Bernoulli beams have been widely used in engineering. They are convenient to describe the bending and stretching deformation of beams, but the shear strain caused by bending in beams is ignored. This simplification is acceptable in the early conceptual design phase and could be easily corrected by using the Timoshenko beams in the finite element verification. In addition, some approximate performance models can help us to quickly estimate the macroscopic performance of the structural network. For the cellular structures, the main factors affecting its macroscopic mechanical properties are the geometric topology and spatial relative density of the cells [22]. The macroscopic equivalent elastic modulus of the structures can be expressed in the form of power exponent, as the following relationship: where E * and ρ * are the macroscopic elastic modulus and macroscopic equivalent density of the structure. E and ρ are the true elastic modulus and true density of the structural component material. The index α mainly depends on the geometric topology of the cell [23]. For a particular geometric topology, the proportional coefficient k depends on the direction of the applied load. The relationship also allows parameters to be calibrated through experimental or simulation results. A typical measure of the lattice shape is its coordination number, which is intuitively represented by the number of poles connected at each node. For high coordination number geometries such as regular octahedral lattice, since there is sufficient connectivity between the nodes, most of load is transmitted through axial tension and compression. The response of the element will mainly be in the axial direction of the rod, so the relative modulus changes almost linearly with the relative density. For low coordination number geometries such as Kelvin lattice [24], the beams transmit loads through bending, and the relative stiffness shows a quadratic relation with relative density [23]. The typical value of α can also be in other non-integer forms. For example, its value for the straight cubic lattice is about 3/2, coupling with two load transfer modes of axial tension and bending. However, Andrews indicates that the size of lattice has a critical effect on the performance of cellular solids [25]. In fact, the parameter that determines whether this approximate model is sufficiently valid is the resolution for a given cells network. The resolution refers to the number of cells in each direction of the structure. When the resolution is greater than 10, the performance of the cellular structure converges to the macroscopic performance given in equation (1). Otherwise, the above approximate model will lose some analysis accuracy. Obviously, our wing is difficult to satisfy this limitation, especially the number of cells in the thickness direction. This means that it is necessary to treat the cellular structure as a rigid frame rather than a continuous material. But the approximate model still has important guiding value for the estimation of performance, helping in pre-selection of materials. We choose the Ultem 2300 after a simple estimate, which is a kind of polyetherimide (PEI) with 30% glass fiber reinforcement, and is developed by General Electric Company.
The analytical model is based on the Euler-Bernoulli beam theory, which is widely used in engineering. Each pole in the lattice is regarded as a beam with 12 DOFs of two nodes, and the overall stiffness matrix is obtained by superposition of each sub-stiffness matrix. Then the displacement of the structure is solved by adding the boundary condition and the load. The solution method is linear, so it is impossible to consider the effects of structural nonlinearities, such as buckling loads and geometric nonlinearities common in large deformation problems. In order to avoid the influence of the singularity of the stiffness matrix during the solution, the row and column where the fixed node is located is crossed out as follows: The subscript S represents a fixed node, and F represents a free node. Thus, u S = 0. The displacement vector u F can be calculated in the block containing the free nodes.

Lattice size
To form the wing skeleton through stacking cellular structures, it is on the surface composed of lattice vertices that the skin will be installed. That means the support for this skin will be different from the traditional wing skin. However, the deformation of the wing skin must be considered in the design stage, because flexible skin is commonly used in morphing wings. Analyzing and controlling the deformation of the skin could help to control the adverse effects of excessive deformation such as bulges or wrinkles. As shown in figure 3, the supported span is equal to the characteristic length of a single cell. Obviously, the larger the lattice size is designed, the larger the skin deformation under the same aerodynamic load will be. But in engineering, a larger lattice size will be more conducive to processing and assembly. Therefore, the value of size should be selected to be a larger value that can limit the deformation to a reasonable range.
For a small piece of skin in the middle of four cells, because its displacements in the x and y directions are both small and not the investigation content, the mechanical model can be simplified to a simply supported thin plate. According to Kirchhoff plate and shell theory in elastic mechanics, σ x , σ y , τ xy are the main stresses, and effects caused by the other three stress components can be ignored. Therefore, it can be assumed that ϵ z = γ xz = γ xz = 0. Only the displacement w in the z direction needs to be considered:   Assuming that aerodynamic load on a small thin plate is uniformly distributed load P 0 , the deflection equation of this simply supported thin plate skin is: where D is bending stiffness of the thin plate: The form of the deflection equation (4) is a convergent infinite series. Since the series converges fast with the increase of the order, the deflection of the simply supported plate can be easily obtained by the previous items. Some extremely soft skins are usually used in morphing wings to ensure adequate deformation tolerance. We plan to use a light and soft polymer material for the skin of the morphing wing, which is called Ultem 1000. Its chemical name is polyetherimide (PEI), whose typical Young's modulus of elasticity is 3400 MPa. The skin is processed into films of thickness of 0.25 mm, and is then incised along the direction of the airflow, to reduce the additional stiffness in the direction of torsion. So the additional torsion stiffness is properly ignored, just as Cramer et al [13] and Jenett et al [15]. considered skins of their morphing wings.
The data of aerodynamic load is generated by CFD calculation, and only the static aerodynamic load distribution at each angle of attack (AOA) is taken. For the NACA0012 airfoil studied in this paper, the maximum lift coefficient AOA is about 8 • ∼ 10 • , so the range of AOA calculated by CFD is −4 • ∼ 12 • , and the gradient is 2 • . According to the aerodynamic load data and the deflection equation (4), the influence of AOA on the maximum deflection of skin can be calculated for each lattice size. That is shown in figure 4.
In the calculation result, it can be seen that when the lattice size is 36 mm or less, the maximum skin deflection does not change much with the angle of attack. In this condition, the skin bulge has less impact on the aerodynamics. Considering the convenience of processing and assembly, the lattice size is decided to be 36 mm finally.

Relative density distribution
It is shown by Jenett that morphing wings could achieve better performance by designing distribution of different performance cells [15]. Therefore, relative density distribution of cells could be subtly considered by optimization to further improve structure performance. Which should be noted is that the different cells considered in this article have different cross-sectional areas only. Considering the actual production, every cell will have exactly the same three-dimensional topology, and the properties of the beams contained in each independent cell are the same too.
Topology optimization is an optimization method that considers the best distribution of materials in space, and is widely used in the field of structural design. The relative density distribution of cellular structure is also a study of the best distribution of materials in a broad sense, so some ideas of topology optimization can be used to solve it. The difference from the general topology optimization problem is that the intermediate density value has no physical meaning in the general topology optimization problem. There are only two states of 'being' and 'nothing' in real space for materials distribution. So we must to make the optimization result converge to the '1/0' state. But in the cellular structure density distribution optimization, the relative density can have intermediate state, so there is no need to make the optimization result converge to the '1/0' state.
The problem of this optimization design can be simply described as: using the regular octahedron lattices model as the optimization object, the optimization parameter is the crosssectional area of each regular octahedron cell, and the optimization goal is to make the torsion deformation under the action of the driving torque closest to the expected torsion angle. Considering the lightweight requirements of the structure, a penalty item about structure weight can be added. The optimization problem studied can be recorded as the following mathematical model: subject to: where u L and u T are the displacements in the z direction of the leading and trailing edge nodes in the wing tip, L chord is the chord length of the wing, θ W is the expected torsion angle, and x is the sequence of the cross-sectional area of the beams in each cell. The two constraints represent the displacementsforce equilibrium equation of the structure under driving load and the minimum cross-sectional area limit of the beams to avoid the singularity of the stiffness matrix. The objective function represents the square error between the wing tip torsion angle and the expected torsion angle, and the weight of the structure ∑ x is punished by the penalty coefficient p. The structural load f is composed of two parts: the aerodynamic load and the driving load. The aerodynamic load is calculated with the target shape, and is transmitted to the structure through only one time of information interaction. The driving load is set as 0.4 N · m to guide the wing to be washed in. Due to symmetry, when the driving torque is set to the opposite direction, the deformation direction of the structure is also the opposite. The specific value of the driving torque is mainly affected by the structural stiffness, and is limited by the output force of the small servo motor. Meanwhile, the required deformation rate and the typical aerodynamic load should also be taken into account. Therefore, its initial value is decided with the engineering experience of designers, and its exact value could be amended by extracting the supporting reaction of the nodes.
For such optimization problems with lots of design variables, the global optimal solution cannot be obtained analytically. To find the numerical solution of the optimization problem, the common methods include gradient-based solvers and probability-based solvers. The gradient-based solvers use the gradient information to iterate repeatedly in the inverse gradient direction to approach the optimal solution, such as batch gradient descent method (BGD) [26,27] and moving asymptote method (MMA) [28]. The probability-based solvers use the probability principle to approach the optimal solution from a random direction on the basis of a large number of iterations, such as genetic algorithm (GA) [29] and simulated annealing algorithm (SA) [30]. Generally, a gradient-based solver requires less calculation, but it is not until the solver could work that the gradient information is provided. Considering the solution efficiency, the gradient-based solver will be taken. According to the finite element equations (2) and crossing out the block where the fixed node is located, the free node displacement can be obtained: The gradient analysis of the objective function is carried out under the framework of the variational principle. Due to the form of the equilibrium equation, the derivative of displacement u with respect to the design variable x is not easy to solve. It is necessary to use the adjoint method to derive and determine the gradient of the objective function. Construct an augmented Lagrangian function: where λ is the Lagrangian multiplier vector related to the equilibrium equation, which is called the adjoint vector. Because of the equilibrium equation,c = c is always valid. According to chain rule, the total derivative of formula (8) is: Considering that the original equation remains valid for any λ, in order to avoid solving the derivative of the displacement with respect to the design variable, a specific λ can be selected to satisfy the following equation: Thus, The gradient of the original objective function is: Since that the design variable x is not explicitly included in the objective function c, ∂c ∂xj = 0. In the boundary conditions of this optimization problem, the external load f is applied to the determined node with a constant value, and is independent of the design variables, so dfF dxj = 0. The influence of the beams cross-sectional area x j on the equilibrium equation only acts on the relative position of the cell j in the stiffness matrix K, so dKFF dxj = (K j ) FF is the sub-stiffness matrix assembled by the expanded stiffness matrix of all the beams in cell j. Among them, only the remaining blocks are working after removing the rows and columns where the fixed nodes are located.
In this optimization problem, due to the lots of the design variables and the striking difference between orders of magnitude of the gradients of the objective function with respect to each design variable, the directly using of BGD method to solve the optimization problem is extremely inefficient and easy to cause diverging. This is always called sparse parameter problem, where we cannot easily address the problem by normalization of design variables. In fact, since cell densities are expressed as dimensionless relative densities, the design variables in this problem have been normalized already. Some studies in the field of machine learning have pointed out that traditional gradient optimizers have poor versatility when dealing with complex optimization problems. The main reason is that the parameter update speed during each optimization iterations cannot be automatically adjusted according to the iteration information, which means the result depends on the selection of initial parameters a lot. The adaptive gradient algorithm has been improved to address this problem [21,31]. The latest algorithm of them is called AdaDelta, which can significantly improve the robustness of the gradient optimizer and has been used to train large-scale neural networks. Therefore, we have tried to solve the structural optimization problem by using AdaDelta method. The core of this algorithm is the update rule of design variables: where η is the specified initial update rate, and ε is a minuscule value added to prevent denominator from being zero. E[g 2 ] t is the cumulative value of the square of gradient in a window period before the current iteration step. However, because it is relatively inefficient to reserve separate memory to store the cumulative value and update the window period repeatedly in each iteration step, the numerical attenuation is used to replace this cumulative effect: The above formula is equivalent to that the accumulation of square of the historical gradient is multiplied by an decay rate ρ, and then (1 − ρ) is used as the weight coefficient of square of the current gradient. A calculation example shows that the AdaDelta method can significantly improve the convergence rate of this optimization problem.

Optimization settings
A morphing wing with the NACA0012 airfoil is then considered in this article. The expected torsion angle θ W is 10 • , and the weight penalty factor p in the objective function (6) is 5 × 10 −10 . The initial update rate of the relative density of all cells is specified as η = 0.1. To make use of the lightweight feature of cellular structures, a kind of super engineering plastics with the elastic model of 9.3 GPa, Ultem 2300, is selected for the wing skeleton. Its ultra-light feature of 1510 kg m −3 density contributes to that this wing skeleton weighs only 155.3 g after being optimized. Figure 5 shows the iterative calculations of the AdaDelta method, where the value of objective function and update rate are respectively represented by the red solid lines and the blue dashed lines. Update rate means the maximum change in the relative density of all cells in the cell structure during each iteration, as given by It can be observed that the objective function converges rapidly as the iteration proceeds. Even though the initial update rate is specified as 0.1, the algorithm adaptively improves the update rate in the first steps of the iterations according to the parameters, the gradient information and the iteration effect. More importantly, when the iteration result is close to the convergence solution, the algorithm automatically reduces the parameter update rate to avoid missing the optimal solution. That will be good news for the engineers, including designers who are not experts in structural optimization, because it means the algorithm can take advantage of its strong adaptability to help them quickly determine an effective update rate for parameters, rather than picking it among a large number of initial update rates that may finally lead to divergence. Because of these advantages, the Adedelta algorithm finally achieves the convergence under the calculation cost of less than ten iteration steps. Surprisingly, with the exception of the initial update rate, the different initial parameter states x t=0 over a relatively large range do not seem to lead to divergence in the results. As a contrast, we also give the iterative calculations of the optimization criterion (OC) method established through gradient descent in figure 5, whose data points are marked as circles. It should be specially explained that because it is difficult to simply adjust the update speed to make the result converge, we also tried to use some techniques to help with convergence. The result curve shown in the figure is generated by the OC method improved by Sigmund [32]. He developed a technique called the mesh-independency filter to improve the convergence of the OC algorithm by convolutional pretreatment of the gradient information. At the same time, an iteration scheme with Lagrange multiplier was adopted to try to reduce the update rate near the optimal solution. However, after trying a series of parameter combinations, we still cannot get a convergent result through these two improvements. As shown in figure 5, the OC algorithm does not use the intermediate information generated during the iteration to improve the update rate, nor does it learn the characteristics of the parameters themselves. This may be the mathematics essence of its difficulty in convergence when faced with complex optimization problems or sparse parameter problems. Sigmund also points out that the convergence of these assistive techniques cannot be proved mathematically, although some improvements have been observed in the field of structural topology optimization. According to the comparison, the AdaDelta algorithm with adaptive update speed is finally adopted to solve the optimization results.

Parameters analysis
In the optimizer settings, there are three main parameters that will determine the optimization effect and convergence speed: initial cell densities x 0 , initial update rate η and decay rate of historical gradient information ρ. We have explored and tested several parameter schemes and obtained some simple rules, which may help researchers to determine the appropriate parameters in their specific optimization problems.
The selection of initial point x 0 does not determine whether the algorithm converges to the correct optimal solution. In the optimization problem we proposed, according to the continuity of the design space and deformation results, as well as the excellent robustness of the adaptive algorithm, when x 0 change within the range of 0.1 to 1.0, the optimization will give the same solution. This conclusion is also applicable to the choice of η and ρ, but the difference is that the result may cannot to ensure convergence under some extreme values.
The selection of all three parameters above affects the convergence rate of the optimizer, where the initial update rate η and the decay rate ρ are the main influencing factors. The mathematical mechanism of this phenomenon is the updating rule of the Adadelta algorithm. If keeping these two parameters at the certain values, the algorithm automatically accelerates the updating of designing variables at the early stage of iteration, so that the effect of the initial point x 0 changing within a small range will be quickly diluted by the adaptive step length.
A large initial update rate means a more aggressive search strategy, which can speed up the process of the optimizer to find the optimal solution, but this may usually cause the optimizer to repeatedly oscillate near the convergence point in the later iterations. Figure 6 shows this phenomenon when the initial update rate η is set as 0.2. Compared with the case of η = 0.1, the optimizer reached the convergence point at a faster speed under this parameter. However, the optimizer missed the convergence point for several times due to the excessive updating rate, and finally consumed more than twice the number of iteration steps to reach the convergence. To strike a balance between search efficiency and convergence speed, we recommend using approximately one-tenth of the initial relative density as the initial update rate.
The function of the decay rate ρ is to adjust the influence ratio of the historical gradient information and the current gradient information on the current update rate. This is also the core reason why the Adadelta algorithm can make effective use of the iterative data itself to guide the next iteration and to obtain higher robustness and convergence speed. When the value of the decay rate is closed to 1, it means that the algorithm is more inclined to guide the updating rate of the iteration through the historical gradient information. On the contrary, when the value of the decay rate is 0, the algorithm degenerates into a complete BGD method, which iterates at a fixed updating rate. The recommended value of ρ will depend strongly on the features of the optimization problem itself. This is because the descent velocity of gradient of each parameter is significantly different in different optimized objects, which makes it difficult to select the matching ρ in a general method. Through a simple experiment, we found that about 0.8 ∼ 0.95 is a relatively good choice in our wing object. As a mathematically unproven recommendation, we prefer to use a decay rate value ρ that is larger but still less than 1, based on the apparent phenomena observed from the simple experiment. This suggestion is supported by the work of the original authors [21], where a decay rate of 0.95 is applied. If it is observed that the optimizer cannot achieve convergence, we recommend using the dichotomy to find a more suitable one. Figure 7 shows a schematic diagram of the optimization result obtained with the above variable density topology optimization method, where the left and top sides represent wing root and leading edge respectively. The gray information of the image indicates the relative size of the cross-sectional area of all the beams in cells, which is relative density thereof. The position of each cell is described by its X and Y location coordinates. For example, the cell in the bottom right corner of the image is located with (6,16). When the optimization result is translated into the actual wing structure, the shape of cells has been adjusted due to the fact that the leading and trailing edges of the wing are too thin to meet the dimension requirements of a complete regular octahedral cell. And the cells with low densities, such as that whose relative densities are less than 0.1, are ignored in the actual structure modeling. Because they provide little stiffness contribution and it is difficult to guarantee their machining accuracy. Figure 8 provides the actual wing structure translated from optimization result, in which the cells are dyed into different colors to show their diverse equivalent densities. More details of the beams in each cell are given in table 1. The x-axis and y-axis are along airflow and wing span, and z-axis is perpendicular to these two directions. It can be seen that the crosssectional area of beams and equivalent relative density of cells are maximized at the thickest part of the airfoil, where spars will be located in tranditional wings. This is consistent with our intuitive prediction.

Finite element analysis and verification
In order to verify the rationality of the optimization results, a finite element model (FEM) was established. Figure 9 shows the numerical simulation result of the deformation in z-axis under the driving load given in design phase, where the deformed structure is colored by the value of the displacement in z-axis. The highest and lowest positions are shown in red and blue respectively in the displacement contour.
The accuracy of the analysis depends largely on the quality of FEM. Figure 10 shows some details of the meshing of FEM. To ensure the sufficient analysis accuracy, it is checked carefully that each beam has been divided into at least four beam elements and each beam element in the FEA model has an appropriate slenderness ratio.When the grids density is increased to about 1.4 times of that of the original grids, no obvious change could be observed in the simulation results of structure deformation.  Furthermore, beam elements used in the FEM are constructed based on the theory of Timoshenko beam, and the 'Large Deflection' switch in Solver Control is turned on. Timoshenko beams help to correct errors caused by neglecting the shear deformation in the Euler-Bernoulli beams. Following the Mechanical User's Guide, Large Deflection is a suggested setting for hyperelastic material models. In fact, when the torsion angle of the wingtip reaches 10 • , it is a reasonable modification to the original linear analysis model to consider the nonlinearity caused by deformation. The geometric nonlinear analysis with large deflection enabled allow designers to expediently improve their analysis accuracy for various morphing wings, especially those made of hyperelastic materials or engineering plastics.
By enable the Large Deflection function, the calculation program will divide the real loads into several load steps and gradually apply them to the structure, instead of loading them all at once [33]. At the end of each load step, the balance of forces is ensured by iteration and the stiffness matrix of the structure is reconstructed. This process can be automatically controlled by the solver, but it is usually impossible to ensure a convergent result. If observing that the results cannot be converged in the calculation, it can be tried to manually divide the load into more steps, or enable the automatic load steps determined by the dichotomy. That means it may be difficult for designers to analyze lots of different structural schemes through a standardized and automated process, and it will take a disappointing amount of time, making the optimization approach less feasible. It is an important reason why the nonlinear structure analysis has not been embedded in optimization yet.
Compared with the case of using the geometric nonlinear analysis including the Large Deflection function, the conventional linear structure analysis saves 60% of the time and 5.4% of the computer memory usage for each time calculation, but the z direction deflection of the structure produces an error of up to 18.6%, resulting in a 17.7% error in the torsion angle.
Comparison of the optimal structure and the original structure shows an enhanced deformability. Under the same driving load M = 0.4 N · m, the optimal structure give a greater torsion angle of 11.9 • . The torsion angle increases by 54.5% compared with that of the original structure of 7.7 • . Since the size of the stringers is not changed, the span-wise bending stiffness of the wing structure is reduced slightly after optimization. This means the optimal structure achieves a significantly reduced torsional stiffness while keeping a certain span-wise bending stiffness, leading to an effective stiffness allocation.
According to the results of the nonlinear finite element analysis ( figure 9), under the original design driving load, the torsion angle generated by the wing structure is 11.9 • . That means there is a 19% positive error compared with the expected torsion angle of 10 • set by the optimization target. But in fact, a positive error meets our estimates. We use an objective function of the composite optimization problem involving the deformation target and the weight target, rather than a simple deformation objective function. That means our design approach provides a compromise optimization result considering the deformation and the weight. To achieve a lower weight of the wing structure, the optimizer further reduces the amount of the material used, which reduces the stiffness of the wing skeleton. So the driven deformation cannot perfectly match the deformation target proposed in the design phase, but slightly larger than the target. If keeping this adjustment, engineers can still achieve the deformation target without increasing the output power and torque of the driving system. So it still measures up to the engineering requirement. Certainly, this effect can also be adjusted easily by selecting different values of the penalty factor p.
Secondly, the Large Deflection function in structural analysis is enabled. Under the deformation of a large scale, the influence of geometric nonlinearity makes the linear structure analysis in the optimization method unable to fully reflect the real deformation of the structure. The FEM is modeled with Timoshenko beams, so it has corrected the analytical model of Euler-Bernoulli beams in structural design. These improvements will also result in slightly different deformation results.
However, due to the efficiency and feasibility of optimization iteration, we still have not used the structural model completely considering geometrical nonlinear effects and the shear deformation in Timoshenko beams in the optimization method, but only use high fidelity FEM for precision correction before detailed design. This technique, where analysis methods with different precision are both used at different stages of design, is called variable-complexity modeling. Variable-complexity modeling has been widely used in the field of structural design, and is considered to be a good choice to strike a balance between computational efficiency and analytical accuracy.
According to the modeling and analysis process, there is another minor source of this error. To match the airfoil shape of the wing, the size in the z direction of each cell in the actual structure has been modified, resulting in a difference between the shape of actual structure and the optimization model. The design presented here shows an extreme situation of a smaller-sized wing. Since there is only one layer of units inside the structure, the shape adjustment of the leading and trailing edges has a great influence on the overall structure. When this method is applied to a structure containing multilayer cells, this influence could be significantly reduced.
To show the deformation of the morphing wing under aerodynamic loads more accurately, static aeroelastic effects of the wing are considered and verified additionally. The looselycoupled method for static aeroelasticity is used in the calculation. In this method, the aerodynamic loads and structural deformation will be iterated and calculated repeatedly until the results converge. Because of the different characteristics between structural mesh and aerodynamic mesh, the transmission of load and displacement information between meshes needs to be realized by interpolation method. For more details on aerodynamic calculations, please refer to section 3.4. In this paper, the radial basis function (RBF) interpolation with the kernel function of Wendland C4 is chosen, which has shown good performance in the aeroelastic problems of structures [34]. It should be pointed out that the RBF interpolation with the kernel function of Wendland C4 can maintain the equivalence of force and torque, and maintain the equivalence of virtual work at the same time. It is the theoretical premise of the reasonability of the interpolation. The result of the deformation of the morphing wing considering static aeroelastic effects is given in figure 11. It can be seen that the coupling effect of aerodynamic loads and structural elasticity makes the deformation of the structure slightly different from that without coupling, but it does not affect the concept and driving ability of the structure on the whole. Considering the balance of computational efficiency and analytical accuracy, this coupling technique is only used for final verification at present.

Driving system
To apply the driving load, a simple driving system is designed and used as an additional component to form an intelligent morphing wing system by added into the wing skeleton. In the driving system, a rotating servo motor is used as the source of driving load, and is installed inside the fuselage. A torsionshaft is directly connected to the servo motor, and is connected to the wing tip through a base with a sleeve. Figure 12 shows the intelligent morphing wing designing scheme with the driving system, where the additional driving system is highlighted in yellow. This driving system works by installing additional components, without changing the skeleton part that has been optimized already. The connection between the various parts of the actuator is realized by some small screws, which are located in the positions of each hole in the figure 13. (The screws themselves are not shown.) Arrows point to the location of the screw connection. And arrows with dotted lines indicate holes that cannot be directly shown due to occlusion.

Aerodynamic performance
The morphing wing could improve its aerodynamic performance of the wing both in small-and large-AOA flight by actively adjusting distribution of the span-wise torsion angle. This ability is demonstrated using computational fluid dynamic (CFD) simulation. The flow around the wing at the airflow velocity of 30 m s −1 is built and calculated in Fluent for the original configuration and the morphing wing. Two improvement z -displacement (m) Figure 11. Deformation of the wing considering static aeroelastic effects. The CFD model contains 1.18 million nodes and 1.92 million cells, and the SST k-ω model is chosen as turbulence model to simulate the flow around the wing in the single side. In particular, the grids near the wing surface is refined. The boundary layer area contains 40 layers of grids, in which the height of the first layer of grids on the wing surface is 0.01 mm. This mesh has been tested for grid independence to ensure that further mesh refinement will not change the result of lift significantly. The aerodynamic grids are shown in figure 14. The boundary conditions are set as figure 15, where c is the chord length of the wing, and b is the span length of the wing. The wall y + (Yplus) of whole field is kept in the range of 0.2 to 2, which is suitable for the SST k-ω model. This means that the quality of mesh is acceptable.
As an available method of lift enhancement, the morphing wing can use a positive tip torsion angle (i.e. the wing is washed in) to increase the effective AOA near the tip when the overall AOA is relatively small. Figures 16(a)     closed to 0, that means the flow does not separate from the surface. This may because the maximum cross-sectional AOA in the model is 12 • , which is still confined by the stalling AOA. Compared with the traditional lift enhancement device, the morphing wing has a larger deformation area, so it can achieve stronger lift enhancement effect. Obviously, when the overall AOA of the wing is large, the lift distribution near the wingtip can also be improved by a negative tip torsion angle (washed out), and the bending moment at wing root can be reduced. Figures 16(c) and (d) show the pathlines coloured by airflow velocity. Similarly, the flow does not separate from the surface. By the wing tip washing out, the lift generated by the wing decreases from 58.52 N to 44.58 N. However, the lift is mainly reduced near the wing tip, resulting in a significant reduction of the bending moment at the root, from 15.49 N · m to 11.03 N · m. By reducing the root bending moment, this morphing wing is at less risk of divergence and damage. Compared with the traditional wings, the morphing wing achieve actively adjustment of lift distribution, enhancing its aeroelastic performance.

Improvements and future applications
The error between the driving deformation and the expected deformation can be further improved. In the optimization model, the penalty factor p can be adjusted to improve the deformation accuracy thereof, but this will weaken the weight reduction effect of the design. In addition, using the cellular structures on larger-sized wings will also help improve this accuracy. Because a structure containing few cells shows a performance close to that of a discrete rigid frame rather than a monolithic material. As the size of the wing and the number of cell layers increase, the thickness of the structure can be adjusted by increasing or decreasing the number of cell layers, or the thickness adjustment can be applied to each cell layer on average, thus reducing the influence of the structure thickness adjustment on the size accuracy of each single cell.
At the same time, the influence of geometric nonlinearity can also be considered in the optimization design to improve the calculation accuracy of the structure solver, and the deformation optimization ability of the optimizer can be improved immediately. However, this choice could only be applicable to some simple structural concepts, and the feasibility of the optimization method will soon become ineffective due to the expensive computational cost as the complexity of the structure increases. The trade-off between computational accuracy and efficiency is a key issue that designers must consider.
By adding an item about deformation under different loads to the objective function, in this method we can also achieve the structural optimization considering multiple loading cases easily through the following form: where the subscript i represents different loading cases, and p m is still a penalty factor of structure weight. If the convenience of the production process is further considered, sometimes we want the cells to have the same relative densities and only optimize whether cells need to be filled in each spatial position. It is also easy to achieve this goal through the optimization method above. This requires a minor change to the optimization model only. Following the idea of the solid isotropic material with penalization model (SIMP) [35], we can penalize the intermediate relative densities of cells, such as the part at 0.2 ∼ 0.8, to get a calculation result that tends to be close to 0 or 1. Specifically, instead of using the actual structural stiffness in the calculation, we replace it with a form with a relative density penalty factor. For example, the following penalty model is adopted to replace the original elastic modulus of the beam material: where E * (x) represents the elastic modulus after punished of material of the beam with the relative density x. x 3 is used as a penalty factor to punish the intermediate relative densities between 0 and 1. As the stiffness provided by the cells of intermediate densities is reduced but their weight are still considered by the weight factor, the optimizer will be more inclined to use the relative densities close to two sides. Due to this change, the analytical formula of the gradient information during optimization should also be adjusted. In the original gradient analytical equation (12), the sub-stiffness matrix (K j ) FF will be multiplied by 3x 2 j , since the elastic modulus E is multiplied by x j 3 here, thus Figure 17 shows the result given by the optimizer, in which the relative densities punishment is considered. The result is consistent with our analysis. The optimizer gives the optimization results of the cellular structure of relative densities close to 0 or 1. In addition, the decay rate ρ is set as 0.95 to help to get a fast convergence. As shown, given the discrete constraint that we can only use the cells with a specified density, this optimization method still provides an effective guide to decisionmaking for engineers. Engineers can assume different target deformations and allowable cell densities and run optimization programs to observe the contribution of cells in various spatial positions to the structure, and then decide whether they need to remove some cells to obtain a lighter structure weight.

Conclusions
We have described a design method and an analysis model for a morphing wing with the cellular structures of non-uniform density, where the structure weight is taken into consideration, and each cell in the structure is carefully considered by the variable density optimization method. We achieve a morphing wing with a skeleton weight of only 155.3 g. The weight of the cells other than stringers and ribs can be reduced by 36.0% compared with the baseline of 242.6 g. The optimal structure shows a torsional deformability increased by 54.5%. By adaptively adjusting the span-wise torsion angle, the aerodynamic or aeroelastic performance of the wing can be further improved, as shown in CFD simulation. Rapid reversible assembly is achieved through stacking cells, helping to save the cost. After optimization, we use a finite element model built to validate the effectiveness of the optimization method and analyze the optimized structure. In the FEA, the influence of geometric nonlinearity in large deformation is considered to correct the errors in linear analysis. The AdaDelta method is used to calculate the optimization results in this work. The method shows strong robustness in solving the optimization problems of cellular structures with large-scale design variables and sparse parameters, and can solve divergence of the optimization results in the traditional gradient algorithm. Some suggestions about parameters selection in the AdaDelta algorithm are given based on a test of parameters.
A FEM is built to validate the effectiveness of the optimization method and analyze the optimized structure. The geometric nonlinear effect and the shear deformation in beams are both considered in the FEM to correct the analytical errors. This principle of variable-complexity modeling allows our optimization method to strike a balance between solution efficiency and analysis accuracy.
Some improvements of this optimization design method are discussed. We give an expression for the optimization model considering multiple loading cases. In addition, we also give an optimization model under the constraint of discrete value of cell density and its calculation case. We believe these improvements could help designers tackle more complex problems with our design method.
The cellular structures with non-uniform density can significantly expand the available design space, and enable further enhancement of performance on the basis of the excellent mechanical properties of cellular structures with uniform density. This work provides engineers a solution of cellular structures with good usability and reliability, and contributes to the field of ultralight deformable components.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).