PHYSICAL REVIEW E 91, 013303 (2015)

Extraction of macroscopic and microscopic adjoint concepts using a lattice Boltzmann method and discrete adjoint approach Mohamad Hamed Hekmat* and Masoud Mirzaei† Center of Excellence for Design & Simulation of Space Systems, K. N. Toosi University of Technology, Tehran, Iran (Received 27 January 2014; revised manuscript received 30 September 2014; published 7 January 2015) In the present research, we tried to improve the performance of the lattice Boltzmann (LB) -based adjoint approach by utilizing the mesoscopic inherent of the LB method. In this regard, two macroscopic discrete adjoint (MADA) and microscopic discrete adjoint (MIDA) approaches are used to answer the following two challenging questions. Is it possible to extend the concept of the macroscopic and microscopic variables of the flow field to the corresponding adjoint ones? Further, similar to the conservative laws in the LB method, is it possible to find the comparable conservation equations in the adjoint approach? If so, then a definite framework, similar to that used in the flow solution by the LB method, can be employed in the flow sensitivity analysis by the MIDA approach. This achievement can decrease the implementation cost and coding efforts of the MIDA method in complicated sensitivity analysis problems. First, the MADA and MIDA equations are extracted based on the LB method using the duality viewpoint. Meanwhile, using an elementary case, inverse design of a two-dimensional unsteady Poiseuille flow in a periodic channel with constant body forces, the procedure of analytical evaluation of the adjoint variables is described. The numerical results show that similar correlations between the distribution functions can be seen between the corresponding adjoint ones. Besides, the results are promising, emphasizing the flow field adjoint variables can be evaluated via the adjoint distribution functions. Finally, the adjoint conservative laws are introduced. DOI: 10.1103/PhysRevE.91.013303

PACS number(s): 02.70.Ns, 05.10.−a

I. INTRODUCTION

In recent years, the use of adjoint methods have become popular in many fields of computational sciences and engineering, particularly shape optimization problems in fluids engineering, such as inverse pressure design in steady potential flows [1], optimum shape design based on the drag minimization in steady [2–4] or unsteady [5] transonic flows, and aerodynamic shape optimization of complex threedimensional (3D) configurations in micro- [6] and large [7,8] scales. These methods have also been used in optimal control theory for controlling physics-based fluid simulations [9–11], uncertainty or sensitivity analysis in many fields of engineering surrounding the whole of the experimental data processing and many computational modeling and process simulation activities [12–14], and data assimilation for oceanographiccoastal applications [15,16] or estimating and controlling models error [17]. Generally, in the adjoint methods, a system is modeled as an algebraic or differential equation, then the adjoint equation of the original problem is derived using the Lagrange multipliers method, and finally both the equations are solved. From the solution of the adjoint equation, the Lagrange multipliers (adjoint variables) are calculated and then the gradients of a cost function with respect to all of the design variables at once. The adjoint methods can be broken down into two categories: continuous and discrete. In the continuous adjoint approach, the adjoint equation is derived by the differential forms of the governing equation and cost function. Variations of the cost function and the governing equation with respect to the state variables and the design variables are combined through the use of adjoint variables. Via

* †

[email protected] Corresponding author: [email protected]

1539-3755/2015/91(1)/013303(15)

gathering the terms associated with the variation of the state variables, the adjoint equation and its boundary conditions are derived. The terms associated with the variation of the design variables produce the cost function gradient vector. The governing equation together with the adjoint equation must finally be discretized using appropriate numerical methods. But, in the discrete version, the adjoint equation is directly extracted from the discrete form of the governing equation and cost function. The discrete adjoint equation is derived by collecting together all the terms multiplied by the variation of the discrete state variables. For more details, the reader is referred to Jameson [1], Nadarajah [2,5], and Giles and Pierce [3]. Although the adjoint methods have a prosperous past in computational fluid dynamics, they have three major drawbacks. First, they may be challenging when we work with the full Navier-Stokes (NS) equations. This leads usually to the complexity both in the mathematical formulation of the adjoint equation and the corresponding boundary conditions. Second, the solution of an adjoint equation is generally a computationally expensive task when the original problem is time dependent [18], since one needs to solve the adjoint equation in reverse time from a final flow solution. Thus the entire flow history should be stored which needs potentially huge memory requirements, and then integrating the adjoint equation using backwards in time also needs powerful processors [3,18]. Finally, the computational costs of the adjoint methods increase proportional to the number of constraints [5]. Hence these methods are not efficient for optimization problems with a high number of constraints. Such problems are the motivation for employing the lattice Boltzmann (LB) method (due to its simplicity, inherent use of immersed boundary methods, parallelizability, and general flexibility) instead of conventional NS-based computational methods in numerical optimization problems using adjoint methods [19]. Tekitek et al. [20] performed the microscopic discrete adjoint

013303-1

©2015 American Physical Society

MOHAMAD HAMED HEKMAT AND MASOUD MIRZAEI

PHYSICAL REVIEW E 91, 013303 (2015)

(MIDA) approach via the LB method. They extracted the MIDA equation based on the LB equation to identify the unknown parameters of the LB method. In addition, Pingen et al. [19,21] explored the ability of the incompressible LB method in topology optimization problems. Their major work was focused on the extraction of the steady MIDA equation and the cost function gradient using the steady LB equation and its implicit solution. Moreover, they examined the performance of mathematical algorithms (e.g., iterative and direct algorithms) to solve these equations implicitly. Ngnotchouye et al. [22] implemented the microscopic continuous adjoint (MICA) method using a compressible LB model. They investigated ability of the compressible LB method in the inverse design of a one-dimensional inviscid and compressible flow field. Recently, Hekmat and Mirzaei [23] introduced the LB-based discrete adjoint method in two versions, macroscopic and microscopic, to decrease the implementation cost of the standard discrete adjoint approach. They focused on the derivation of the adjoint equations and investigated the numerical efficiency and computational cost of the two proposed approaches in the optimization problems of unsteady flow fields. Even though the LB-based adjoint approaches proposed in the past few years have tackled flow optimization problems successfully, their deficiency in flow sensitivity analysis problems is yet a major challenge. It is due to the inability in providing a definite framework such that the complete attainment to all sensitivities of a cost function to flow field variables requires an extensive algebraic handling. This lack is related to solving the adjoint equation in a microscopic level and performing the flow sensitivity analysis in a macroscopic one. This incompatibility decreases the efficiency of the LB-based adjoint method in complicated sensitivity analysis problems. In this research, the mesoscopic inherent of the LB method is utilized for addressing that problem. In this regards, the discrete adjoint method is used in two macroscopic and microscopic versions based on the LB method. In the macroscopic discrete adjoint (MADA) method, the macroscopic properties and the governing equations (conservative laws) of the flow field are used to access the MADA equation. In another one, the MIDA equation is achieved by considering the microscopic properties and the governing equations (LB equation). Then, the results of the implementation of the two methods for an inverse design problem of flow field are reported and the accuracy of the cost function gradients obtained from them is evaluated in comparison with the results of the central finite difference approach. However, the main focus of the study is on the remarkable relationships between the macroscopic and microscopic adjoint variables. In this regard, the inherent of the adjoint variables is investigated. Also, like the conservative laws in the LB method which express the relations between the macroscopic and microscopic properties of the flow field, the comparable relations between the corresponding adjoint ones are generally derived. This accomplishment can facilitate the sensitivity analysis procedure of flow field variables using the MIDA method. In other words, similar to the numerical procedure of the LB method, in which the LB equation is first solved in the microscopic level and then the distribution functions are reached and, finally, the flow field variables are calculated through the conservative laws, the MIDA equation is solved in the same level to obtain the adjoint distribution

functions and, lastly, the flow field adjoint variables can be evaluated by the adjoint conservative laws. Furthermore, due to the similarity between the MIDA and LB equations, the similar numerical method can be used to solve them. This paper is organized as follows. At first, Sec. II gives an overview of the LB method for the dynamics of incompressible fluid flows. The general formulation of an optimization problem is introduced in Sec. III. Section IV describes the derivation procedure of the general MADA and MIDA equations, briefly. In Sec. V some remarks about the present numerical implementation are provided. Finally, Sec. VI begins with validating the implementation and also includes a comprehensive discussion of the obtained numerical results from the MADA and MIDA methods. II. LATTICE BOLTZMANN METHOD

In recent years, there has been rapid progress in developing the LB method for solving a variety of fluid dynamic problems such as numerical simulation of internal and external flows [24–26], turbulent flows [27], flows through porous materials [28], fluid-structure interaction [29], convection heat transfer [30,31], non-Newtonian flows [32,33], multiphase flows [34,35], and further aerodynamic analysis of flows [36,37] and newly numerical optimization [19–23]. Unlike the conventional numerical methods for fluid flow problems which are based on the discretized continuous macroscopic (NS) equations, the LB method is based on the microscopic kinetic (LB) equation [38]. There are two ways to extract the LB equation. The first is to use the lattice gas automata (LGA), proposed by Frisch et al. [39,40] and Wolfarm [41], and replacement of the average of the Boolean occupation number with the particle distribution function and also negligence of individual particle motion and particle-particle correlations in the kinetic equations [38,42]. The second is direct extraction of the LB equation from the Boltzmann transport equation with Bhatnagar-Gross-Krook (BGK) approximation [43–45]. For a system without an external force, the Boltzmann-BGK equation can be written as 1 ∂f + c · ∇f = (f eq − f ), ∂t λ

(2.1)

where f = f (x,c,t) is the density distribution function, c is the particle velocity vector, x is the location vector in physical space, t is time, and the dot (·) denotes the scalar product. Here, the local equilibrium distribution function is denoted by f eq and λ is called the relaxation time. The right-hand side of the above equation represents the BGK collision operator [46]. Now, Eq. (2.1) can be first discretized with respect to velocity using a finite set of velocity vectors ci without affecting the conservative laws [43,44] and, then, with respect to space and time using a Euler time step in conjunction with an upwind spatial discretization, as shown by Sterling and Chen [45]: fi (x + ci t,t + t) = (1 − β)fi (x,t) + β fi eq (x,t). (2.2) In the above equation, ci is the particle velocity vector belonging to some discrete set depending on the lattice arrangement, fi is the distribution function associated with the ith discrete velocity ci , fi eq is the corresponding local

013303-2

EXTRACTION OF MACROSCOPIC AND MICROSCOPIC . . .

PHYSICAL REVIEW E 91, 013303 (2015)

equilibrium distribution function in the discrete velocity space, β = t/λ is the inverse of the dimensionless relaxation time, ci t is the lattice spacing, and t is the time step. To facilitate the implementation of the LB method, it is common to chose time step and lattice spacing equal to 1 (c = 1). The beauty of Eq. (2.2) lies in its simplicity and can be applied for many physics by simply specifying a different equilibrium distribution function. Moreover, adding a source (force) term to the above equation is straightforward. For a system with an external force, the LB equation (2.2) can be rewritten as the following vector form: eq

Ft+1 = (1 − β)Ft + β Ft + V ,

(2.3)

where Ft ≡ F (x,t) = [f0,t , . . . ,fm −1,t ]T ∈ Rm is the distribution functions (microscopic properties) vector at time t and V = [V0 , . . . ,Vm−1 ]T ∈ Rm is the force vector. m refers to the number of discrete velocities in lattice space that depends on the lattice arrangement. In an incompressible flow, ideally, so that the divergence of the velocity is zero, the density is assumed equal to the initial density, i.e., ρt = ρ0 . Therefore, assuming an order O(Ma2 ) density fluctuation, the equilibrium distribution function can be derived using Taylor series expansion of the MaxwellBoltzmann equilibrium distribution and ignoring O(Ma3 ) terms or higher, as shown by He and Lou [43]:

surface points) and, finally, a set of constraints that must be satisfied during the optimization process such as the flow field governing equation. On the whole, an optimization problem in a 2D unsteady and incompressible flow field can be expressed in the following general form: minI (θt ,κ) such that Rt (θt ,κ) = 0,

αi,t = 1 +

3 (ci · ut ) 9(ci · ut )2 3u2t + − , c2 2c4 2c2



(y)

Moreover, macroscopic pressure pt ≡ p(x,t) is related to density using the ideal gas equation of state: pt = cs2 ρt ,

(2.6)



I (θt ,κ) = where

tf 

It (θt ,κ) =

t=0

 x

t=0

Itx (θt ,κ),

(3.3)

x

IV. DISCRETE ADJOINT APPROACH BASED ON THE LATTICE BOLTZMANN METHOD

Generally, there are two ways to extract the adjoint equation. The first way is using the Lagrange multiplier viewpoint. In this scheme, the Lagrange multipliers are the adjoint variables and are introduced into an augmented cost function [1]. In a recent paper on development of discrete adjoint approach based on the LB method by the authors [23], the MADA and MIDA equations were derived in detail using this scheme. The second way is derivation of the adjoint equation using the duality viewpoint [47]. In this section, the mathematical derivation of the MADA and MIDA equations and the cost function gradient vector based on the LB equation is briefly presented using the similar procedure reported by Giels and Pierce [3]. As a brief introduction to the duality viewpoint, suppose matrix A and vectors g and c are known. One may wish to evaluate gT b

Optimization problems are made up of three basic ingredients. First, a cost function which we want to minimize or maximize. For instance, in the inverse design of a flow or shape, the total deviation of the current and target field pressure distribution might be minimized or in an aerodynamic shape optimization problem, drag, lift, or their ratio can be considered as the cost function. Second, a set of unknown variables that are referred to as design variables which affect the value of the cost function and we try to find the values of them that minimize or maximize the cost function. For example, in the shape inverse design or aerodynamic shape optimization problems, design variables are the parameters that define the shape (such as

tf  

is performed over all the lattice points of .

where cs is the lattice speed of sound. III. GENERAL OPTIMIZATION PROBLEM STATEMENT



wherein It is the cost function at time t. Using the numerical integration rule and supposing that x = t = 1, the cost function (3.2) can be evaluated as

(2.4)

T and ρt ≡ ρ(x,t) are the where ut ≡ u(x,t) = [u(x) t ,ut ] macroscopic velocity vector and density, respectively. Also, wi is the lattice weight parameter in the ith lattice direction that depends on the lattice arrangement. The macroscopic properties such as density ρt and velocity vector ut can be evaluated by taking statistical moments of the distribution function, leading to the following equations:  1  ρt = fi (x,t), ut = ci fi (x,t). (2.5) ρt i i

(3.1)

where θt is the r-dimensional state variables vector at time t, κ is the q-dimensional design variables vector, Rt is the governing equation of the flow field, and I is the value of cost function at the time interval ≡ [0,tf ] that can be written in the general form as    It (θt ,κ) dt = Itx (θt ,κ)dx dt, (3.2) I (θt ,κ) =

eq

fi,t = ρ wi αi,t ,

t ∈ ,

θt ∈R r ,κ∈R q

such that Ab = c.

(4.1)

The direct way to solve the main problem (4.1) is first evaluating b and then g T b. As an alternative, the following problem can be considered as a dual problem of (4.1): such that AT s = g.

sT c

(4.2)

Note that s c = s (Ab) = (A s) b = g b. The advantage of resorting to solving the dual problem becomes obvious when there is an aggregate of vectors bi and ci . In this case, clearly the dual problem is simpler to solve in ultimately evaluating the quantity g T b. Now, recall the present optimization problem (3.1). To solve this problem using the adjoint method, the principle step is evaluating the gradient vector of the cost function with respect to the design variables vector. Since the state variables vector

013303-3

T

T

T

T

T

MOHAMAD HAMED HEKMAT AND MASOUD MIRZAEI

PHYSICAL REVIEW E 91, 013303 (2015)

depends on the design variables vector implicitly, the local gradient vector of cost function (3.3) at time t can be evaluated as follows: Gtx = ∇κ Itx =

∂ItxT ∂I T ∂θt + tx , t = 1, . . . ,tf . ∂κ ∂θt ∂κ

(4.3)

The main point in the adjoint method is eliminating the dependency of cost function gradient vector (4.3) on the changes of the state variables vector by introducing the adjoint equation. The governing equation in Eq. (3.1) can be declared based on the LB method in a general compact vector form as follows: R = − ( ,κ) = 0,

( ,κ) = [ 1 (θ0 ,κ), . . . , tf (θtf −1 ,κ)]T ,

(4.4)

where = [θ0 , . . . ,θtf ]T and R = [R1 , . . . ,Rtf ]T . Differentiating constraint (4.4) gives a linear constraint on the derivative matrix ∂ /∂κ. Thus, to remove the contributions of δ from the change associated directly with δκ in the variation of the cost function, the main problem can be considered as follows:   ∂ItxT ∂ ∂ ∂ ∂ such that I − = . (4.5) ∂ ∂κ ∂ ∂κ ∂κ As it was demonstrated in the above, the dual problem of (4.5) can be considered as   ∂ T ∂Itx T ∂ such that I − , (4.6) =  ∂κ ∂ ∂ where  = [ψ1 , . . . ,ψtf ]T is the dual or adjoint variables vector. Using the constraint in problem (4.6), the adjoint terminal condition ∂Itf x ψ tf = (4.7) ∂θtf

are corresponding to the flow field variables (macroscopic (y) (1) (3) properties) ρt , u(x) t , and ut . Therefore, ψt to ψt can be now referred to as the macroscopic adjoint variables. The reason for this name will be explained later. Accordingly, Eq. (4.8) can here be called the MADA equation [23]. B. Microscopic discrete adjoint approach

In this approach, the LB equation (2.3) is considered as the constraint and consequently t ≡ (x,t) = [ 0,t , . . . , m−1,t ]T ∈ Rm is defined based on Eq. (2.3). Moreover, the distribution functions (microscopic properties) vector Ft is chosen as the state variables vector θt and ψt ≡ ψ(x,t) = [ψ0,t , . . . ,ψm −1,t ]T ∈ Rm . Here, the adjoint variables ψ0 to ψm−1 are corresponding to the distribution functions (microscopic properties) f0 to fm−1 in various lattice directions. Therefore, ψ0 to ψm−1 can be now referred to as the microscopic adjoint properties or the adjoint distribution functions. The reason for this name will be explained later. Correspondingly, Eq. (4.8) can here be called the MIDA equation [23]. V. NUMERICAL IMPLEMENTATION

In this section, we briefly summarize the key components of the numerical implementation of the proposed optimization methods. This includes introducing the considered flow field, boundary conditions, and lattice arrangement, selecting the cost function and design variables, and determining the sensitivities of the state variables and cost function using both the MADA and MIDA approaches. Then, the optimization algorithm is introduced and, after that, the procedure of numerical solution of MADA and MIDA equations and extension of the proposed approaches to more realistic situations are briefly described. Finally, comments on the cost of the present methods are provided.

and the unsteady discrete adjoint equation ψtT

∂ItxT T ∂ t+1 = + ψt+1 , t = tf − 1, . . . ,1 ∂θt ∂θt

A. Boundary conditions and lattice arrangement

(4.8)

are attained. Finally, the second term on the right-hand side of Eq. (4.3) can be replaced with the dual form given in Eq. (4.6): ∂ t ∂ItxT + ψtT , t = 1, . . . ,tf . (4.9) ∂κ ∂κ According to Eq. (4.9), Gtx is independent of δθ and consequently, for a large number of design variables, the gradient vector can be evaluated just via one flow field solution in addition to one adjoint solution in each optimization cycle. Gtx = ∇κ Itx =

In this research, two presented approaches are implemented for a 2D unsteady laminar Poiseuille flow in a periodic channel with uniform and constant body forces. For this purpose, the periodic boundary condition is applied to nodes on the left-inlet and right-outlet boundaries of the fluid. Moreover, the no-slip boundary condition is performed in nodes on the upper and lower solid boundaries by the complete bounce-back scheme. Moreover, to simulate the flow field by the LB method, the two-dimensional, nine-velocity D2Q9 lattice model [48] is used (see Fig. 1). 6

2

5

A. Macroscopic discrete adjoint approach

Here, the conservation equations (2.5) are chosen as the constraints in the optimization problem. Thus t ≡ (x,t) = (2) (3) 3 [ (1) t , t , t ] ∈ R and its components are determined via the mass, x- momentum, and y-momentum conservative laws, respectively. Moreover, in this approach, the state variables vector θt is the flow field variables (macroscopic properties) vector ωt . Also, the adjoint variables vector is defined as ψt ≡ ψ(x,t) = [ψt(1) ,ψt(2) ,ψt(3) ] ∈ R3 so that its components

0

3

7

4

1

8

FIG. 1. (Color online) Lattice arrangement of the D2Q9 model.

013303-4

EXTRACTION OF MACROSCOPIC AND MICROSCOPIC . . .

PHYSICAL REVIEW E 91, 013303 (2015)

These velocity vectors ci = [ci,x , ci,y ]T can be written as ⎧ ⎪ ⎨(0,0), ci = (cos[(i − 1)π/2], sin[(i − 1)π/2])c, ⎪ ⎩√ 2(cos[(i − 5)π/2 + π/4], sin[(i − 5)π/2 + π/4])c,

For the D2Q9 model, the lattice weight parameters are as ⎧ i = 0, ⎨4/9, i = 1,2,3,4, (5.2) wi = 1/9, ⎩1/36, i = 5,6,7,8. For low Mach number (Ma  1) or low-density fluctuation (ρt /ρt  1) flows, it is possible to recover the macroscopic incompressible NS equations from the LB equation using the Chapman-Enskog expansion and choosing the kinematic viscosity as [48] ν=

x 2 1 (2 − β) . 6β t

(5.3)

In this study, the flow inverse design problem is employed to implement the proposed MADA and MIDA approaches. Hence the local cost function at time t can be described as 2 σ 1 (5.4) Itx = ωt − ωtd  + |κ|2 , 2 2 (y)

T 3 where ωt = [ρt ,u(x) t ,ut ] ∈ R is the flow field variables (macroscopic properties) vector at time t and ωtd is the desired flow field variables vector at time t. The first term in the cost function measures the distance of ωt and ωtd and the second term is a regularization term with σ  0 (typically σ ∈ [10−4 ,10−2 ]), which leads to improved smoothness properties of the optimization problem for σ > 0 [49]. In addition, the external body forces applied to the fluid are regarded as the design variables. Thus the following expression is substituted into Eq. (2.3):

Vi = cs2 (ci · κ).

=

C. Sensitivity analysis for the present problem

In this subsection, using the considered test case, the procedure of analytical evaluation of the macroscopic and microscopic adjoint variables and the cost function gradients are presented.

(5.1)

i = 5,6,7,8.

∂ ∂ f0 (xi,j ,t + 1) + f1 (xi,j ,t + 1) + · · · ∂ρt ∂ρt ∂ ∂ + f7 (xi,j ,t + 1) + f8 (xi,j ,t + 1), (5.6) ∂ρt ∂ρt

where xi.j ≡ (xi ,yj ) is the spatial discretized form of the location vector x for the lattice node (i, j ) th. By expanding of the LB equation (2.3) for each lattice node (i, j ) and with considering the equilibrium distribution function (2.4), all the distribution functions at time t + 1 can be rewritten in terms of the flow field variables at time t. Then by differentiating both sides of the obtained equation with respect to ρt , ∂fi,t /∂ρt are obtained. For example, ∂f8,t /∂ρt can be evaluated as ∂ eq ∂f8 (xi,j ,t + 1) = β f (xi−1,j +1 ,t) ∂ρt ∂ρt 8 = βw8 α8 (xi−1,j +1 ,t)|i=2,...,k

j =1,...,l−1 ,

∂f8 ∂ eq (x1,j ,t + 1) = β f (xk,j +1 ,t) ∂ρt ∂ρt 8 = βw8 α8 (xk,j +1 ,t)|j =1,...,l−1 ,

(5.7)

∂f8 ∂ eq (xi,l ,t + 1) = β f (xi,l ,t) ∂ρt ∂ρt 6 = βw6 α6 (xi,l ,t)|i=1,...,k ,

(5.5)

In the above equation, κ = [κ (x) ,κ (y) ]T ∈ R2 is the 2D body forces (design variables) vector. It can be shown that the existence of the force term (5.5) in the LB equation results in the recovery of the NS equations in the presence of the body force κ [50].

i = 1,2,3,4,

∂Itx /∂ωt , ∂ t+1 /∂ωt , ∂Itx /∂κ, and ∂ t /∂κ, which are obtained based on analytically derived expressions. Differentiating the cost function (5.4) with respect to the flow field or design variables, one can easily determine the sensitivity of Itx with respect to ωt or κ. For evaluating the derivative ∂ t+1 /∂ωt , first, it is proper to write it in a matrix form. For example, consider the first array of this matrix, i.e., ∂ (1) t+1 /∂ρt . According to the mass conservative law, we have 8

∂ (1) ∂ρt+1 ∂  t+1 = = fi ∂ρt ∂ρt ∂ρt i=0

√ Furthermore, the speed of sound is c/ 3 for the D2Q9 model. B. Cost function and design variables

i = 0,

where k and l are the number of lattice nodes along the horizontal and the vertical axis, respectively. To obtain the derivative ∂ t /∂κ using the LB method, it is appropriate to write it in a matrix form. Now, consider the (x) first array of the above matrix (∂ (1) t /∂κ ). According to the mass conservative, we have 8

 ∂ (1) ∂ρt ∂ t = = fi ∂κ (x) ∂κ (x) ∂κ (x) i=0 =

1. Sensitivity analysis in MADA approach

In MADA approach, solving the sensitivity equation (4.9) needs the evaluation of four partial derivatives: 013303-5

∂ ∂ f0 (xi,j ,t) + (x) f1 (xi,j ,t) + · · · ∂κ (x) ∂κ ∂ ∂ (5.8) + (x) f7 (xi,j ,t) + (x) f8 (xi,j ,t). ∂κ ∂κ

MOHAMAD HAMED HEKMAT AND MASOUD MIRZAEI

PHYSICAL REVIEW E 91, 013303 (2015)

Taking the derivative from both sides of Eq. (2.3) with respect to κ (x) , the derivatives of the distribution functions fi,t with respect to κ (x) will be available. Note that ∂ ∂ eq eq f0 (xi,j ,t − 1) = f (xi,j ,t − 1) = · · · (x) ∂κ ∂κ (x) 1 ∂ eq f (xi,j ,t − 1) = 0. (5.9) = ∂κ (x) 8 (x)

For example, ∂f8,t /∂κ is evaluated as     ∂f8 ∂f8  (xi,j ,t) = (x1,j ,t)  (x) (x) ∂κ ∂κ i=2,...,k j =1,...,l−1 j =1,...,l−1   ∂f8 = − (x) (xi,l ,t) ∂κ i=1,...,k = cs2 .

eq

The derivative ∂Ft /∂Ft can be evaluated using the chain rule. It should be noted that eq

∂fi,t ∂ρt

eq

= wi αi,t ,

∂ItxT ∂I T ∂ωt = tx , ∂Ft ∂ωt ∂Ft

∂ρt = [1,1,1,1,1,1,1,1,1]T , ∂Ft

(5.12)

∂u(x) −u(x) t t = [1,1,1,1,1,1,1,1,1]T ∂Ft ρt 1 + [0,1,0,−1,0,1,−1,−1,1]T , ρt ∂Ft

,

(y) ∂ut

= wi ρt

[1,1,1,1,1,1,1,1,1]T

ρt 1 + [0,0,1,0,−1,1,1,−1,−1]T . ρt

(5.14)

∂αi,t (y)

∂ut

.

(y)

ut

can be

i = 0,2,4, i = 1, i = 3, i = 5,8, i = 6,7, (5.18)

∂αi,t (y)

∂ut

⎧ 0, ⎪ ⎪ (y) ⎪ ⎪ ⎨3 + 9ut , (y) (y) 9ut , = −3ut + (−3) +   ⎪ i (x) ⎪ ⎪3 + 9 u(y) , t  − (−1) ut ⎪  ⎩ (y) , (−3) + 9 ut − (−1)i u(x) t

i = 0,1,3, i = 2, i = 4, i = 5,6, i = 7,8. (5.19)

D. Numerical solution of the MADA and MIDA equations

The solving procedure of the LB and MADA and MIDA equations are slightly different. In the LB method, the flow field variables vector is evaluated by solving the discrete LB equation at the time interval [1, tf ] as forward in time. While, in the adjoint method, the adjoint variables vector is evaluated by solving the MADA or MIDA equation at the same interval as backward in time. Figure 2 shows the evaluating procedure of the adjoint variables as backward in time using the evaluated flow field variables as forward in time.

(5.13)

Taking the derivative of both sides of Eq. (2.3) with respect to Ft , ∂ t+1 /∂Ft is achieved: ∂ t+1 ∂Ft

∂u(x) t

The derivative of αi,t with respect to u(x) t and computed as follows: ⎧ 0, ⎪ ⎪ (x) ⎪ ⎪ ⎨3 + 9ut , ∂αi,t (x) 9u(x) = −3ut + (−3) + t , (x)  ⎪ (y)  (x) ∂ut ⎪ ⎪ 3 + 9 u − (−1)i ut , t ⎪  (x) ⎩ (y)  (−3) + 9 ut − (−1)i ut ,

(5.11)

where ∂Itx /∂ωt can be determined as before. The derivatives (y) of ρt , u(x) t , and ut with respect to Ft are obtained using the conservative laws (2.5):

=

= wi ρt

eq

∂fi,t

(5.10)

To compute the microscopic adjoint variables vector using the MIDA equation, at first, the derivatives ∂Itx /∂Ft and ∂ t+1 /∂Ft should be determined. Since the flow field variables are functions of the distribution functions at time t, thus the derivative ∂Itx /∂Ft can be evaluated, using the chain rule, as follows:

(y) −ut

∂u(x) t

∂αi,t

(5.17)

2. Sensitivity analysis in MIDA approach

(y) ∂ut

∂fi,t

E. Optimization algorithm

After calculation of the gradient vector, the values of the design variables can be changed using an optimization algorithm. The steepest descent algorithm was adapted to treat the design variables towards the optimum values. In that algorithm, the design variables vector κ at the nth optimization cycle can be updated as κ n+1 − κ n = −c ∇κ I,

 ∂  eq (1 − β)Ft + β Ft + V = ∂Ft

(5.20)

Solution of the Conservation and the LB Equations

eq

where ⎡ ∂f

0,t

∂f

⎢ 0,t ∂Ft ⎢ . = ⎢ .. ⎣ ∂Ft ∂f8,t ∂f0,t

··· ..

.

···

⎤

⎡

(5.15)

1

0

⎥ ⎢0 .. ⎥ ⎢ = ⎢. . ⎥ ⎦ ⎢ ⎣ .. ∂f

···

1 .. .

···

∂f8,t

0

···

∂f0,t ∂f8,t

8,t

0

..

0

1

2

tf

1

tf

1

2

tf 1

tf

Initial Conditions

.

⎤ 0 0⎥ ⎥ ⎥ .. ⎥ . (5.16) .⎦ 1

Adjoint Terminal Condition

∂F ∂Ft +β t , = (1 − β) ∂Ft ∂Ft

Solution of the MADA/MIDA Equation

FIG. 2. (Color online) Flow field variables evaluated and collected as forward in time and used later for evaluating the macroscopic adjoint variables as backward in time.

013303-6

EXTRACTION OF MACROSCOPIC AND MICROSCOPIC . . .

PHYSICAL REVIEW E 91, 013303 (2015)

where c is the step length. Using a linear search method, the step length is chosen such that the maximum reduction in the cost function is achieved [51]. F. LB-based vs NS-based adjoint approach

The investigation of the derivation procedure of the MADA and MIDA equations shows that the simplicity of the LB equation, as an alternative to the NS equations, can facilitate the process of extracting the discrete adjoint equation and so the implementation cost of the method is lower. Furthermore, comparing the discrete LB equation with the MADA and MIDA one, it is observed that, as was expected, theses equations are quite similar to the discrete LB equation and they have the inherent aspects of the original LB equation (e.g., linearity, simplicity, parallelizability, etc.). Therefore, it is expected that the computational cost for solving this equation is nearly equal to that of the LB equation [23]. G. Extension of the proposed approaches to more realistic situations

The proposed LB-based adjoint methods are not limited to the type of the working fluids or flow fields. The present approaches are only optimization techniques that their dependence on the considered flow field is only summarized in extraction of the adjoint equations and computation of the adjoint variables. So, if every flow field is also selected and the goal is to minimize a cost function, the present methods can be easily performed by extracting the regarded adjoint equation (according to the similar mathematical procedure given in Sec. IV) and evaluating the adjoint variables and gradients (according to the comparable process provided in Sec. V C). Table I is a representation of the optimization procedure using the proposed adjoint approach based on the LB method. VI. RESULTS AND DISCUSSION

In this section, the results of the MADA and MIDA methods are presented. Various test cases have been used to implement the two proposed approaches. Because the similarity of the main results, only an elementary test case, described earlier in Sec. V, is presented. Here, all simulations of the flow field were performed using a 50 × 40 lattice and a fluid with kinematic viscosity of 1 (giving λ = 3.5) and initial density of 1. It should be noted that all of the quantities are in lattice units. A. Verification of the flow field and the adjoint solutions

Since the flow field solution performs a major role in the calculation of the cost function, the accuracy of the flow solution should be validated. The validation of the existing LB solver is based on a comparison with the numerical solutions of the NS equations. Figure 3 gives a comparison of the x-velocity profile for the present test case obtained from the LB solution and that of the NS at various times. It can be seen that the present results are in good agreement with the NS solutions. Moreover, before beginning the discussion, it is necessary that the adjoint gradients are validated. In order to verify the accuracy of the gradients calculated by the proposed adjoint approaches, they are compared with gradients obtained from

TABLE I. Optimization process using the proposed adjoint approaches based on the LB method. Step 1: Define optimization problem Select the cost function and design variables vector Select the state variables vector (the governing equation) Step 2: Extract the main equations Extract the MADA or MIDA equation Extract the adjoint terminal and boundary conditions Extract the cost function gradient vector Extract the partial derivatives of: the governing equation the cost function with respect to the state and design variables vector based on the LB method Until reach to convergence of optimization algorithm Step 3: Solve the LB equation as forward in time Set the initial conditions Until reach to terminal time: Compute the equilibrium distribution function Perform collision Perform streaming Apply the boundary conditions Compute the new distribution functions Compute the new state variables Store the new state variables Step 4: Solve the MADA or MIDA equation as backward in time Set the adjoint terminal condition Compute the partial derivatives of the cost function with respect to the state variables at terminal time Apply the adjoint terminal and boundary conditions Compute the macroscopic-microscopic adjoint variables at terminal time Until reach to initial time Compute the partial derivatives of the cost function with respect to the state variables Compute the partial derivatives of the governing equations with respect to the state variables Solve the MADA or MIDA equation Store the new adjoint variables Step 5: Evaluate gradient vector Compute the partial derivatives of the cost function with respect to the design variables Compute the partial derivatives of the governing equations with respect to the design variables Step 6: Modify design variables

the central finite difference (CFD) method. Here, the fluid has an initial parabolic velocity profile with maximum value of 0.1 along the horizontal axis together with a vertical velocity equal to zero and is driven by a constant volume force κ throughout time t = 1, . . . ,100. Also, we consider the states of the flow driven by body force κ d = [0.01, 0.002] as the desired ones. To select an appropriate step size (ε) and compute the accurate gradients using the CFD method, several step sizes were examined. Table II shows the relative error of the gradients calculated by the MADA method with respect to those of the CFD method for various step sizes, at κ = [−0.02,−0.002]. The results illustrate the growth of subtractive cancellation errors for ε  10−10 and the truncation error for ε  10−5 . Also, a good conformity can be seen between the results of

013303-7

MOHAMAD HAMED HEKMAT AND MASOUD MIRZAEI 40

Channel Width

t=50 t=55

PHYSICAL REVIEW E 91, 013303 (2015) Present Solver NS Solver

t=45

30 20 10

t=15 t=20 t=25 t=30 t=35 t=40

t=70 t=65 t=60

-0.015

-0.01

-0.005 x-Velocity

0

0.005

(a)

Channel Width

40 30 20 10

Present Solver NS Solver t=30 t=25 t=20 t=15 t=10 t=5

t=60 t=65

t=40 t=35

0

Present Solver NS Solver

0.01

0.02

0.03 0.04 x-Velocity

t=250 t=300 t=350 t=400

t=55 t=50 t=45

0.05

t=50 t=100 t=150 t=200

0.06

0

0.02

0.04 x-Velocity

(b)

0.06

0.08

(c)

FIG. 3. (Color online) Comparison of the x-velocity profile obtained from the present LB solution and the NS results for an unsteady laminar (x) = −0.0018; Poiseuille flow inside of a 2D periodic channel, at various times. (a) With initial velocity u(x) 0 = 0.01 and external body force κ (x) (x) (b) without initial velocity and with external body force κ = 0.01; (c) with initial velocity u0 = 0.1 and without external body force.

these two methods for 10−9  ε  10−7 . In Fig. 4, we show that the proposed methods are able to evaluate exactly the gradients of the cost function.

t = 1 to t = 100. In this case, the flow states driven by κ d = [−0.025,0] are regarded as the desired ones. Figure 5 shows the macroscopic adjoint states in the width of the channel obtained from the MADA method. The figures on the left column show the evolution of the macroscopic

B. Macroscopic adjoint properties vs adjoint distribution functions

nd

2 Design Variable

10 10−4 10−5 10−6 10−7 10−8 10−9 10−10 10−11 10−12

−2

3.63 × 10 1.09 × 10−3 4.33 × 10−5 2.71 × 10−7 8.94 × 10−8 8.35 × 10−9 8.02 × 10−9 6.22 × 10−8 9.35 × 10−7 7.42 × 10−6

G(y) relative error −1

1.73 × 10 3.37 × 10−2 6.29 × 10−5 8.78 × 10−7 1.00 × 10−9 2.29 × 10−10 9.42 × 10−10 9.61 × 10−9 1.20 × 10−7 6.58 × 10−6

0.002

0.004

Present MADA Approach Present MIDA Approach CFD Method

-0.04 -0.05

0 -0.06 -0.07

-5

-0.08 -0.09

-10

nd

−3

G(x) relative error

0

st

ε

-0.002

2 Component of Gradient Vector

TABLE II. Relative error of the gradients calculated by the MADA method (|(GCFD − GMADA )/GMADA |) at κ = [−0.02,−0.002]. Optimal κ d = [0.01,0.002].

-0.004 5

1 Component of Gradient Vector

The MADA and MIDA methods were implemented for the problem described in Sec. V. Here, the fluid has an initial parabolic x-velocity profile with maximum value of 0.01 and a zero y velocity and is affected by κ = [−0.0018,0] from

-0.1 -0.04

-0.02

0

0.02

0.04

0.06

st

1 Design Variable

FIG. 4. (Color online) Comparison of the first component of the gradient vector at κ (y) = 0.003 and the second component at κ (x) = 0.03 calculated by the MADA and MIDA methods with those of the CFD approach. Optimal κ d = [0.01,0.002]. 013303-8

EXTRACTION OF MACROSCOPIC AND MICROSCOPIC . . .

PHYSICAL REVIEW E 91, 013303 (2015)

Channel Width

40 t=95

30

t=98

Final time (t=100) t=99 t=1 t=10

20 t=90

10

t=80 t=70

0.99

Channel Width

40

0.995

1 Density

1.005

Final time (t=100)

1.01

t=30

-0.5 -0.4 -0.3 -0.2 -0.1 Macroscopic Adjoint Variable 1

Initial time (t=0)

0

0.02

Final time (t=100) t=99 t=97

t=10 t=5 t=1

30

0

20 t=93 t=85

10 t=60

-0.02

-0.01 x-Velocity

0

0.01

0

1 2 3 4 Macroscopic Adjoint Variable 2

5

6

Channel Width

40 30 t=90 t=85 t=60

20

Final time (t=100) t=99 t=95

10 t=15

-1.0x10-08

-5.0x10-09

0.0x10+00 y-Velocity

5.0x10-09

1.0x10-08

-0.02

-0.01 0 0.01 Macroscopic Adjoint Variable 3

0.02

FIG. 5. (Color online) Evolution of the flow field variables (macroscopic properties) obtained from LB method and the corresponding macroscopic adjoint variables obtained from the MADA approach at κ = [−0.0018,0]. Optimal κ d = [−0.025,0]. (y)

properties ρt , u(x) from t = 1, . . . ,100. Computed t , and ut through the MADA equation (4.8), the macroscopic adjoint states ψt(1) , ψt(2) and ψt(3) are shown on the left column of Fig. 5. It is interesting to observe the similar trends of the macroscopic adjoint variables as backward in time. As can be seen, the evolution of the first adjoint variable has two extension-retraction trends. From t = 100 to t = 1, ψt(1) grows along the negative direction until it reaches to a peak and retracts back to zero. Then, it extends slightly along the positive direction and reduces again. But, the evolution of the second and third adjoint variables has one extension-retraction behavior. ψt(2) extends along the positive direction until it reaches a peak and then reduces to zero. ψt(3) grows along the negative and positive direction until it reaches a peak and retracts back to zero. By looking at the construction of the MADA equation, ψt(1) can be viewed as a combination of the (y) (2) changes rate of Itx , ρt+1 , u(x) t+1 , and ut+1 with respect to ρt . ψt can be regarded as a combination of the changes rate of Itx , (y) (x) (3) ρt+ 1 , u(x) t+1 , and ut+1 with respect to ut , and similarly, ψt , (y) with respect to ut .

The obtained results from the adjoint sensitivity analysis demonstrate that, for all the times, the sensitivities of Itx and (y) ut+1 with respect to ρt are equal to zero and that of ρt+1 with respect to ρt is positive and constant but, the change rate of u(x) t+1 with respect to ρt decreases smoothly, as a parabola, from positive values (for the initial times) to negative values (for the intermediate and final times). So, according to the MADA equation, as backward in time, the growth in the evolution of ψt(1) at the final times is due to the increasing positive and high values of ψt(2) and, the increasing negative sensitivity of u(x) t+1 to ρt , along the positive direction. In addition, the extension at the initial times is due to the reducing positive and small values of ψt(2) and, the increasing positive sensitivity of u(x) t+1 to ρt , along the positive direction. Similarly, the evolution trend of ψt(2) can be explained with noting that the sensitivity of Itx to u(x) t reduces as a parabola from positive values (for (y) t = 100) to zero (for t = 1) smoothly, while ρt+1 and ut+1 have no sensitivity and u(x) t+1 has a positive and constant sensitivity . Likewise, the variations rate of Itx and with respect to u(x) t (y) ρt+1 are equal to zero and that of ut+1 is positive and constant,

013303-9

MOHAMAD HAMED HEKMAT AND MASOUD MIRZAEI

PHYSICAL REVIEW E 91, 013303 (2015)

(y)

with respect to ut but, the change rate of u(x) t+1 with respect to (y) ut is fully nonlinear such that, as backward in time, it extends on the upper half of the channel along the negative direction and on the lower half of one along the positive direction. So, the observed behavior in the evolution of ψt(3) is only affected (y) from the nonlinear sensitivity u(x) t+1 to ut . Figure 6 shows the microscopic adjoint states in the width of the channel obtained from the MIDA approach. The figures on the left column show the evolution of the microscopic

Channel Width

40 Final time (t=100)

properties f0 to f8 . Computed through the MIDA equation (4.8), the microscopic adjoint states ψ0 to ψ8 are shown on the left column of Fig. 6. In the D2Q9 model, the lattice directions and the corresponding velocity vectors have been considered such that they are symmetric with respect to the x and y axes. For example, the lattice directions 1 and 3 and the corresponding velocities have symmetry about the y axis, and similarly, 2 and 4, about the x axis, 5 and 8 about the x axis, 6 and 7 about the x axis, 5 and 6 about the y axis, and finally

Initial time (t=0)

Final time (t=100) t=99 t=97

t=10 t=5 t=1

30 20

t=90 t=80

10 t=65

0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 Distribution Function 1

0

1 2 3 4 Adjoint Distribution Function 1

40

Final time (t=100) t=99

Channel Width

t=98 t=95

30

t=1 t=10

20 10

t=90 t=95 Final time (t=100)

Initial time (t=0)

40 Initial time

t=90 t=80

t=10

t=30

0.105 0.11 Distribution Function 2

t=70

0.11 0.115 0.111

-0.5

t=30

-0.4 -0.3 -0.2 -0.1 Adjoint Distribution Function 2

0

0.10

0.2 0.02

Final time (t=100) t=99 t=97

Final time (t=100)

(t=0)

Channel Width

5

30

t=10 t=5 t=1

20 t=90

10

t=80 t=75

0.1

Channel Width

40 30

Final time (t=100)

0.15 Distribution Function 3

t=30

t=95 t=90

-6

0.2

-5

-4 -3 -2 -1 Adjoint Distribution Function 3

0

1

Final time (t=100)

t=10 Initial time (t=0)

t=99 t=98 t=95

t=1 t=10

20 t=90

10

t=80 t=70 t=30

0.105 0.11 Distribution Function 4

0.11 0.115 0.111

-0.5

-0.4 -0.3 -0.2 -0.1 Adjoint Distribution Function 4

0

0.10

0.2 0.02

FIG. 6. (Color online) Evolution of the distribution functions (microscopic properties) obtained from the LB method and the corresponding microscopic adjoint variables obtained from the MIDA approach at κ = [−0.0018,0]. Optimal κ d = [−0.025,0]. 013303-10

EXTRACTION OF MACROSCOPIC AND MICROSCOPIC . . .

Channel Width

40

PHYSICAL REVIEW E 91, 013303 (2015)

t=1

Final time (t=100) t=99

t=10

30

t=97

t=5 t=1

20

t=90

10 t=80

Final time (t=100)

t=75

0.01

0.015

0.02 0.025 0.03 Distribution Function 5

0.035

Channel Width

40

0.04

0

5

Final time (t=100) t=99 t=97

Final time (t=100)

30

t=10 t=5 t=1

20 10

t=90 t=80 t=70

t=1

0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 -7 Distribution Function 6 40

-6

-5 -4 -3 -2 -1 Adjoint Distribution Function 6

0

1

Final time (t=100) Final time (t=100)

Channel Width

1 2 3 4 Adjoint Distribution Function 5

t=99 t=97

30

t=10 t=5 t=1

20 t=90

10

t=1

t=80 t=70

0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 -7 Distribution Function 7

-6

40

0

1

Final time (t=100)

Final time (t=100)

Channel Width

-5 -4 -3 -2 -1 Adjoint Distribution Function 7

t=1

t=99

t=10

30

t=97

t=5 t=1

20

t=90 t=80

10

t=75

0.01

0.015

0.02 0.025 0.03 Distribution Function 8

0.035

0

0.04

1 2 3 4 Adjoint Distribution Function 8

5

FIG. 6. (Color online) (Continued.)

7 and 8 about the y axis. Therefore, according to the above description, we expect that the distribution functions may also have been such comparable similarities. Taking a meticulous look at Fig. 6, similar correlations between the distribution functions can be seen clearly. For example, from a qualitative point of view, the distribution functions f1 , f5 , and f7 are, actually, the rotated distribution functions f3 , f6 , and f8 about the y axis, respectively. Moreover, the rotation of f4 , f8 , and f7 about the x axis leads to the evolution of f2 , f5 , and f6 , both qualitatively and quantitatively.

Comparing the evolution of f0 to f8 to that of ψ0 to ψ8 in Fig. 6, one can see an interesting relationship between the results, as expected. According to Fig. 6, it should be noted that all the similar correlations, observed earlier, for the distribution functions f0 to f8 are observed for the corresponding adjoint distribution functions ψ0 to ψ8 . In addition to the points mentioned above, the extension-retraction behaviors, similar to those previously observed and described about the macroscopic adjoint properties, are well followed for the adjoint distribution functions.

013303-11

MOHAMAD HAMED HEKMAT AND MASOUD MIRZAEI 40

(

)/(

40

) -MIDA -MADA

t=55

20 t=45 10

t=35 t=25 t=99 Final time (t=100)

0

0.2

t=90 t=10 t=5 t=82 t=1 t=75 t=95 t=98

0.4 0.6 0.8 1 Normolized Adjoint Density

-0.02

0

FIG. 7. (Color online) Comparison of the normalized adjoint density obtained from the MADA and MIDA methods at κ = [−0.0018,0]. Optimal κ d = [−0.025,0].

Another fascinating issue is the relationships between the macroscopic adjoint properties and the adjoint distribution functions. Figure 7 shows the comparison of the macroscopic adjoint property ψt(1) normalized by ψt (1) max with the corresponding microscopic variable ψi normalized by (ψi )max , for t = 1, . . . ,100. It can be seen that the results from the both MADA and MIDA methods are entirely the same. For that reason, this property can be referred to as the normalized adjoint density. Similarly, Figs. 8 and 9 compare the normalized macro(3) (3) scopic properties ψt(2) /ψt (2) max and ψt /ψt max with the analogous microscopic variables (cix ψi )/(cix ψi )max and (ciy ψi )/(ciy ψi )max , obtained from the MADA and MIDA methods, respectively. According to these figures, an excellent conformity can be seen between the results of these two methods. Subsequently, these normalized properties can be seen as the normalized adjoint momentums. The present observations and results are promising that the macroscopic adjoint properties can likely be evaluated throughout the computational domain via the adjoint distribution functions. As a supplementary argument, an important and challenging question in the field of the adjoint method which can always be raised is what are the physical concepts of the adjoint variables? For example, by integrating from the flow field variables, the conservative laws, mass, momentum, and energy conservation are extracted. Here, a questionable issue may be brought up whether such laws exist for the adjoint variables. According to the MIDA equation, using the MADA equation and the conservative laws (2.5), and also considering cix ciy = cix = ciy = 0, we can generally derive the

Channel Width

40

Final time (t=100) t=99 t=97 t=93

30 20 t=30

t=20 t=10 10 t=1

(2)

( cix

0

( ciy i)/( ciy i)Max-MIDA (3)

30

Channel Width

Channel Width

/

PHYSICAL REVIEW E 91, 013303 (2015)

(2)

/ Max-MADA i)/( cix i)Max-MIDA

0.2 0.4 0.6 0.8 x-Normolized Adjoint Momentum

1

FIG. 8. (Color online) Comparison of the x-normalized adjoint momentum obtained from the MADA and MIDA methods at κ = [−0.0018,0]. Optimal κ d = [−0.025,0].

30 20

t=90 t=80 t=60 t=50

/

(3) Max

Final time (t=100) t=99 t=95

t=40

10

-MADA

t=20

-1

-0.5 0 0.5 y-Normolized Adjoint Momentum

1

FIG. 9. (Color online) Comparison of the y-normalized adjoint momentums obtained from the MADA and MIDA methods at κ = [−0.0018,0]. Optimal κ d = [−0.025,0].

adjoint mass and momentum conservative laws as follows:  ψi (x,t), ψ (1) = 1 i

ψ (2) = 2



cix ψi (x,t),

(6.1)

i

ψ (3) = 3



ciy ψi (x,t),

i

where 1 1 1

1 =  , 2 =  2 , 3 =  2 . i1 i cix i ciy

(6.2)

For the D2Q9 model 1 = 1/9 and 2 = 3 = 1/6. In the above laws, the macroscopic adjoint variables ψ (1) , ψ (2) , and ψ (3) can be named the adjoint density, x-adjoint momentum, and y-adjoint momentum, respectively. Moreover, ψi can be referred to as the adjoint distribution function along ith direction of lattice lines. For a detailed derivation of the adjoint conservative laws, the reader is referred to the Appendix. Finally, in Table III, the obtained LB-based adjoint concepts are summarized and compared with the corresponding NS-based ones. VII. CONCLUSIONS

In the present research, the discrete adjoint approach using the LB method was developed for optimization problems of unsteady flow fields. This development led to the MADA and MIDA approaches. The derivation procedure depicted that the proposed approaches could be performed similarly for any optimization problem with the corresponding cost function and design variables vector. Moreover, this approach was not limited to flow fields and could be employed for both steady and unsteady flows. First, the general MADA and MIDA equations and the cost function gradient vector were derived mathematically using the discrete LB equation. Also, for an elementary case, the analytical evaluation of the macroscopic and microscopic adjoint variables and the cost function gradients was presented. Finally, to interpret the macroscopic and microscopic adjoint variables in lattice space, the numerical implementation of the two present approaches for the inverse design of a 2D unsteady laminar Poiseuille flow in a periodic channel with constant body forces was carried out.

013303-12

EXTRACTION OF MACROSCOPIC AND MICROSCOPIC . . .

PHYSICAL REVIEW E 91, 013303 (2015)

TABLE III. New adjoint concepts in lattice space. Lattice space

Adjoint lattice space ∂I T

eq

T ∂ t+1 ψtT = ∂Ftxt + ψt+1 ∂Ft Microscopic discrete adjoint equation

Ft+1 = (1 − β)Ft + β Ft + V Lattice Boltzmann equation f0 to fm−1 : Distribution functions (microscopic properties)  ρ = i fi (x,t) Mass conservative law  ρu(x) = i cix fi (x,t) x-momentum conservative law  ρu(y) = i ciy fi (x,t) y-momentum conservative law

ψ0 to ψm−1 : Adjoint distribution functions (microscopic adjoint properties)  ψ (1) = 1 i ψi (x,t) Adjoint mass conservative law ψ (2) = 2 i cix ψi (x,t) x-adjoint momentum conservative law  ψ (3) = 3 i ciy ψi (x,t) y-adjoint momentum conservative law ψ (1) : Adjoint density ψ (2) : x-adjoint momentum ψ (3) : y-adjoint momentum

ρ : Density ρu(x) : x momentum ρu(x) : y momentum

The results showed that there are similar trends in the evolution of the macroscopic adjoint variables. These similarities were rationalized via the construction of the MADA equation and using the adjoint sensitivity analysis. Moreover, all the similar correlations, observed for the distribution functions, were observed for the corresponding adjoint ones. Besides, the results indicated that the normalized adjoint density and momentums

obtained from the MADA and MIDA methods were the same for all of the time. These observations suggested that the macroscopic adjoint properties can be evaluated throughout the computational domain via the microscopic ones. Similar to the conservative laws in the LB method, we generally derived the adjoint conservative laws that state relations between the macroscopic and microscopic adjoint ones.

APPENDIX: DERIVATION OF THE ADJOINT CONSERVATIVE LAWS

Consider the MIDA equation and its corresponding microscopic vector t . For t = tf , the terminal microscopic adjoint condition can be extended using the chain rule, as follows: (y)

ψ i,tf =

(x) ∂Itf x ∂Itf x ∂ρtf ∂Itf x ∂utf ∂Itf x ∂utf = + (x) + (y) . ∂fi,tf ∂ρtf ∂fi,tf ∂utf ∂fi,tf ∂utf ∂fi,tf

(A1)

The derivatives of the macroscopic properties with respect to the distribution functions can be evaluated by the conservative laws (2.5). Moreover, according to the terminal adjoint condition (4.7), the derivatives of the cost function with respect to the macroscopic properties are equal to the macroscopic adjoint ones. So, Eq. (A1) can be rewritten as follows: ψ i,tf =

∂Itf x = ψt(1) + cix ψt(2) + ciy ψt(3) . f f f ∂fi,tf

(A2)

Finally, using the above equation and considering cix ciy = cix = ciy = 0, we get    (1) (2) (3) i ψi,tf i cix ψi,tf i ciy ψi,tf , ψtf =  2 , ψtf =  ψtf =  . 2 1 i i cix i ciy

(A3)

In addition, for t = tf −1 , the MIDA equation can be extended using the chain rule and considering the conservative laws (2.5), as follows: ⎤ ⎡ (y) (x) ∂Itf −1 x ∂ρtf −1 ∂Itf −1 x ∂utf −1 ∂Itf −1 x ∂utf −1 ⎦ + (x) + ψ i,tf −1 = ⎣ (y) ∂ρtf −1 ∂fi,tf −1 ∂ut −1 ∂fi,tf −1 ∂ut −1 ∂fi,tf −1 



f

f

(x) (y) ∂ j,tf ∂ρtf −1 ∂ j,tf ∂utf −1 ∂ j,tf ∂utf −1 + + ψj,tf (x) + ψj,tf (y) ψj,tf ∂ρtf −1 ∂fi,t −1 ∂ut −1 ∂fi,tf −1 ∂ut −1 ∂fi,tf −1 j f f f ⎡ ⎤ ∂Itf −1 x ∂Itf −1 x ∂Itf −1 x =⎣ + cix (x) + ciy (y) ⎦ ∂ρtf −1 ∂utf −1 ∂utf −1

 ∂ j,tf ∂ j,tf ∂ j,tf + + cix ψj,tf (x) + ciy ψj,tf (y) . ψj,tf ∂ρtf −1 ∂ut −1 ∂ut −1 j f

013303-13

f

(A4)

MOHAMAD HAMED HEKMAT AND MASOUD MIRZAEI

PHYSICAL REVIEW E 91, 013303 (2015)

Using the final expression in Eq. (A4), we have ⎡ ⎤

  ∂Itf −1 x  ∂Itf −1 x ∂Itf −1 x ∂ j,tf ∂ j,tf ∂ j,tf ⎣ ψ i,tf −1 = + cix (x) + ciy (y) ⎦ + + cix ψj,tf (x) + ciy ψj,tf (y) ψj,tf ∂ρtf −1 ∂ρtf −1 ∂utf −1 ∂utf −1 ∂utf −1 ∂utf −1 i i i j           ∂Itf −1 x ∂Itf −1 x ∂Itf −1 x = 1 + cix + ciy (y) ∂ρtf −1 ∂u(x) ∂utf −1 tf −1 i i i  

 

 

       ∂ i,tf ∂ i,tf ∂ i,tf + 1 ψi,tf cix ψi,tf (x) ciy ψi,tf (y) + + . (A5) ∂ρtf −1 ∂ut −1 ∂ut −1 i i i i i i f

f

Since cix and ciy are equal to zero, the above equation can be simplified as      ∂Itf −1 x  ∂ i,tf ψ i,tf −1 = 1 + ψi,tf . ∂ρtf −1 ∂ρtf −1 i i i

(A6)

Substituting Eq. (A2) into (A6) leads to the following equation:         ∂Itf −1 x ∂ ∂ ∂ ψ i,tf −1 = 1 + ψt(1)

i,tf + ψt(2) cix i,tf + ψt(3) ciy i,tf . f f f ∂ρ ∂ρ ∂ρ ∂ρtf −1 i t −1 t −1 t −1 f f f i i i i

(A7)

According to the MADA equation, the expression in the bracket on the right-hand side of Eq. (A7) is equal to the first macroscopic adjoint variable. Therefore,  ψi,tf −1 (1) ψtf −1 = i . (A8) i1 Similarly, using the final expression in Eq. (A4), we have ⎡ ⎤

   ∂Itf −1 x ∂I t x ∂I ∂ j,tf ∂ j,tf ∂ j,tf tf −1 x f −1 2 2 ⎣cix cix ψ i,tf −1 = + cix (x) + cix ciy (y) ⎦ + + cix ψj,tf (x) + cix ciy (y) cix ψj,tf ∂ρtf −1 ∂ρtf −1 ∂utf −1 ∂utf −1 ∂utf −1 ∂utf −1 i i i j           ∂Itf −1 x ∂Itf −1 x ∂Itf −1 x 2 = cix + cix + cix ciy (y) ∂ρtf −1 ∂u(x) ∂utf −1 tf −1 i i i  

 

 

       ∂ i,tf ∂ i,tf ∂ i,tf 2 + cix ψi,tf cix ψi,tf (x) cix ciy ψi,tf (y) + + . (A9) ∂ρtf −1 ∂ut −1 ∂ut −1 i i i i i i f

f

Since cix and cix ciy are equal to zero, the above equation can be simplified as      ∂Itf −1 x  ∂ i,tf 2 cix ψ i,tf −1 = cix + ψi,tf (x) . ∂u(x) ∂utf −1 tf −1 i i i Substituting Eq. (A2) into (A10) leads to the following equation:      ∂Itf −1 x ∂  ∂  ∂  (1) (2) (3) 2 cix ψ i,tf −1 = cix + ψtf

i,tf + ψtf cix i,tf + ψtf ciy i,tf . ∂ρtf −1 ∂u(x) ∂u(x) ∂u(x) tf −1 i tf −1 i tf −1 i i i

(A10)

(A11)

According to the MADA equation, the expression in the bracket on the right-hand side of Eq. (A11) is equal to the second macroscopic adjoint variable. Therefore, the above equation is finally rearranged as follows:  cix ψi,tf −1 (2) ψtf −1 = i  2 . (A12) i cix Moreover, the above procedure can be performed to achieve the following equation:  ciy ψi,tf −1 (3) ψtf −1 = i  2 . i ciy 013303-14

(A13)

EXTRACTION OF MACROSCOPIC AND MICROSCOPIC . . .

PHYSICAL REVIEW E 91, 013303 (2015)

Therefore, by performing the similar mathematical operations, it can be found that the following adjoint conservative laws are generally satisfied for all the times:    ψi (x,t), ψ (2) = 2 cix ψi (x,t), ψ (3) = 3 ciy ψi (x,t), (A14) ψ (1) = 1 i

i

i

where 1

1 =  , i1

2 = 

[1] A. Jameson, Sci. Comput. 3, 233 (1989). [2] S. Nadarajah and A. Jameson, AIAA paper 2000-0667, 1 (2000). [3] M. B. Giles and N. A. Pierce, Flow Turbul. Combust. 65, 393 (2000). [4] M. H. Hekmat, M. Mirzaei, and E. Izadpanah, J. Mech. Sci. Technol. 23, 2479 (2009). [5] S. Nadarajah, Ph.D. thesis, Stanford University, 2003. [6] O. Tonomura, M. Kano, and S. Hasebe, Comput. Aided Chem. Eng. 28, 37 (2010). [7] P. E. V. Jacques and R. P. Dwight, Comput. Fluids 39, 373 (2010). [8] J. E. Hicken and D. W. Zingg, AIAA J. 48, 400 (2010). [9] R. D. Joslin, M. D. Gunzburger, R. A. Nicolaides, G. Erlebacher, and M. Yousuff Hussaini, AIAA J. 35, 816 (1997). [10] A. McNamara, A. Treuille, Z. Popovi´c, and J. Stam, ACM Trans. Graph. 23, 449 (2004). [11] J. B. Freund, J. Sound Vib. 330, 4114 (2011). [12] D. G. Cacuci, M. I. Bujor, and I. M. Navon, Sensitivity Uncertainty Analysis: Applications to Large-scale Systems (CRC, Boca Raton, 2005). [13] B. F. Sanders and N. D. Katopodes, J. Eng. Mech. 126, 909 (2000). [14] L. M. Arriola and J. M. Hyman, Comput. Sci. Eng. 9, 10 (2007). [15] F. Fang, C. C. Pain, M. D. Piggott, G. J. Gorman, and A. J. H. Goddard, Int. J. Numer. Methods Fluids 47, 995 (2005). [16] F. Fang, M. D. Piggott, C. C. Pain, G. J. Gorman, and A. J. H. Goddard, Ocean Model. 15, 39 (2006). [17] D. Furbish, M. Y. Hussaini, F. X. Le Dimet, and P. Ngnepieba, Tellus A, Dyn. Meteorol. Oceanogr. 60, 979 (2008). [18] K. Mani, Ph.D. thesis, The University of Wyoming, 2009. [19] G. Pingen, A. Evgrafov, and K. Maute, Comput. Fluids 38, 910 (2009). [20] M. M. Tekitek, M. Bouzidi, F. Dubois, and P. Lallemand, Comput. Fluids 35, 805 (2006). [21] G. Pingen, A. Evgrafov, and K. Maute, Struct. Multidisc. Optim. J. 34, 507 (2007). [22] J. M. T. Ngnotchouye, M. Herty, S. Steffensen, and M. K. Banda, Comput. Appl. Math. 30, 399 (2011). [23] M. H. Hekmat and M. Mirzaei, Adv. Mech. Eng. 2014, 1 (2014). [24] D. Yu, R. Mei, L.-S. Luo, and W. Shyy, Prog. Aerosp. Sci. 39, 329 (2003). [25] Y. Li, Y.-L. Duan, Y. Guo, and R. Liu, Appl. Math. Mech. 30, 445 (2009).

1 , 2 i cix

3 = 

1 . 2 i ciy

(A15)

[26] M. Mirzaei and A. Poozesh, Phys. Rev. E 87, 063312 (2013). [27] K. N. Premnath, M. J. Pattison, and S. Banerjee, Phys. Rev. E 79, 026703 (2009). [28] M. A. A. Spaid and F. R. Phelan, Phys. Fluids 9, 2468 (1997). [29] D. R. J. Owen, C. R. Leonardi, and Y. T. Feng, Int. J. Numer. Meth. Eng. 87, 66 (2011). [30] N. I. Prasianakis and I. V. Karlin, Phys. Rev. E 76, 016702 (2007). [31] H. Chen, I. Goldhirsch, and S. A. Orszag, J. Sci. Comput. 34, 87 (2008). [32] S. Gabbanelli, G. Drazer, and J. Koplik, Phys. Rev. E 72, 046312 (2005). [33] G. Pontrelli, S. Ubertini, and S. Succi, J. Stat. Mech. (2009) P06005. [34] X. Shan and H. Chen, Phys. Rev. E 47, 1815 (1993). [35] T. Wang and J. Wang, Phys. Rev. E 71, 045301 (2005). [36] J. Zhang and X.-Y. Lu, Phys. Rev. E 80, 017302 (2009). [37] H. Aono, A. Gupta, D. Qi, and W. Shyy, AIAA paper 2010-4867, 1 (2010). [38] G. McNamara and G. Zanetti, Phys. Rev. Lett. 61, 2332 (1988). [39] U. Frisch, B. Hasslacher, and Y. Pomeau, Phys. Rev. Lett. 56, 1505 (1986). [40] U. Frisch, D. d’Humi`eres, B. Hasslacher, P. Lallemand, Y. Pomeau, and J.-P. Rivet, Complex Syst. 1, 649 (1987). [41] S. Wolfram, J. Stat. Phys. 45, 471 (1986). [42] W. Gladrow and A. Dieter, Lattice-Gas Cellular Automata and Lattice Boltzmann Models: An Introduction (Springer, Heidelberg, Germany, 2000). [43] X. He and L.-S. Luo, Phys. Rev. E 55, R6333(R) (1997). [44] X. He and L.-S. Luo, Phys. Rev. E 56, 6811 (1997). [45] J. D. Sterling and S. Chen, J. Comput. Phys. 123, 196 (1996). [46] P. L. Bhatnagar, E. P. Gross, and M. Krook, Phys. Rev. 94, 511 (1954). [47] M. B. Giles and N. A. Pierce, AIAA paper 97-1850, 182 (1997). [48] Y. H. Qian, D. D’Humi´eres, and P. Lallemand, Europhys. Lett. 17, 479 (1992). [49] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich, Lecture Notes of the Autumn School Modeling and Optimization with Partial Differential equations, Hamburg, Germany, 2005 (unpublished). [50] Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E 65, 046308 (2002). [51] J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed. (Springer, New York, 2006).

013303-15