1 Introduction

As a promising approach for simulating complex fluid dynamic problems, the lattice Boltzmann method (LBM) has been adopted successfully in many different flow problems, such as incompressible flow, porous-media flow, and multi-phase flow^[1-4]. However, most studies focused on the LBM calculation of the flow field, and the LBM has been rarely applied to the inverse design problems.

The adjoint method, which was presented by Jameson^[5] in 1980s, is one of the most widely adopted airfoil optimization approaches in recent years. This approach can dramatically reduce the cost of the computation since the computational expense incurred in the calculation of the complete gradient is effectively independent of the number of design variables. During the last two decades, this method has been developed rapidly. For example, there have been many studies including the design of airfoils, wings, and the whole aircraft using continuous and discrete adjoint methods^[6-9]. Anderson and Venkatakrishnan^[10] applied this approach on unstructured grids. Some other studies adopted this approach in the optimization of aerodynamic noise and fluid structure interaction.

Although the adjoint approach and the LBM are both research hot spots in recent years, there are only few cross-over studies which combine them together. Pingen et al.^[11], Evgrafov et al.^[12], and Pingen et al.^[13-14] studied the topological optimization based on the adjoint approach and the LBM. Yaji et al.^[15] also studied the level set optimization method based on the LBM. Ngnotchouye et al.^[16] studied the one-dimensional Euler equation based on the LBM and adjoint approach. Unfortunately, this method cannot be generalized to two-dimensional problems due to the difficulties of boundary treatment.

To the authors' knowledge, the LBM based adjoint approach has not been successfully applied to study aerodynamic shape optimization design problems. This study is the first effort to study the design problem with the combination of the LBM and the adjoint method. There will be clear advantages if the LBM is adopted in the optimal design process. They mainly manifest in the unity of the governing equation. By changing the lattice Boltzmann model, different physical problems can be solved. The optimal design approach has stronger reliability and transportability. Nevertheless, we have to face the singularity problem of the adjoint equation boundary conditions in the design process, which may limit the application of the LBM in optimal design to some extent.

The paper is organized as follows. Section 2 briefly describes the LBM theory and the basic model. Section 3 presents the LBM on a body-fitted grid, and the basic theory of the GILBM (generalized form of interpolation supplemented LBM) is introduced. Section 4 gives the optimization method based on the LBM and the adjoint approach in detail, and a test case is given in Section 5 to show feasibility of the new method. Conclusions are given in Section 6.

2 LBM 2.1 Governing equation and basic model

The governing equation of the LBM, the lattice Boltzmann equation (LBE), is a special discrete formulation of the Boltzmann equation. The LBE is described as

(1)

where i represents the particles with different discrete velocities, f_i is the distribution function, c_i=(c_i1, c_i2) is the discrete velocity, Ω_i is the collision operator, τ is the relaxition time, and f_i^eq is the equilibrium distribution function for the discrete velocity c_i.

The evolution equation can be described as

(2)

The LBE consists of collision and advection steps as follows:

(3)

(4)

where f_i^* is the post-collision distribution function.

In this paper, the nine-velocity model^[17] is adopted as the discrete velocity model. The particle velocity vector set is

(5)

The equilibrium distribution functions for the discrete velocities c_i are defined as

(6)

The constants w_i and c_s are

The particle equilibrium distribution functions need to satisfy the following relations:

(7)

(8)

(9)

where u is the fluid velocity vector, u²=u· u, ρ is the fluid density, p is the pressure, and δ_ab is the Kronecker delta.

According to the collision invariants, we have

(10)

(11)

The macro stress can be calculated through the gas state equation: p = ρc_s².

2.2 Boundary treatment

In this paper, the non-equilibrium extrapolation method which was proposed by Guo et al.^[18] is adopted to treat the boundary conditions. The basic idea of this method is to decompose the distribution function at the boundary nodes into the equilibrium and non-equilibrium parts, and then approximate the non-equilibrium part with a first-order extrapolation of the non-equilibrium part of the distribution at the neighbouring fluid nodes.

3 GILBM

The traditional way for solving the LBM is based on the Cartesian grid. The adoption of Cartesian grid makes the LBM have the advantages of easy implementation and good parallel performance. However, it brings a lot of disadvantages to the LBM, such as huge computation, unadjustable calculation accuracy, and difficulty in dealing with the complicated boundary. In order to solve these problems, many researchers have studied the LBM on non-uniform mesh. Among methods that have already been presented, the GILBM^[19] is the most representative one, which combines the LBM with body-fitted meshes.

Based on the computational method in the traditional finite volume method, the concept of computational planes is introduced in the GILBM. The collision and the advection happen on the physical planes, and calculate on the computational planes. The governing equations (3) and (4) on orthogonal coordinates are transformed into generalized coordinates. The transformation of the collision term equation into generalized coordinates is simple, because only local information of the current grid node is required. Thus, by replacing x and y in Eq. (3) to ξ and η, the collision term calculation on generalized coordinates is obtained as

(12)

After the collision step, the advection term calculation is performed. Transformation of the advection term equation requires some discussion. Figure 1 shows the discrete velocities on the physical and computational planes. They are constant vectors on the physical plane, but not constant on the computational plane. Through the coordinate transformation, the discrete velocities e_i are defined as

(13)

where J is the Jacobian transformation and given as .

After that, integration of the contravariant velocity over time step ∆t would be performed to calculate the position where the distribution function comes from.

The transformed advection equation becomes

(14)

4 Continuous adjoint approach based on GILBM 4.1 Derivation of continuous adjoint equation

In the process of design optimization, an appropriate form of the objective function must be chosen. In this paper, the objective function is defined as below to carry on the theoretical derivation,

(15)

where p_d is the desired pressure.

A variation in the airfoil shape will cause a variation in the objective function,

(16)

As mentioned above, the governing equation in the computational plane for the steady flow can be described as

(17)

where e_i1 and e_i2 present the collision velocities in the computational plane.

Using the Lagrange multiplier method to add the governing equation to the variation of the objective function as a constraint, we have

(18)

where λ_i is the Lagrange multiplier.

Dealing with the third term of Eq. (18), we have

(19)

(20)

According to the Gauss theorem, we have

(21)

where n_i is the component of a unit vector normal to the boundaries. If the mesh is assumed as O-type mesh, there is no boundary integral in the η-direction. On the airfoils n₁ = 0 and n₂ = −1, we have

(22)

Then,

(23)

In order to eliminate the influence of flow field changes on objective function gradients, we make the above formula be 0,

(24)

Equation (24) is the adjoint equation for the lattice Boltzmann equation. Adding the time term and transforming it to the physical plane, we have

(25)

The gradient of the objective function is

(26)

4.2 Singularity in derivation of boundary conditions and its solution

According to the gas state equation, we have

(27)

Introducing this equation to Eq. (26), we have

(28)

Let the coefficient of δf_i be 0,

(29)

When e_i2=0, the right-hand of Eq. (29) is zero, while the left-hand would not always be zero. The singularity problem appears. Therefore, in this way, we cannot handle the boundary conditions of adjoint equation.

Take into account that the LBM is a mesoscopic method. Its macroscopic physical character is a combination of distribution functions. We hope to find another way to represent the pressure to solve the problem. After study, under the assumption of a first approximation for the equilibrium distribution function on the boundary, we propose a new method for boundary treatment which can avoid the problem mentioned above. Specific processes are as follows.

According to the relations, we have

(30)

(31)

(32)

We also know that u=0 and v=0 on the airfoil. Thus,

(33)

After the flow becomes steady, we assume

(34)

Therefore,

(35)

Then, the variation of the objective function is

(36)

We have the boundary condition

(37)

Let

Then, the boundary condition is

(38)

The variation of the objective function is

(39)

4.3 Numerical implementation

The accurate solution of the adjoint equation is very critical in the optimization procedure. According to the definition of collision operator, the adjoint equation can be written as follows:

(40)

If we define

(41)

then the adjoint equation changes into

(42)

Comparing Eq. (1) and Eq. (42), it is obvious that they have the same formulation. There-fore, we can use the same calculation method to solve the adjoint equation.

4.4 Implementation of optimization procedure

The procedures can be summarized as follows:

Step 1  Solve the flow field governing equation via the GILBM.

Step 2  Solve the adjoint equation for adjoint variables.

Step 3  Disturb the mesh, and then calculate the gradient of the objective function.

Step 4  Use the optimization algorithm to obtain a new profile.

Step 5  Return to Step 1 until the convergence criterion is achieved.

5 Test case

In the test case, the NACA0012 airfoil is modified to achieve the pressure distribution of the NACA0010 airfoil. These two airfoils are both symmetrical airfoils, and the maximum relative thickness differs 2%. Figure 2 is the mesh used for optimization procedure.

The computations are performed at a free stream velocity of 0.1c, and the Reynolds number is 100. There are 18 points up and down along the surface of the airfoil chosen as the design variables. Hicks-Henne bump functions are employed to modify the airfoil shape. In order to validate the accuracy of gradient value by the approach mentioned above, the results are compared with the gradient information obtained by the traditional finite difference method. As can be seen from Fig. 3, the results of adjoint method show consistence with the finite difference method.

After 75 design cycles, the solution converges. Figure 4 is the objective function value changing in the design process. We can see that the value stays steady at last. The objective function reduces from 4.48× 10^-6 to 4.75× 10^-8, and the shape and pressure distributions of the final airfoil match with the target very closely, as we can see from Figs. 5 and 6. Figures 7 and 8 show the contours of distribution functions f_i and Lagrange multipliers λ_i after ten design iteration steps.

6 Conclusions

This paper studies the LBM based on body-fitted grid, and combines this method with the adjoint approach for optimizing the aerodynamic characters of the airfoil. For a given form of the objective function, the theoretical derivation of adjoint equation and its boundary conditions has been given. The adjoint equation of LBE is solved by the LBM. Although the boundary treatment proposed in this paper is under the first approximation to the distribution function, the optimized case results show that the approach is practical and effective.

Through theoretical derivation, numerical implementation, and the test case, this paper proposes a new quick and efficient optimization design method.