Next Article in Journal
Sequence Spaces and Spectrum of q-Difference Operator of Second Order
Next Article in Special Issue
Optimization Design of Unbonded Areas Layout in Titanium Alloy Laminates for Fatigue Performance
Previous Article in Journal
Key Component Capture and Safety Intelligent Analysis of Beam String Structure Based on Digital Twins
Previous Article in Special Issue
On the Indispensability of Isogeometric Analysis in Topology Optimization for Smooth or Binary Designs
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Aeroelastic Topology Optimization of Wing Structure Based on Moving Boundary Meshfree Method

1
Institute of Unmanned System, Beihang University, Beijing 100191, China
2
School of Aeronautic Science and Engineering, Beihang University, Beijing 100191, China
3
Zhejiang Key Laboratory of General Aviation Operation Technology, General Aviation Research Institute of Zhejiang Jiande, Jiande 311612, China
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(6), 1154; https://doi.org/10.3390/sym14061154
Submission received: 17 May 2022 / Revised: 30 May 2022 / Accepted: 1 June 2022 / Published: 3 June 2022

Abstract

:
The increasing structural flexibility of large aircraft leads to significant aeroelastic effects. More efficient topology optimization techniques are required for the design to further take advantage of aeroelasticity and obtain lightweight structures. This paper proposes a moving boundary meshfree topology optimization that combines the Galerkin method of weighted residuals and non-uniform rational B-splines (NURBS). The solution domain is described by the control points of NURBS and its property is calculated adaptively with an integration subtraction technique. The minimal compliance is searched for using the globally convergent method of moving asymptotes (GCMMA) by designing the locations of control points as subject to volume and flux constraints. The method is first applied to a typical two-dimensional design example with symmetric boundary conditions. The results show that the shape constraints can be conveniently applied, and smoother boundaries are obtained with fewer parameters. Then, a three-dimensional wing structure with asymmetric boundary conditions is optimized. A three-dimensional flight load that combines the high-order-panel and meshfree methods is employed to calculate the elastic loads and update asymmetric external loads during the optimization process. The designed wing satisfies engineering requirements and the presented method can solve the practical topology optimization problems of three-dimensional structures.

1. Introduction

Structural optimization is an efficient approach to seeking the optimal structure distribution in engineering design. Structural optimization in engineering can usually be divided into three categories: size optimization, shape optimization, and topology optimization. The topology optimization of aircraft is generally applied in the conceptual stage and it has more design parameters and a larger design space than size optimization and shape optimization [1]. Therefore, topology optimization plays an essential role across the whole design and the optimization results are the basis of subsequent designs.
Topology optimization was developed first in aerospace engineering due to the strict requirements of high-performance and lightweight structures. Researchers and airline industries have invested much effort in structural optimization. The topology optimization of the full-scale wing structure of a Boeing 777 adopted super-large-scale parallel computing and obtained a wing design scheme similar to bionic structures [2], and the weight of a single wing was reduced by 2% to 5%. Groen [3] proposed a novel method to obtain a near-optimal frame structure based on the solution of a homogenization-based topology optimization model. Ghasemi [4] proposed a high-performance density-based topology optimization tool for laminar flows with a focus on 2D and 3D aerodynamic problems via the OpenFOAM software. Munk [5] developed a design methodology that reduces the maximum stress in a component for a given applied load and weight constraint, and the sustainability of the aircraft was increased. Palacios [6] investigated the use of density-based topology optimization for the aeroelastic design of very flexible wings using a nonlinear aeroelasticity solver and a gradient-based approach. The structural weight was reduced while maintaining the lift. Li [7] presented a multiscale optimization method to solve size distribution problems in conformal lattice materials, and the design resulted in mechanical experiments on specimens fabricated by 3D printing. López [8] compared the structural designs derived from deterministic topology optimization and reliability-based topology optimization and proved that including uncertain data in the topology optimization can help to reduce the weight of the component.
The current topology optimization methods mainly include the level set method (LSM) [9], evolutionary structural optimization (ESO) [10], and solid isotropic material with penalization (SIMP) [11]. All the above methods are based on the finite element method (FEM) and SIMP is widely employed in the topology modules of commercial FEM software systems. The increasing demands of the low flexibility structures of modern aircraft lead to significant structural deformation and non-negligible aeroelastic effects [12], which are a challenge for traditional topology optimization. For FEM-based optimization methods, if the boundary of the solution domain is changed, matched meshes are needed and the global stiffness matrix should also be partially corrected. The above process is complicated and, once aeroelasticity is considered, the topology optimization becomes tedious and time-consuming.
The analysis domain of meshfree methods is represented only by arbitrarily discrete field nodes. It is possible to overcome the limitation of tensor product mesh in complex topology problems using isogeometric analysis [13]. Belytschko [14] published a comprehensive review article detailing meshfree methods in 1996, which is considered the beginning of meshfree methods as an independent branch of mechanics. Although they started late, meshfree methods have developed various forms with differing discretization and approximation. The essential software framework developed gradually [15]. So far, meshfree methods have been applied widely in astrophysics, microparticles, hydrodynamics, solid mechanics, fracture mechanics, explosion mechanics, thermodynamics, biomechanics, and electromagnetic mechanics [16,17]. The proposal of the element-free Galerkin method (EFG) [18] is important progress, and it can be combined well with the approximation techniques such as the radial point interpolation method (RPIM) [19], Kriging method [20], reproducing kernel particle method (RKPM) [21], etc. The EFG greatly extends the application range and has become a very competitive meshfree method.
This paper presents an efficient technique named moving boundary meshfree topology optimization that can solve complex topology optimization problems in aircraft wing structures. The Galerkin method of weighted residuals and RPIM are employed to discretize the system equations. The control points of non-uniform rational B-splines (NURBS) are employed to exactly describe the analysis domain. The domain property is calculated adaptively through integration subtraction and the globally convergent method of moving asymptotes (GCMMA) is used to search for the optimal topology structure subject to the volume and flux constraints. The two-dimensional beam proposed by Messerschmitt-Bolkow-Blohm (the MBB beam) is a typical example of topology optimization and it is used as the first design object for the proposed method. National Aeronautics and Space Administration (NASA) provide many common research models (CRMs), and a three-dimensional wing structure with a high aspect ratio is selected as the second design object. The three-dimensional flight load is used to calculate the elastic load and update the external loads of the topology optimization. The proposed optimization framework is expected to be a reference for the structural design and aeroelastic optimization of large aircraft.

2. Moving Boundary Meshfree Method

2.1. Meshfree RPIM Method

For a linear elastic problem in domain Ω of the boundary Γ , the governing equation is described as [22]:
L T σ + b = 0 in Ω n σ = t ¯ in Γ t u = u ¯ in Γ u
where L T is the differential operator; σ is the stress vector; and b is the external force vector. t ¯ and u ¯ are the applied stress and displacement vectors on the boundary, respectively. Using the Galerkin method of weighted residuals, the general weak form of Equation (1) is given by
Ω ( L δ u ) T D ( L u ) d Ω Ω δ u T b d Ω Γ t δ u T t ¯ d Γ = 0
where D is the elastic matrix. Unlike FEM meshes, the meshfree method employs scattered field nodes to discretize the complex solution domain Ω , as shown in Figure 1.
The information of an unknown node is interpolated by the field nodes of the support domain around the unknown node. RPIM is an efficient and stable interpolation method [23] and is used to approximate the field function. The value of the field nodes in the support domain is written as:
U ˜ s = U s 0 = R 0 P m P m 0 a b = G a 0
where R 0 is the radial basis matrix related to the distance of the Euclidean distance of the field points. P m is the polynomial matrix, and a and b are the undetermined coefficient matrices. The unknown node values in the support domain are approximated by:
u ( x ) = R ( x ) a + P ( x ) b = R ( x ) P ( x ) G 1 U ˜ s = Φ ( x ) U ˜ s
where Φ ( x ) is the shape function of RPIM that represents the influence of each field node in the support domain to the unknown node.
Combined with RPIM, the solution domain is approximated by finite field nodes and the global integration can be further discretized to the Gaussian integral of several background meshes. The integral of the governing equation is usually transformed into the parameter domain for calculation, which is written as:
Ω f ( x ) d Ω = k = 1 n c Ω k f ( x ) d Ω = k = 1 n c i = 1 n k ω i f ( x i ) J i k
where Ω k is the k-th background mesh; n c is number of background meshes; n k is the number of Gaussian integral points of the k-th background mesh; ω i is the weighted factor of the i-th Gaussian integral point; and J i k is the Jacobian matrix of the k-th background mesh at the i-th point. The background meshes are divided conveniently by homotopic mapping after deploying field nodes. Therefore, Equation (2) can be written in matrix form:
K U = F
where K is the stiffness matrix, and U and F are the vectors of displacement and applied load, respectively. The deformation, stress, strain, and other responses can be solved by applying the boundary condition.

2.2. NURBS

NURBS curves can efficiently describe arbitrary free-form shapes and conic sections, which makes NURBS the basis for all standard geometric exchange formats [24,25]. NURBS is a linear combination of B-spline basic functions and the B-spline curve C ( u ) is defined as:
C ( u ) = i = 0 n N i , p ( u ) P i a u b
where u is the parameter of the knot vector between a and b and P i are the control points to determine the shape of the curve. N i , p ( u ) is the i-th B-spline basic function of degree p , which has the important properties of nonnegativity, unit partition, local support, and differentiability.
The NURBS surface is defined using the tensor-production method and the generalized forms of the NURBS curve and surface with normalized parameters are written as:
C ( u ) = i = 0 n N i , p ( u ) ω i P i i = 0 n N i , p ( u ) ω i 0 u 1 S ( u ) = i = 0 n j = 0 m N i , p ( u ) N j , q ( v ) ω i , j P i , j i = 0 n j = 0 m N i , p ( u ) N j , q ( v ) ω i , j 0 u , v 1
where ω is the corresponding weight. The shapes of the domain boundary and the holes inside the boundary can be described precisely by the control points of NURBS.

2.3. Integration Subtraction Technique

The integration is performed in the parameter domain conveniently compared with the physical domain. The stiffness of the j-th background element in the i-th hole is given by
K j = k = 1 n g ω k ( B k T D B k ) J k
where B k is the strain matrix and D is the material constant matrix; n g is the number of Gaussian integral points in the current background element; ω k is the integral weight; and J k is the Jacobian matrix.
The shape and inner holes of the solution domain change during the topology optimization. To solve this problem, the process in conventional methods involves deleting extra meshes, repairing gaps, and reconstructing boundary meshes. It is complex and time-consuming to reconstruct the meshes and the stiffness matrix when changing the inner boundary shape or adding new holes, even though the background meshes are easy to generate. An integration subtraction technique is proposed and the stiffness matrix of the predesigned domain is assumed to be K 0 ; Equation (6) can then be rewritten as
( K 0 K i ) U = F
As shown in Figure 2, when a hole moves from C 1 (red) to C 2 (blue), the homotopic mesh is conveniently applied between the changed regions. Instead of deleting Ω C 1 and adding Ω C 2 , the new global stiffness matrix is obtained easily by adding the stiffness K ( C 2 C 1 ) of the changed region Ω ( C 2 C 1 ) , and the accuracy is satisfied. More details and the numerical results of the integration subtraction are introduced in authors’ previous research [26].

3. Moving Boundary Optimization

Moving boundary optimization (MBO) refers to the presented topology optimization based on the moving boundary meshfree method. The parametric modeling, sensitivity analysis, and optimization strategy are the key points in topology optimization.

3.1. Topological Structure Based on NURBS

The solution domains of the meshfree method are described by several NURBS closed curves. Figure 3 depicts the variation of the solution domain and its boundary from the k-th generation to the k + 1-th generation. The domain Ω 2 stays the same and the domain Ω 3 is a new hole generated during the optimization.
The shape of the NURBS is determined by the coordinates of the control points, so the optimization problem is to minimize the structural compliance by searching for the optimal control points subject to topological constraints.
min P C ( P ) = U T K U s . t .   K U = F V = 1 V 0 i = 1 N Ω i d Ω V ¯ ψ i = j = 1 n Γ j r i j r i j 2 n d Γ j δ , i = 1 , 2 , , n 0
where C ( P ) is the structural compliance; V V ¯ indicates the volume constraint; and ψ i < δ is the flux constraint. The design variables increase or decrease with the variation of the holes and the boundaries.

3.2. Sensitivity Analysis

The positions of control points are the design variables, so the sensitivity analysis involves calculating the derivatives of responses to the coordinates of control points.

3.2.1. Volume Sensitivity

The volume sensitivity of a two-dimensional problem is written as
V Γ P i x = 1 2 j = 1 n g i w j f j P i x g ( u j ) g ( u j ) f j P i x V Γ P i y = 1 2 j = 1 n g i w j f ( u j ) g j P i y g j P i y f ( u j )
where
f j P i x = f j P i x ω P i x ω P i x = N i , p ( u j ) N i , p ( u j ) / ω ( u j ) ω ( u ) ω ( u j ) ω i = g j P i y f j P i x = f j P i x ω P i x ω P i x = N i , p ( u j ) ω ( u j ) ω i = g j P i y
In Equation (12), P i ( x i , y i , ω i ) is a control point of the NURBS curve Γ ( u ) ( Γ ( u ) = ( f ( u ) , g ( u ) ) ) and V Γ is the volume surrounded by Γ ( u ) .

3.2.2. Compliance Sensitivity

The calculation of compliance sensitivity can be divided into two aspects: curve homotopy mapping and point homotopy mapping. Point homotopy requires calculation of the partial derivatives of the strain matrix B , which require calculation of the second derivatives of the shape function. The interpolation accuracy of the second derivatives is low and point homotopy mapping is not recommended. Curve homotopy only needs the integral along the curve and avoids the calculation of partial derivatives, so curve homotopy mapping is more accurate and convenient. Similar to the moving boundary method, the varied stiffness of the global solution domain is approximately regarded as the stiffness of the changed region. When the i-th control point moves Δ P i x along the x-direction, the j-th integral point moves Δ x ( u j ) and the varied stiffness can be simplified to the integral along the curve.
Δ K = Ω B T D B d x d y = Δ Ω B T D B ( g ( u ) ) d u d x Γ B T D B g ( u ) Δ x ( u ) d u = k = 1 n k j = 1 n g B j T D B j ω j ( g ( u j ) ) Δ x ( u j ) = k = 1 n k j = 1 n g B j T D B j ω j ( g ( u j ) ) f j P i x Δ P i x
where n k is the segment number between [ u i , u i + p + 1 ) . The approximate value of the compliance sensitivity is written as
K P i x Δ K Δ P i x = k = 1 n k j = 1 n g B j T D B j ω j ( g ( u j ) ) f j P i x C P i x = U T K P i x U

3.2.3. Topological Sensitivity

The topological sensitivity is the influence of a new hole on the objective and constraints, and the topological sensitivity of a field node is defined as
δ T V = V r = 2 π r δ T C i = U s T B i T D B i U s
where δ T V is the topological sensitivity of the volume constraint and r is the radius of the new hole. δ T C i is the topological sensitivity of compliance of the i-th field point and U s are the displacements of the field points in the support domain surrounding the i-th field point.
When the average value of the topological sensitivities of a field point is larger than a specific threshold value, a new hole is generated at the current field point.

3.3. Optimization Flowchart

The GCMMA is utilized to search for the optimal topology structure [27], and the standard form of a general nonlinear optimization problem is written as
min   f 0 ( x ) + a 0 z + i = 1 m ( c i y i + 1 2 d i y i 2 ) s . t .   g i ( x ) a i z y i 0 i = 1 , , m x j min x j x j max j = 1 , , n y i 0 , z 0 i = 1 , , m
where a i , b i , c i , d i are the coefficients; x is the vector of design variables; f 0 ( x ) is the objective function; g i ( x ) are the constraint functions, including volume and flux constraints, in the topology optimization; and y i , z are the additional design variables.
The optimization flowchart using MBO is shown in Figure 4 and the main procedures include (1) initialization of the design variables and solution domain, (2) construction of GCMMA form of the topology problem, (3) analysis of the solution domain, (4) resolution of the dual problem of GCMMA and the increase of new holes and design variables, and (5) convergence judgment.

3.4. Topology Optimizaiton of MBB Beam

3.4.1. Design Strategy

The MBB beam is a common model for topology optimization. The dimensionless beam with symmetric, simply supported ends is shown in Figure 5 and the unit load is applied at the middle of the top edge.
The constraint is that the volume ratio is less than 50%. The initial number of the control points is 36 and the two coordinates of a control point correspond to two design variables. The control points guarantee that the current solution domain must be within the initial boundary. When the volume ratio is less than 60%, the boundary optimization is completed and the inner structure optimization is carried out. Three holes with radii of 1 are generated at the appropriate positions according to the topological sensitivity. The circumference of each hole is divided into 12 segments and another 12 control points are added to the design variables to optimize the inner boundary, as shown in Figure 6. The initial circular hole is tiny and almost equal to a point. The circular hole gradually develops into an arbitrary shape according to the movement of the control points.
In addition, SIMP and LSM are also performed subject to the same constraints to verify the effectiveness of MBO. The design variables of these two methods are the densities of 2400 elements, and the minimal density is 0.001.

3.4.2. Design Results

The convergence processes of SIMP, LSM, and MBO are compared in Figure 7.
SIMP has the fastest convergence speed, while the convergence speeds of MBO and LSM are similar. The optimal structural compliances and volume ratios of SIMP, LSM, and MBO are compared in Table 1.
The optimal structural compliance of SIMP is the largest. The optimal volume ratio of LSM is larger than the constraint, and the volume ratio cannot be further reduced by adjusting parameters. The MBB beam structure optimization processes in the three methods are shown in Figure 8. SIMP, LSM, and MBO show their particular variation trends in the optimization process, and the final topological structures are similar.
According to the optimization results for the MBB beam, MBO can obtain better structural compliance than SIMP and a better volume ratio than LSM, so MBO shows an advantage in topology optimization.

3.4.3. Influence of Design Variables and Constraint

SIMP and LSM both have the problem of mesh dependence, so the influence of the number of design variables on the three optimization methods should be studied. The mesh density is increased by two times. The numbers of finite elements for SIMP and LSM are increased to 9600 and the initial number of control points for MBO is increased to 68. The optimization results are shown in Table 2.
The optimal topological structures are compared in Figure 9.
The optimization results for SIMP are quite different from the results for the original meshes. Holes appear in the solid beam structure of the original model. Compared with the results for LSM with the original meshes, the optimization is insufficient and many small beams appear, which can impact engineering applications. The optimization results for MBO show little difference from the previous results.
These three optimizations were performed on the same personal computer and the calculation times are compared in Table 3. For an n-dimensional problem, when the mesh density increases by two times, the numbers of design variables for SIMP and LSM increase to 2 n times and the number of design variables for MBO increases to 2 n 1 times. More design variables mean more time for the sensitivity analysis. Judging from the design results and the calculation time, SIMP and LSM are greatly influenced by the mesh density. With the increase in design variables, MBO gradually shows its advantages in optimization efficiency.
The influence of constraint on the three optimization methods is examined next. The optimizations are carried out subject to the constraints of a 60% and 40% volume ratio, respectively. Except for the volume constraint, the object function, design variables, and analysis conditions are the same as those in Section 3.4.1. The optimization results are shown in Table 4.
The optimal topological structures are compared in Figure 10.
The optimization results show that SIMP and LSM are greatly influenced by constraint and the number of design variables. MBO is more stable and has better boundary quality and less structural compliance.

3.5. Discussion

The SIMP, LSM, and MBO methods are compared according to the topology optimization results for the two-dimensional MBB beam. Then, the optimization stability of these three methods is preliminarily explored by adjusting the number of design variables and the volume constraint.
SIMP and LSM are limited by the mesh density and the volume constraint. Different conditions may lead to different results and the structural boundary may appear serrated, which requires manual smoothing. On the other hand, MBO is almost uninfluenced by them and has better optimization stability in simple two-dimensional problems.
MBO comprises shape optimization and topology optimization. For the MBB optimization, MBO can obtain a lighter result with smoother boundaries, which is more suitable for practical manufacturing.
It is more time-consuming for MBO to solve the problem with fewer design variables. However, with the increase in design variables, MBO needs much fewer parameters to describe the structure than FEM-based optimization, which helps MBO save much time in the sensitivity analysis.
The result of MBO can be pre-designed manually. In addition to automatically judging the hole numbers or the hole positions with the algorithm, they can also be specified manually. If the MBB beam is optimized using a symmetric four-hole topological structure, the optimum has a compliance of 48.28 and a volume ratio of 49.9%. The topological structure result is shown in Figure 11.

4. Topology Optimization of CRM Wing Structure

The topology optimization of the CRM wing structure is performed under asymmetric boundary conditions. The wing has a fixed support root and the external loads of the topology optimization are developed in the static aeroelastic analysis. The aeroelastic analysis of the wing structure needs to be analyzed repeatedly because its topological structure varies in each optimization generation.

4.1. Static Aeroelasticity

The static aeroelasticity is analyzed with an improved three-dimensional flight load that combines the high-order panel method and the meshfree method. The high-order panel method is employed for aerodynamic analysis and the meshfree method is employed for static analysis. The difference between the improved and basic three-dimensional flight load methods is that the static analysis is changed from FEM to the meshfree method. The method contains two iterative loops: the inner loop concerns the structural deformation convergence and the outer loop concerns the trim parameter convergence. More details about the basic three-dimensional flight load are introduced in the authors’ previous research [28,29]. The flowchart is shown in Figure 12.
The fight speed is 0.785 Ma at an altitude of 10,000 m and the dynamic pressure is 13,360 Pa. The static aeroelasticity is analyzed at a fixed angle of attack of 2°, so only the inner loop of the three-dimensional flight load is executed. After the static aeroelastic equilibrium is reached, the aerodynamic loads on the high-order aerodynamic panels are applied as the external loads to the wing structure to carry out the topology optimization.
The optimization object is the CRM wing structure, which is a common research model for aerodynamic and aeroelastic research, as shown in Figure 13. The half span is 29.38 m, the root chord is 11.87 m, and the taper ratio is 4.34. The sweep angle of the leading edge is 37.23° and the dihedral angle is 4.85°. The material of the wing is 7075 aluminum alloy, the elastic modulus is 71.7 GPa, the Poisson ratio is 0.33, and the material density is 2700 kg/m3.
The structural model is adjusted to satisfy the requirements of topology optimization. The wing shape remains unchanged and the inner space is filled. Some ribs are reserved to maintain the shape and transfer flight loads. The stiffness matrix of a three-dimensional wing structure is extremely huge, so, limited by the calculation conditions, only the red inner wing segment is optimized and the green segment is used for aeroelastic analysis.

4.2. Design Strategy

The design object is expanded from a two-dimensional beam to a three-dimensional wing. SIMP and MBO are carried out to compare the optimization results. The aerodynamic shape of the CRM wing remains unchanged and only the inner topology structure is optimized.
In the MBO method, the inner holes of the wing structure are controlled by the NURBS surfaces and the hole shapes vary with the coordinates of the control points. A sphere hole, like the earth, is described by 2 poles, 12 longitude lines, and 5 latitude lines. The distribution of the 62 control points is shown in Figure 14a and the initial spherical hole shape with a radius of 1 described by these control points is shown in Figure 14b.
The sphere centers are can be automatically generated by the algorithm or specified manually. In this research, six initial sphere centers are specified at the inner segment of the wing, as shown in Figure 15. The initial sphere holes in Figure 15 are a little exaggerated to display them more clearly. They are tiny and almost point-like. A sphere hole gradually develops into an arbitrary shape with the movement of its control points. The moving direction of a control point is determined by its compliance sensitivity and volume sensitivity. Therefore, the number of design variables is 372 in total.
The optimization flowchart is similar to that described in Section 3.3 and the external loads used in the topology optimization are updated by the three-dimensional flight load in each generation, as shown in Figure 16. The object is the minimal structural compliance and the volume ratio constraint is 50%.

4.3. Design Results

The design results for MBO and SIMP are compared in Table 5.
The final topological structure optimized by MBO is shown in Figure 17a and the basic configuration extracted from the optimization result is shown in Figure 17b. The MBO result tends toward an integral panel wing structure and retains the thick skin. The interior is a honeycomb-like support structure and the support beams are perpendicular to the rib, which is conducive to improving the structural stiffness. The configuration is consistent with the form of the integral panels commonly used in aircraft wing surfaces, and the boundary is relatively smooth.
The final topological structure optimized by SIMP is shown in Figure 18a and, with regard to the beam positions of the CRM wing structure, a double-beam configuration of the inner wing segment with a 50% volume ratio is designed, as shown in Figure 18b.
In the SIMP result, the wing structure changes from a double-beam configuration to a single-beam configuration, and holes appear in the middle of the front beam. The structural compliance of the CRM structure is 9.865 × 107 and the volume ratio is 50.44%. The results show that topology optimization can help reduce structural compliance.

4.4. Discussion

MBO is preliminarily applied to a three-dimensional CRM wing structure and, as a comparison, SIMP is also performed with the same conditions.
The optimization results of MBO and SIMP have similar levels of compliance, but the structural boundary of MBO is smoother and flatter. The clear structure of MBO is more suitable for engineering practice than SIMP.
Compared with SIMP, MBO has poor computational robustness and unstable results. The calculation codes of the meshfree methods need to be further improved and optimized.
For the problem of a three-dimensional complex structure, some particular boundaries, edges, and corners with large curvatures are difficult for the current NURBS curves and surfaces to describe. They should be further designed in the subsequent shape and size optimization.

5. Conclusions

A moving boundary meshfree topology optimization that combines the Galerkin method of weighted residuals and NURBS is presented in this paper.
The proposed method was first applied to the two-dimensional MBB beam. Compared with the results of the FEM-based methods, the moving boundary optimization could obtain more than 3.5% compliance reduction subject to the same constraints. The shape constraints were conveniently applied and the optimized boundary was smoother. The proposed method only needed 2% of the design variables of the FEM-based methods or fewer, which helps reduce the time taken for sensitivity analysis and improves the optimization efficiency. Furthermore, the moving boundary method was almost unaffected by mesh density, the number of design variables, and constraints. The optimization stability makes it more flexible in dealing with two-dimensional problems.
Then, the proposed method was applied to the inner segment of a three-dimensional CRM wing structure, considering the aeroelastic effect. The moving boundary optimization needed 2% of the design variables of the FEM-based method and obtained the prospective compliance subject to the same constraint. The optimal wing structure was clear and close to the double-beam configuration in practical engineering. The proposed method was still efficient for the three-dimensional wing structure.
As a new topology optimization method, there is still a big gap between the proposed method and commercial FEM-based software in terms of robustness and applicability with complex structures. Some regions with large curvatures need to be described more exactly in the meshfree method. The program codes need to be further improved and optimized. The moving boundary technique has great research value and broad development prospects. Our future research will focus on the above key points.

Author Contributions

Conceptualization, X.W., Z.W. (Zhiqiang Wan) and S.Z.; methodology, X.W. and Z.W (Zhiqiang Wan).; software, S.Z.; validation, S.Z.; formal analysis, X.W.; investigation, Z.W. (Zhiqiang Wan); resources, S.Z.; data curation, X.W.; writing—original draft preparation, X.W.; writing—review and editing, S.Z. and Z.W. (Zhiqiang Wan); visualization, X.W.; supervision, Z.W. (Zhiqiang Wan); project administration, X.W.; funding acquisition, Z.W. (Zhi Wang). All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Zhejiang Key Laboratory of General Aviation Operation Technology (General Aviation Institute of Zhejiang JianDe), no.: JDGA2020-4.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhu, J.H.; Zhang, W.H.; Xia, L. Topology Optimization in Aircraft and Aerospace Structures Design. Arch. Computat. Methods Eng. 2016, 23, 595–622. [Google Scholar] [CrossRef]
  2. Aage, N.; Andreassen, E.; Lazarov, B.S.; Sigmund, O. Topology optimization in aircraft and aerospace structures design computational morphogenesis for structural design. Nature 2017, 550, 84–86. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Larsen, S.D.; Sigmund, O.; Groen, J.P. Optimal truss and frame design from projected homogenization-based topology optimization. Struct. Multidisc. Optim. 2018, 57, 1461–1474. [Google Scholar] [CrossRef] [Green Version]
  4. Ghasemi, A.; Elham, A. Efficient multi-stage aerodynamic topology optimization using an operator-based analytical differentiation. Struct. Multidisc. Optim. 2022, 65, 130. [Google Scholar] [CrossRef]
  5. Munk, D.J.; Miller, J.D. Topology Optimization of Aircraft Components for Increased Sustainability. AIAA J. 2022, 60, 445–460. [Google Scholar] [CrossRef]
  6. Gomes, P.; Palacios, R. Aerostructural topology optimization using high fidelity modeling. Struct. Multidisc. Optim. 2022, 65, 137. [Google Scholar] [CrossRef]
  7. Wu, T.Y.; Li, S. An efficient multiscale optimization method for conformal lattice materials. Struct. Multidisc. Optim. 2021, 63, 1063–1083. [Google Scholar] [CrossRef]
  8. López, C.; Baldomir, A.; Hernández, S. The relevance of reliability-based topology optimization in early design stages of aircraft structures. Struct. Multidisc. Optim. 2018, 57, 417–439. [Google Scholar] [CrossRef]
  9. Zhu, B.L.; Zhang, X.M.; Fatikow, S. Structural topology and shape optimization using a level set method with distance-suppression scheme. Comput. Methods Appl. Mech. Eng. 2015, 283, 1214–1239. [Google Scholar] [CrossRef] [Green Version]
  10. Yin, L.; Zhang, F.; Deng, X.W.; Wu, P.; Zeng, H.X.; Liu, M. Isogeometric Bi-Directional Evolutionary Structural Optimization. IEEE Access 2019, 7, 91134–91145. [Google Scholar] [CrossRef]
  11. Jantos, D.R.; Riedel, C.; Hackl, K.; Junker, P. Comparison of thermodynamic topology optimization with SIMP. Contin. Mech. Thermodyn. 2019, 31, 521–548. [Google Scholar] [CrossRef]
  12. Cabral, M.V.; Marques, F.D.; Ferreira, A.J.M. Nonlinear supersonic post-flutter response of two-bay composite laminate curved panels. Compos. Struct. 2022, 286, 115128. [Google Scholar] [CrossRef]
  13. Zhang, H.J.; Wang, D.D. Reproducing kernel formulation of B-spline and NURBS basis functions: A meshfree local refinement strategy for isogeometric analysis. Comput. Methods Appl. Mech. Eng. 2017, 320, 474–508. [Google Scholar] [CrossRef]
  14. Belytschko, T.; Krongauz, Y.; Organ, D.; Fleming, M.; Krysl, P. Meshless methods: An overview and recent developments. Comput. Methods Appl. Mech. Eng. 1996, 139, 3–47. [Google Scholar] [CrossRef] [Green Version]
  15. Hsieh, Y.M.; Pan, M.S. ESFM: An Essential Software Framework for Meshfree Methods. Adv. Eng. Softw. 2014, 76, 133–147. [Google Scholar] [CrossRef]
  16. Chen, J.S.; Hillman, M.; Chi, S.W. Meshfree Methods: Progress Made after 20 Years. J. Eng. Mech. 2017, 143, 04017001. [Google Scholar] [CrossRef] [Green Version]
  17. Garg, S.; Pant, M. Meshfree Methods: A Comprehensive Review of Applications. Int. J. Comput. Methods 2018, 15, 1830001. [Google Scholar] [CrossRef]
  18. Shen, Q.; Ding, R.; Huo, Y.B. Numerical solution of the quasistatic contact problem with the Tresca friction in elastic-viscoplastic materials by the element-free Galerkin method. Eng. Anal. Bound. Elem. 2021, 132, 202–220. [Google Scholar] [CrossRef]
  19. Li, Y.; Liu, G.R. An element-free smoothed radial point interpolation method (EFS-RPIM) for 2D and 3D solid mechanics problems. Comput. Math. Appl. 2019, 77, 441–465. [Google Scholar] [CrossRef]
  20. Dehghan, M.; Narimani, N. The element-free Galerkin method based on moving least squares and moving Kriging approximations for solving two-dimensional tumor-induced angiogenesis model. Eng. Comput. 2020, 36, 1517–1537. [Google Scholar] [CrossRef]
  21. Li, X.L.; Li, S.L. A fast element-free Galerkin method for the fractional diffusion-wave equation. Appl. Math. Lett. 2021, 122, 107529. [Google Scholar] [CrossRef]
  22. Fung, Y.C.; Tong, P.; Chen, X.H. Classical and Computational Solid Mechanics, 2nd ed.; World Scientific: Singapore, 2017; pp. 520–546. [Google Scholar]
  23. Sivaram, S.A.; Vinoy, K.J. Inverse Multiquadric Radial Basis Functions in Eigenvalue Analysis of a Circular Waveguide Using Radial Point Interpolation Method. IEEE Microw. Wirel. Compon. Lett. 2020, 30, 537–540. [Google Scholar] [CrossRef]
  24. Gupta, A.; Verma, S.; Ghosh, A. Static and dynamic NURBS-based isogeometric analysis of composite plates under hygrothermal environment. Compos. Struct. 2022, 284, 115083. [Google Scholar] [CrossRef]
  25. Montemurro, M.; Refai, K. A Topology Optimization Method Based on Non-Uniform Rational Basis Spline Hyper-Surfaces for Heat Conduction Problems. Symmetry 2021, 13, 888. [Google Scholar] [CrossRef]
  26. Liu, Y.Z.; Wan, Z.Q.; Yang, C.; Wang, X.Z. NURBS-Enhanced Meshfree Method with an Integration Subtraction Technique for Complex Topology. Appl. Sci. 2020, 10, 2587. [Google Scholar] [CrossRef] [Green Version]
  27. Li, X.Y.; He, L.L. Shape Optimization Design for a Centrifuge Structure with Multi Topological Configurations Based on the B-Spline FCM and GCMMA. Appl. Sci. 2020, 10, 620. [Google Scholar] [CrossRef] [Green Version]
  28. Liu, Y.Z.; Zhu, S.Y.; Wan, Z.Q.; Yang, C. A High Efficiency Aeroelastic Analysis Method based on Rigid External Aerodynamic Force and Elastic Correction by High-Order Panel Method. In Proceedings of the 55th AIAA Aerospace Sciences Meeting, Grapevine, TX, USA, 9–13 January 2017. [Google Scholar]
  29. Wang, X.Z.; Wan, Z.Q.; Liu, Z.; Yang, C. Integrated optimization on aerodynamics-structure coupling and flight stability of a large airplane in preliminary design. Chin. J. Aeronaut. 2018, 31, 1258–1272. [Google Scholar] [CrossRef]
Figure 1. Discretization of solution domain.
Figure 1. Discretization of solution domain.
Symmetry 14 01154 g001
Figure 2. Integration subtraction when a hole moves from C 1 (red) to C 2 (blue).
Figure 2. Integration subtraction when a hole moves from C 1 (red) to C 2 (blue).
Symmetry 14 01154 g002
Figure 3. Solution domains and boundaries of the k-th and k + 1-th generations.
Figure 3. Solution domains and boundaries of the k-th and k + 1-th generations.
Symmetry 14 01154 g003
Figure 4. Flowchart of topology optimization using MBO.
Figure 4. Flowchart of topology optimization using MBO.
Symmetry 14 01154 g004
Figure 5. Initial structure of MBB beam.
Figure 5. Initial structure of MBB beam.
Symmetry 14 01154 g005
Figure 6. Distribution of control points.
Figure 6. Distribution of control points.
Symmetry 14 01154 g006
Figure 7. Convergence process: structural compliance (left) and volume ratio (right).
Figure 7. Convergence process: structural compliance (left) and volume ratio (right).
Symmetry 14 01154 g007
Figure 8. Optimization processes of SIMP, LSM, and MBO.
Figure 8. Optimization processes of SIMP, LSM, and MBO.
Symmetry 14 01154 g008
Figure 9. Optimal topological structures of SIMP, LSM, and MBO.
Figure 9. Optimal topological structures of SIMP, LSM, and MBO.
Symmetry 14 01154 g009
Figure 10. Optimal topological structures of different constraints.
Figure 10. Optimal topological structures of different constraints.
Symmetry 14 01154 g010
Figure 11. Optimal topological structure with four holes.
Figure 11. Optimal topological structure with four holes.
Symmetry 14 01154 g011
Figure 12. Meshfree method-based three-dimensional flight load.
Figure 12. Meshfree method-based three-dimensional flight load.
Symmetry 14 01154 g012
Figure 13. Structural model of CRM wing.
Figure 13. Structural model of CRM wing.
Symmetry 14 01154 g013
Figure 14. Control points along longitude and latitude lines (a) and initial sphere hole (b).
Figure 14. Control points along longitude and latitude lines (a) and initial sphere hole (b).
Symmetry 14 01154 g014
Figure 15. Initial positions of the spherical hole centers.
Figure 15. Initial positions of the spherical hole centers.
Symmetry 14 01154 g015
Figure 16. Flowchart of aeroelastic topology optimization.
Figure 16. Flowchart of aeroelastic topology optimization.
Symmetry 14 01154 g016
Figure 17. Topological configuration (a) and extraction of structural feature (b).
Figure 17. Topological configuration (a) and extraction of structural feature (b).
Symmetry 14 01154 g017
Figure 18. SIMP result (a) and designed CRM wing structure (b).
Figure 18. SIMP result (a) and designed CRM wing structure (b).
Symmetry 14 01154 g018
Table 1. Optimization results.
Table 1. Optimization results.
SIMPLSMMBO
Compliance49.3547.3947.43
Volume ratio50.0%50.4%50.0%
Table 2. Optimization results (increased design variables).
Table 2. Optimization results (increased design variables).
SIMPLSMMBO
Compliance50.0149.3848.06
Volume ratio50.0%50.1%49.6%
Table 3. Optimization generation and calculation time.
Table 3. Optimization generation and calculation time.
Coarse MeshFine Mesh
VariableTime (s)GenerationVariableTime (s)Generation
SIMP240010.7319600699.1100
LSM240017.0599600293.948
MBO72309.0100136326.0100
Table 4. Optimization results (different constraints for the volume ratio).
Table 4. Optimization results (different constraints for the volume ratio).
ConstraintSIMPLSMMBO
Compliance40%59.4657.2758.98
Volume ratio40.0%39.9%39.9%
Compliance60%43.5242.5042.28
Volume ratio60.0%60.4%59.9%
Table 5. Optimization results for MBO and SIMP.
Table 5. Optimization results for MBO and SIMP.
Design VariableComplianceVolume RatioGeneration
MBO3729.426 × 10749.57%65
SIMP18,0009.420 × 10750.00%27
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wang, X.; Zhang, S.; Wan, Z.; Wang, Z. Aeroelastic Topology Optimization of Wing Structure Based on Moving Boundary Meshfree Method. Symmetry 2022, 14, 1154. https://doi.org/10.3390/sym14061154

AMA Style

Wang X, Zhang S, Wan Z, Wang Z. Aeroelastic Topology Optimization of Wing Structure Based on Moving Boundary Meshfree Method. Symmetry. 2022; 14(6):1154. https://doi.org/10.3390/sym14061154

Chicago/Turabian Style

Wang, Xiaozhe, Shanshan Zhang, Zhiqiang Wan, and Zhi Wang. 2022. "Aeroelastic Topology Optimization of Wing Structure Based on Moving Boundary Meshfree Method" Symmetry 14, no. 6: 1154. https://doi.org/10.3390/sym14061154

APA Style

Wang, X., Zhang, S., Wan, Z., & Wang, Z. (2022). Aeroelastic Topology Optimization of Wing Structure Based on Moving Boundary Meshfree Method. Symmetry, 14(6), 1154. https://doi.org/10.3390/sym14061154

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop