Topology optimization using the lattice Boltzmann method for unsteady natural convection problems

This paper proposes a density-based topology optimization method for natural convection problems using the lattice Boltzmann method (LBM). As the LBM can be developed as a completely explicit scheme, its attractive features over the traditional ones, such as the finite element method, are (1) suitability for solving unsteady flow problems and (2) scalability for large-scale parallel computing. We develop an LBM code for solving unsteady natural convection problems and provide its sensitivity analysis based on the so-called adjoint lattice Boltzmann method. Notably, the adjoint equation is derived from the discrete particle velocity Boltzmann equation and can be solved similarly to the original LBM concerning unsteady natural convection problems. We first show that the proposed method can produce similar results to the previous work in a steady-state natural convection problem. We then demonstrate the efficacy of the proposed method through 2D numerical examples concerning unsteady natural convection. As a large-scale problem, we tackle a 3D unsteady natural convection problem on a parallel supercomputer. All the developed codes written in C++ are available at https://github.com/PANFACTORY/PANSLBM2.git.


Introduction
Natural convection is a fluid phenomenon driven by temperature-dependent buoyancy force. This phenomenon is widely known as a heat transfer mechanism that does not require cooling fans and plays a vital role in designing highperformance cooling devices such as heat sinks in electronic systems. Since natural convection is often beyond the extent of human intuition when considering complex systems, research on mathematical optimization concerning natural convection has been studied for more than a half-century. Elenbaas (1942) discussed a maximum diffusion problem of natural convection in vertical plates. Bahadur and Bar-Cohen (2005) summarized various performance indices of pin fin heat sinks. Yu et al. (2011) proposed a surrogate-based optimization method for designing radial heat sinks. Although these traditional methods-the so-called sizing optimization-can provide improvements compared with existing devices, the obtained results are strongly dependent on the initial guesses in general and, therefore, limited in terms of the degree of design freedom.
Topology optimization (Bendsøe and Kikuchi 1988) is one of the methodologies in structural optimization to generate novel solutions under a high degree of design freedom. One of the most attractive features of topology optimization is that an optimized configuration can be generated from a fixed design domain even if the initial guess is quite far from an optimum. So far, topology optimization has been applied to various structural optimization problems and demonstrated its applicability on multiphysics, multiscale, practical engineering applications, and so on (Bendsøe and Sigmund 2003). Borrvall and Petersson (2003) pioneered the methodology of topology optimization for fluid problems, in which dissipation energy minimization problems are formulated under the Stokes flow assumption. After that, this has been applied to optimization problems concerning various fluid regimes, such as Navier-Stokes laminar flow (Gersborg-Hansen et al. 2005;Olesen et al. 2006), forced convection (Yoon 2010;Matsumori et al. 2013), chemical transport systems (Okkels and Bruus 2007;Chen et al. 2019), turbulence (Kontoleontos et al. 2013;Dilgen et al. 2018). Additionally, machine learning-based frameworks have been recently proposed for solving complex flow optimization problems (Yaji et al. 2022;Hammond et al. 2022) As mentioned in a recent review paper , the application to natural convection problems has been limited in the research community of topology optimization. One of the representative difficulties is that natural convection requires dealing with the two-way coupling; fluid and temperature fields interfere, whereas forced convection is the one-way coupling. This fact implies that the numerical simulation of natural convection is relatively more challenging and requires higher computational cost under the same degree of freedom as a case of forced convection. Alexandersen et al. (2014) pioneered topology optimization for 2D natural convection problems. After that, more detailed investigation on parameters of the 2D natural convection problems (Li et al. 2022b), large-scale 3D problems (Alexandersen et al. 2016), simplified models (Asmussen et al. 2019;Pollini et al. 2020), multiobjective problems (Ramalingom et al. 2019), and a combination with an adaptive mesh refinement scheme (Li et al. 2022a) have been reported. However, it should be noted that most previous works focus on steady-state natural convection. At the same time, the unsteady state is more practical in considering the natural convection phenomena. Furthermore, it may have applicability for designing practical and interesting optimization problems such as rapid cooling system design. Coffin and Maute (2016) proposed topology optimization concerning unsteady natural convection and discussed the specific feature of the optimized design in the unsteady case comparing with the steady one. Besides, topology optimization of unsteady natural convection has been successfully applied to designing energy storage systems (Pizzolato et al. 2017;Vogel and Johnson 2019;Xie et al. 2019;Zhao et al. 2020). These previous works have been limited in 2D problems due to the massive computational cost for sensitivity analysis, as mentioned in Coffin and Maute (2016). Such a full-scale 3D topology optimization may be realized by using a large-scale parallel computer. However, since the previous works mainly use implicit schemes such as the finite element method (FEM), it is challenging to exploit its parallel performance compared to explicit schemes that do not require matrix operations. On the other hand, since an explicit scheme generally requires a lot of time steps compared to an implicit scheme, the memory requirement may be strict because the adjoint method for time-dependent problems requires preserving the state variable field for each time step. Hence, it cannot be concluded that an explicit scheme is more efficient than an implicit scheme, but we believe that a topology optimization approach based on an explicit scheme is a promising option for solving unsteady natural convection problems if a parallel computer could be used.
Among the many explicit schemes, the lattice Boltzmann method (LBM) (Chen and Doolen 1998;Aidun and Clausen 2010;Krüger et al. 2017) is one of the representative schemes of the computational fluid dynamics (CFD). The LBM is inspired by the analogy of molecular fluid mechanics and can compute fluid regimes from fluid particles' so-called velocity distribution functions governed by the Boltzmann equation. In general, the LBM does not require solving the Poisson equation that is the majority of the computational cost in the other standard CFD methods based on a segregated approach in terms of the fluid velocity and pressure. Hence, the LBM can be coded as a simple numerical algorithm that realizes high-performance parallel computing. Besides, many researchers have reported the applicability of the LBM for computing natural convection fields (Dixit and Babu 2006;Mohamad and Kuzmin 2010;Li et al. 2016). Pingen et al. (2007) have opened the door of topology optimization using the LBM and first indicated that the optimized configurations are similar to those of the case using the FEM in pressure drop minimization problems. After that, their method has been applied to Tesla valve design (Pingen et al. 2008), scalar transport problems (Makhija et al. 2012) and so on. However, since these previous works use the traditional discrete adjoint sensitivity analysis, large-scale and asymmetric matrix needs to be treated due to the feature of velocity distribution functions. Yaji et al. (2014) proposed a topology optimization method based on the adjoint lattice Boltzmann method (ALBM) (Krause et al. 2013), in which not only the forward problem but the adjoint problem are solved using the LBM. The basic idea of the ALBM is that the adjoint equation is derived in the continuous fashion from the continuous Boltzmann equation. Then the adjoint equation is discretized using the LBM. Hence, the adjoint problem can be solved by a completely explicit scheme without any matrix operations. This method has been applied to pressure drop minimization problems (Yaji et al. 2014;Xie et al. 2021), convection problems (Yaji et al. 2016;Dugast et al. 2018;Luo et al. 2022), chemical transport problems (Dugast et al. 2020), and unsteady flow problems (Chen et al. 2017;Yaji et al. 2018;Nguyen et al. 2020). Besides, there is another way to formulate the adjoint problem using the LBM, in which the original discrete equation-lattice Boltzmann equation-is used for deriving the discrete adjoint equation (Tekitek et al. 2006). This method has also been applied to steady-state problems (Liu et al. 2014;Łaniewski-Wołłk and Rokicki 2016) and unsteady problems (Nørgaard et al. 2016(Nørgaard et al. , 2017. Yonekura and Kanno (2015) proposed a radical method based on a one-shot approach, in which solving the adjoint problem can be avoided under a low Reynolds number condition. Notably, these LBM-based topology optimization methods are suitable for large-scale parallel computing as the adjoint problem has the feature of the explicit scheme in the LBM. In addition, except Yonekura and Kanno (2015), these methods can be naturally applied to unsteady problems by a slight change from the steady-state case, since the LBM is originally a numerical scheme for solving time-dependent flows.
In this study, based on Yaji et al. (2014Yaji et al. ( , 2016, we propose a topology optimization method for unsteady natural convection problems. We derive the adjoint problem in a general formulation, in which the ALBM based on discrete velocity Boltzmann equations is employed for reflecting the LBM boundary conditions to the adjoint problem (Yaji et al. 2016). It should be noted that the ALBM in this study is a semi-continuous adjoint method concerning continuous time and space and discrete particle velocities. Then, we provide the concrete optimization algorithm as the numerical implementation. In numerical examples, we validate the developed code by solving a steady-state case (Alexandersen et al. 2014) as a benchmark. We then tackle two cases of unsteady natural convection problems. The objective function of the first example is a time-averaged temperature at a heated boundary, and that of the second is a weighted sum of time-averaged flow rate and the deviation for designing a natural convection-driven pump. As a large-scale problem, we tackle a 3D unsteady natural convection problem on a parallel supercomputer. Lastly, we summarize what we found in this study as the conclusion.

Topology optimization
Here, we briefly review the basic concept of topology optimization in a general expression. Let us consider a structural optimization problem to determine the boundary of a design domain Ω . The basic concept of topology optimization is the introduction of a fixed design domain D ⊂ ℝ d (d: spatial dimension) that includes the original domain, i.e., Ω ⊆ D . Then, a characteristic function Ω ∶ D → {0, 1} is introduced to replace the original optimization problem with a material distribution problem in D. Herein, the characteristic function is defined as: where x represents a position in D. Since the discrete nature of Ω makes it hard to solve such a 0-1 optimization problem directly, relaxation or regularization techniques are typically essential. One of the most popular methods is the so-called density approach (Bendsøe and Sigmund 2003), whose basic concept is to replace Ω with a continuous function, i.e., ∶ D → [0, 1] , where is termed a design variable field. Using , a general form of the topology optimization problem based on the density approach can be expressed as follows: where J is an objective functional, and G is a constraint functional. Besides, U represents the state variable field that is typically given as the solution of partial differential equations for physical fields. In this study, U corresponds to the fluid velocity, pressure and temperature in natural convection systems. After some discretization treatments, such as the FEM for each field, the optimization problem (2) can be solved using a gradient-based optimizer.

Governing equations
In this study, the state variable field U in (2) is given as the solution of an unsteady natural convection system. The time interval is defined as I = [t 0 , t 1 ] where t 0 and t 1 are the initial and final times, respectively. We consider an analysis domain O = D ∪ Ω non with D ∩ Ω non = � in which Ω non represents the non-design domain. We assume that the analysis domain is composed of the fluid domain Ω f and the solid domain Ω s . In addition, let us consider that the design variable field is defined as = 0 for Ω s and = 1 for Ω f , respectively. Note that the grayscale region, 0 < < 1 , is regarded as a porous medium based on the density approach.
Consider the following natural convection system based on the Boussinesq approximation for the fluid velocity u ∶ Ω f × I → ℝ d , pressure (divided by the fluid density), p ∶ Ω f × I → ℝ and temperature T ∶ O × I → ℝ in the fluid domain: (1) where is the kinematic viscosity, e is the unit vector of the gravitational direction, g is the gravitational acceleration, is the volumetric expansion coefficient, T ref is a reference temperature which depends on the analysis problem, and K f is the thermal diffusivity of the fluid. Additionally, the temperature T satisfies the following equation in the solid domain: where K s is the thermal diffusivity of the solid.
Since our concern is to determine the distribution of Ω f and Ω s in D, the governing equations in (3-6) need to be expanded so that the state variable fields belong to both domains using the design variable field . For this, we employ a Brinkman approach (Borrvall and Petersson 2003), in which u belongs to not only the fluid domain but the solid domain, by replacing (4) with Herein, is a design-dependent parameter expressed using , as follows: where the parameter ̄ has a sufficiently large positive number for expressing the solid domain in D, and q > 0 is a tuning parameter to control the convexity of . As the similar way in (8), the governing equations for T in (5) and (6) are combined for expressing the temperature field by a single equation in D, as follows: where the design-dependent thermal diffusivity K is defined as Herein, q K > 0 is a parameter for controlling the convexity of K . Note that these governing equations can be solved under appropriate boundary conditions.

Lattice Boltzmann method (LBM)
In this study, the lattice Boltzmann method (LBM) is used to obtain the macroscopic variable fields discussed in Sect. 2.2. The basic idea of LBM is that the fluid regime is represented as an aggregation of fictitious particles, which makes it possible to obtain macroscopic variable fields such as the fluid velocity, pressure and temperature from the moments of the velocity distribution functions. Based on the basic idea of the LBM (Krüger et al. 2017), we consider a modeled thermal fluid composed of identical particles whose velocities are restricted to a finite set of Q vectors, c 0 , c 1 , ⋯ , c Q−1 ∈ ℝ d , where the number of Q depends on the lattice model discussed later.
Let us now discuss how to solve the governing equations in (3, 7 and 9) by the LBM. We assume that the fluid pressure and velocity are given by the moments of the velocity distribution functions f i ∶ O × I → ℝ with i = 0, 1, ⋯ , Q − 1 , as follows: where is the fluid density. Similarly, the temperature is given by using another velocity distribution function Note that although (13) needs not to be defined using the same Q with (11 and 12), we use the same value in this study for brevity. From the Bhatnagar-Gross-Krook (BGK) model in the kinetic theory, f i and g i respectively satisfy the following discrete velocity Boltzmann equations: Here, both equations are nondimensionalized as with the previous work (Yaji et al. 2016), and Sh is the Strouhal number given by the ratio of reference values of the flow speed and the particle speed. Besides, f and g are dimensionless parameters of the same order as the Knudsen number. Note that the design-dependent terms in (14 and 15) are equal to zero in Ω non . In addition, f eq i and g eq i represent the local equilibrium distribution functions given by Sh The parameter w i in (14, 16 and (17) denotes a weight in each lattice grid and its value is given in each lattice model. In this study, two kinds of lattice models are used with the nine-velocity (D2Q9) model for 2D problems and the fifteen-velocity (D3Q15) model for 3D problems. Figure 1 shows the schematic of the lattice models and the exact values of the weights w i are given by for D2Q9 model, and for D3Q15 model.
Consider a discretized equations for (14 and 15), in terms of the position x and the time t using a lattice spacing Δx and a time step Δt . In the following, we assume Δx = Δt = 1 . Based on the discretization way in (Inamuro et al. 2002), we obtain the so-called lattice Boltzmann equations for f i and g i as follows: 8,9,10,11,12,13,14) (20) .
Besides, f = f ∕Δx and g = g ∕Δx are the dimensionless relaxation time.
From the asymptotic expansion theory in the LBM, the dynamic viscosity and the thermal diffusivity K in the LBM can be given as follows: Herein, the relaxation times f and g in (20 and 21) are determined so that (22 and 23) are satisfied. Besides, g is the function of from the definition of K in (10).

Initial and boundary conditions
In order to solve the lattice Boltzmann equations in (20 and 21), appropriate initial and boundary conditions are needed for the thermal-fluid problem.
As for the initial condition, we assume that the initial values of f i and g i are equal to those of f eq i and g eq i , as follows: On the other hand, as for the boundary condition, many equations are proposed for thermal-fluid problems. Let us consider how to determine the boundary condition in the D2Q9 model in Fig. 1a as an example, in which f i and g i ( i = 2, 5, 6 ) are unknown at the lower boundary in the y-direction. As for the no-slip condition, we apply the socalled bounce-back condition: In addition, as for the free-slip condition, we apply the mirror-reflection condition: Δx.
(24) As for the prescribed temperature and heat flux conditions, we use the following equations (Inamuro et al. 2002;Yoshino and Inamuro 2003): where T ′ is given by for the constant temperature condition and for the constant heat flux condition. Herein, T w is the prescribed temperature at the intended boundary, and u y and q y are the velocity and heat flux in the y-direction, respectively. Herein, the heat flux q can be calculated by Note that the boundary conditions for the other boundaries and the case of the D3Q15 model in Fig. 1b can be treated as the similar way discussed above.

Adjoint lattice Boltzmann method
Here, we discucss a strategy based on the adjoint lattice Boltzmann method (ALBM) (Krause et al. 2013) to derive design sensitivities for our topology optimization problems. The basic idea of the ALBM is that the adjoint equation is derived in the continuous fashion from the continuous Boltzmann equation and then is discretized using the LBM. Hence, the discretized adjoint equation can be solved in an explicit manner without any matrix operations. Although the original ALBM is supposed to use a fully continuous Boltzmann equation for the adjoint analysis, we use the discrete velocity Boltzmann equation for reflecting the practical LBM boundary conditions to the adjoint equation (Yaji et al. 2016). In the following, we discuss the derivation procedures for the adjoint equation based on the ALBM in the natural convection system described in Sects. 2.2-2.4.
Assuming that the objective functional J is given by the following multiple integrals with respect to the analysis domain O or its boundary O and time interval I: represent the general expressions of the functions to be integrated. Let j O and j O be functions of (f i , g i ) in this study. Following the ALBM procedure, the Lagrangian L given by J and the governing equations in (14) and (15) can be expressed as The sensitivity of L with respect to the design variable field is identical to that of J since the governing equations must be zero in the analysis domain. The functional derivative of L with respect to is given by where f i = ( f i ∕ ) and g i = ( g i ∕ ) . Besides, n represents the outward normal vector on the boundary, and m = ∑ Q−1 i=0 w i c ifi . Let us now define the macroscopic adjoint variable fields using f i and g i : Using (35)-(38), f eq i and g eq i in (34) that correspond to the local equilibrium distribution functions can be written as It should be noted that although the derivation of the Lagrangian L contains its partial derivative with respect to in general, we drop it for simplicity because J does not explicitly contain in this paper. In addition, boundary condition terms should be included to L, but it becomes complicated and depends on the lattice model. Hence, we also drop the boundary condition terms for simplicity.
To determine f i and g i as the adjoint variable field from (34), we consider the following equations: Note that the adjoint equations closely resemble the original discrete velocity Boltzmann equations in (14) and (15) and therefore can be solved using the LBM, as follows: Similar to treating the lattice Boltzmann equation, the initial and the boundary conditions are needed for solving (43 and 44). Due to the nature of the time-dependent optimization problem, the above adjoint equations need to be solved backward in time from the final time t 1 to the initial time t 0 in the original time interval I . These terminal conditions can be derived from the fifth and eighth terms of the right-hand side in (34), as follows: Besides, the boundary conditions can be determined from the second, fourth, sixth and nineth terms of the right-hand side in (34). Here, we notice that the velocity distribution functions for particles flowing outside on the prescribed boundary must be applied since the direction of the fictitious particles for the ALBM is opposite to those for the LBM. In the following, we will discuss briefly how to determine the boundary conditions for the 2D problem by using the D2Q9 model as an example. Assuming that the boundary is located at the lower side in the y-direction, which is the same manner to solve the lattice Boltzmann equation. For the case of the no-slip condition, the following equation can be derived by applying the bounce-back condition to (34) as For the case of free-slip condition, the following equation can be derived by considering the mirror-reflect condition as In addition, as for the constant temperature condition, the following equations can be derived by substituting (28) with (29) into (34) as t).
As for the constant heat flux condition, the following equations can be derived by substituting (28) with (30) into (34)  as Consequently, the sensitivity L � = J � can be obtained from Note that the sensitivity for a constraint, G ′ , can also be obtained from the same way by replacing J with G, where j O and g O are respectively replaced by j O and g O in (32). Besides, when setting a prescribed heat flux condition, the gradient of the thermal diffusivity K must be considered in the sensitivity equation. Hence, when the value of the prescribed heat flux on the boundary is not equal to zero, the following term must be added into (51).
where Γ q is the boundary at which the heat flux q n = q ⋅ n is applied.

Dimensionless parameters
The previous study (Alexandersen et al. 2014) reported that the balance between the convection-dominated heat transfer and the diffusion-based heat transfer would affect the optimized shape significantly in topology optimization for natural convection problems.
In this study, the Rayleigh number and the Prandtl number are introduced to describe the thermal-flow state. The definitions of these dimensionless parameters are given as follows: (49) g 4 =g 7 =g 8 = − 1 6 4g 2 +g 5 +g 6 − u x 1 + 3u y g 5 −g 6 . (50) for the Rayleigh number which describes the ratio between the buoyancy force and the thermal diffusivity, and for the Prandtl number which describes the relative spreading of viscous and thermal effects. In the above, ΔT is the reference temperature difference, L is the reference length scale, and is the kinematic viscosity of the fluid.

Optimization algorithm
The numerical procedure of our proposed method is described as follows: Step 1 The design variable field is initialized on the analysis domain discretized using square lattice mesh.
Step 2 The state variable fields , u , T, q are computed using the LBM.
Step 3 The values of the objective functional J and the constraint functional G are computed.
Step 4 The adjoint variable fields ̃ , j , T , q are computed using the ALBM.
Step 5 The design sensitivities J ′ and G ′ are computed from the state and adjoint variable fields.
Step 6 The design variable field is updated using the method of moving asymptotes (MMA) (Svanberg 1987).
Step 7 The procedure returns to the Step 2 of the iteration loop until the following criterion is met: where superscript k represents the number of iterations carried out during the process. In this study, the criterion is set to opt = 0.01 for all numerical examples in Sect. 4. As for the MMA, it is widely known that the algorithm is suited for nonlinear optimization problems with a large number of design variables and few constraints. In this study, the move limit is set to 0.2, and the otherwise, the default values are set in our calculation.

Approximate procedure for a steady-state problem
Although the main focus of this study is unsteady natural convection problems, we will solve a steady-state natural convection problem for the verification of the proposed framework in Sect. 4.1. It is known that the computational cost can be reduced in LBM-based topology optimization for steady-state problems (Yaji et al. 2014;Liu et al. 2014). The basic idea of this reduction technique is that the converged state and adjoint variable fields at each optimization step are used as an initial value for solving the time evolution equations in the next optimization step. As for the state variable field, we use the following convergence criteria for the steady-state condition: where u and q represent the criteria for the steady-state condition of fluid velocity and heat flux, respectively. Hereinafter, the state variable fields for the density, velocity, temperature and heat flux under the steady-state condition are denoted as * , u * , T * , q * , respectively. At the first optimization step ( k = 0 ), (24) and (25) are used as initial conditions, and the second step onwards, the converged values of * , u * , T * , q * at the former step are used as the initial conditions for the current step. In other words, we use the following equations as the initial conditions: As for the adjoint variable fields, * , u * , T * , q * are used for (41) and (42), and the adjoint variable fields are calculated until they converge. The convergence criteria for the steadystate condition of the adjoint variable fields are defined as follows: where j and q represent the criteria for the steady-state condition of j and q , respectively. Note that since the direction of the time evolution of the ALBM is opposite with the LBM, the convergence criteria in (60) and (61) are different from (56) and (57), in terms of the time step. Using the similar notation with the state variable fields, the adjoint variable fields under the steady-state condition are denoted by ̃ * , j * , T * , q * , respectively. Besides, at the second step onwards in the optimization process, ̃ * , ũ * , T * , q * calculated at the former step are used as the final time values that corresponds to the initial condition for the adjoint equations. Hence, we employ the following equations for f i and g i at t = t 1 : Finally, the sensitivities for the steady-state problem can be expressed as

Approximate procedure for a periodic problem
In Sect. 4.3, we will investigate an unsteady natural convection problem, in which the state variable fields are periodically steady-state. Ideally, the state variable fields must be recorded for all time steps, including the developing phase that does not achieve the periodically steady-state, and the adjoint equations are then solved using the results of the state variable fields. However, it results in prohibitive memory requirements in large-scale problems.
To suppress the amount of required memory, we divide I into two intervals, namely, the developing interval I 0 = [t 0 , t † ) and the periodically steady-state interval I 1 = [t † , t 1 ] . In I 0 , the state variable fields are stored at the final time step only, as our interest is periodically steadystate. On the other hand, all the state variable fields are stored for solving the adjoint equations in I 1 . In practice, we skip to calculate the state variable fields in I 0 and their initial conditions at (k + 1)-optimization step are defined as those of the final time at k-optimization step. That is, each initial condition for the state variable fields are given by The advantage of using (65) and (66) is that the state variable fields are expected to rapidly achieve periodically steady-state, as long as drastically topological changes are disallowed, in comparison with the case in which the state variable fields are initialized every optimization step. Note that I 0 with t † = t 1 is used for computing the state variable field at the initial optimization step ( k = 0 ). Additionally, the adjoint variable fields are computed in the same manner discussed above, and their initial conditions can be given by Herein, (67) and (68) are used in Sect. 4.3 and are given under the assumption that drastically topological changes are disallowed. In addition, we define t † as a reasonable value for each optimization problem by preliminary computation.

2D steady-state problem
We first conduct a verification check of the proposed method through comparison with results of the previous research on a 2D heat sink design problem concerning steady-state natural convection (Alexandersen et al. 2014). Fig. 2 shows the dimensions and boundary conditions for the analysis domain. Herein, the reference length is defined as the height of the analysis domain, namely, L = 160Δx . Since the proposed method is developed for unsteady flow problems, we suppose that the natural convection is driven by a heat source and finally converges to the steady-state flow. The objective functional is defined as the average temperature at the heated boundary Γ SC : Besides, we consider a volume constraint in terms of the solid region, given by where we set v s = 0.5 that corresponds to impose the 50% allowable solid volume as a constraint.
The heat flux at Γ SC is set to q 0 = 1.0 × 10 −2 . The Prandtl number Pr and the reference temperature T ref are set to 6.0 and 0, respectively. The thermal conductivity of the solid is set to K s = 10K f , where that of the fluid is given by K f = ∕Pr from (54). The kinematic viscosity is set to = 0.1 for all numerical examples in this study. The initial configuration is defined as a uniform distribution, = 0.5 . In addition, the parameters u , q , ũ , q relating the convergence criteria are all set to 1.0 × 10 −6 . It should be noted that the parameters mentioned above are different from the previous work (Alexandersen et al. 2014) since the LBM used in this study is difficult to realize the same , K f and K s settings due to the limitation of the BGK model. Other riched LBM models, such as the multi-relaxation time (MRT)-LBM (Li et al. 2016), may achieve those, but it is beyond the scope of this paper.
To avoid generating the checkerboard pattern, we employ the density filter (Bourdin 2001), in which the filter radius is set to 2.4Δx in this study. Furthermore, we use a Heaviside projection (Wang et al. 2011) to remove the grayscale generated by the density filter, where the threshold parameter is set to 0.5 and the sharpness parameter is doubled every 100 optimization steps from 1 to 16.
As mentioned in the previous work (Alexandersen et al. 2014), a branched tree-like structure is obtained as the optimized design under a small Ra setting; on the other hand, a simple structure composed of smaller and shorter branches is obtained under a larger Ra setting. This is because that conduction and convection are respectively dominant when low and high Ra settings. In fact, as shown in Fig. 4, it can be confirmed that the effect of natural convection in terms of the magnitude of the fluid velocity is stronger when setting larger Ra. Table 1 shows the objective functional values of the optimized designs under different Ra settings for a crosscheck. The result indicates that each optimized design is reasonable as they perform better than the others for its particular Ra settings. Consequently, it can be confirmed that the proposed framework can provide appropriate solutions under the investigated steady-state condition.

2D transient problem
As the second numerical example, we tackle a 2D heat sink design problem concerning transient natural convection, which is regarded as an expansion of the steady-state case discussed in Sect. 4.1. In Sect. 4.2, we demonstrate that the optimized designs considering the undeveloped phase of natural convection achieve different shapes in comparison with those of the steady-state case. We consider that the objective functional for the transient problem is defined as the time-averaged one of (69): where, the time interval is set to I = 0, t 1 , in which we investigate the effect of the optimized results with respect to t 1 . Otherwise, the dimensions and parameter settings are same with those of the steady-state problem in Sect. 4.1. Note that, for forced convection problems, Zeng et al. (2020) demonstrated a transient formulation is effective for instantaneous chip cooling problems rather than a steady-state formulation. Figure 5 shows the optimized configurations, and Fig. 6 shows the velocity magnitude and temperature for three cases of t 1 . Here, the Rayleigh number is set to Ra = 2.0 × 10 5 . It can be confirmed that the temperature field concentrically propagates along with the undeveloped velocity field when a small t 1 setting. On the other hand, when a large t 1 setting, the temperature field propagates to the upper region with the developed velocity field. Notably, although Ra = 2.0 × 10 5 in Fig. 5 is a convection-dominant setting for the steady-state case disucussed in Sect. 4.1, conduction-dominant structure that has many needle-like structure is obtained under a small t 1 setting. This is because the objective functional in this optimization problem is defined as the time integral in terms of 0, t 1 . Therefore, the undeveloped phase occupies the majority of the interval if a small t 1 setting. Table 2 shows the objective functional values of the optimized configuration under different t 1 settings for a  crosscheck. As can be confirmed in Table 2, each optimized configuration is reasonable as they perform better than the others for its particular t 1 settings. Figure 7 shows the optimized configurations for three cases of the Rayleigh number Ra = 5.0 × 10 4 , 1.0 × 10 5 and 2.0 × 10 5 , respectively. Here, t 1 is set to 1.5 × 10 5 for both cases. Similar to the case in the steady-state problem, the optimized configuration in the transient problem appears to have large convection cells when setting the larger Ra.
Let us discuss the dependency of the optimized configuration with respect to t 1 and Ra from the perspective of the dominant heat transfer mechanism for the transient problems. For this, we introduce the cumulative norm of flow velocity, U cn , given by We use U cn to evaluate the magnitude of convection in terms of the accumulation of the velocity norm at each optimized configuration. Figure 8 shows the relationship between t 1 and U cn for three Ra settings, including the optimized shapes at each t 1 setting. As discussed in Sect. 4.1, since convection plays a significant role in the heat transfer when larger t 1 and Ra settings that correspond to a larger U cn , it can be confirmed that the optimized shapes tend to be simpler. Interestingly, around U cn = 100 in Fig. 8, it can be confirmed that both optimized configurations look similar, as can be also seen in Figs. 5b and 7b. This result indicates that U cn can be a dominant factor to determine the optimized shape for any t 1 and Ra settings.
It should be noted that the optimized designs in Figs. 5c and 7c contain small grayscale islands. Their effect on performance is negligible since the relative error between the optimized design and its binarized one is 0.45%.
To further investigate the obtained result on a different perspective, we conduct a high-fidelity analysis using the bodyfitted mesh on the boundary of the optimized designs. Figure 9 shows comparison results on fluid velocity and temperature between the high-fidelity and the original analyses at the final time step in the case of t 1 = 1.0 × 10 5 . It can be observed that the distributions are similar to each other. On the other hand, the relative infinity norms based on each original result, in terms of fluid velocity and temperature, are 16.3% and 5.56%, respectively. Table 3 shows the crosscheck of the objective functional values of the high-fidelity analysis for Fig. 5. The  Table 2, is 6.02% between the high-fidelity and the original results.
Since the exact boundary is given in the body-fitted model, the heat transfer coefficient h can be precisely obtained from the high-fidelity analysis. As a useful dimensionless number using h, let us now check the Biot number: where k is the thermal conductivity of the solid, and h is calculated as its average in terms of time between 0 and 1.5 × 10 5 . According to the previous work on the optimal shape of a fin in natural convection (Alexandersen and Sigmund 2021), the length of the optimal fin at smaller Bi is longer. The Biot numbers of the optimized design for Fig. 5 are Bi = 0.2903, 0.3351, 0.7550 , respectively. Hence, it can be confirmed that our results also have the relationship between Bi and the fin length corresponding to the perimeter of the solid.

2D periodic problem
As the third example, we treat a 2D micropump design problem. Figure 10 shows the layout and boundary conditions of the micropump design problem, in which the fluid motion is driven by the heated boundary condition on the lower left side. A similar design setting has already been investigated by Alexandersen et al. (2014) for a steady-state problem; on where, T l represents the temperature on the lower right side, and T h (> T l ) represents the amplitude of the temperature variation. Besides, t p denotes a time for the one cycle. In this study, T l and T h are set to 0 and 1, respectively. Note that the reference temperate difference is defined as ΔT = T h − T l = 1 for calculating Ra. The reference temperature T ref in (7) is given by T h + T l ∕2 = 1∕2 . In addition, the values of Prandtl number and Rayleigh number are set to Pr = 1.0 and Ra = 1.0 × 10 3 , respectively. Let us introduce the following objective functional for designing an unsteady micropump under the periodically steady-state interval I 1 = [t † , t 1 ] discussed in Sect. 3.3, as follows: where, and represents a weight parameter varying from 0 to 1. Note that (76, 77) has the negative sign to replace the maximization problem with the minimization one. Here, the first term on the right-hand side in (75) represents the time-averaged mass flow through the boundary Γ MF , and the second term represents the standard deviation of the mass flow. Therefore, minimizing (75) intends to enhance the mass flow as well as to reduce its fluctuation.
Since it is difficult to identify a periodic steady-state condition, we used a fixed t † = 1.0 × 10 6 determined by our preliminary investigation for this numerical example. As shown in Fig. 11, it can be confirmed that the mass flow j t achieves the periodic steady-state.
As a constraint, we introduce where the upper bound of the fluid volume rate is set to v f = 0.5 . Besides, the initial configuration is set to = 0.5. In contrast to the heat sink design problem discussed in Sects. 4.1 and 4.2, we do not use the filter and projection scheme to the miropump design problem, as there is no chance to occur checkerboard pattern in the obtained designs. On the other hand, to reduce the intermediate relative densities, the parameter q in (8) is initially set to 1.0 × 10 −6 and multiplied by 10 at every 100 optimization steps or until the convergence criterion in (55) is satisfied. This procedure is repeated until q = 1.0 × 10 −2 and (55) are met. Figure 12 shows the comparison result of the optimized configurations for three cases of the weight parameter, i.e., = 0.25, 0.50, 0.75 . Here, the time of the one cycle is set to t p = 5.0 × 10 4 . The influence of on the obtained configurations can be interpreted as follows. When the parameter is a relatively small value, J in (75) tends to emphasize the standard deviation of the mass flow across Γ MF . Therefore, more solid parts should be located around the heated boundary to reduce the fluctuation of the mass flow through Γ MF . On the other hand, when the parameter is a relatively large value, J in (75) tends to emphasize the amount of the mass flow across Γ MF . In this situation, the flow channel must be heated quickly to boost the upwind flow by natural convection to enhance the amount of the mass flow through Γ MF . Consequently, the optimized shape of the flow channel becomes symmetrical about the vertical midplane. Figure 13 represents the variation of the mass flow j t ( ) , in which = t − t † is defined, observed at the cross-section Γ MF in the optimized configuration for the one cycle t p . In addition, the dashed lines in Fig. 13 represent the time-average value of j t ( ) for the three cases and are given by As shown in Fig. 13, the oscillation j t ( ) and each average value j t increase as the weight parameter increases. Figures 14 and 15 show variations of the temperature distribution with the velocity vectors in the optimized configurations for two cases of weight parameter = 0.25 and 0.75, respectively. For the case of the lower weight = 0.25 , the flow channel is located to circumvent the heated boundary. This is because the main purpose of the case of = 0.25 is to reduce the variation of j t , which corresponds to avoiding the region of the large fluctuated T. Consequently, almost a uniform norm of the flow velocity can be observed throughout the channel, as shown in Fig. 14. On the other hand, for the case of the higher weight = 0.75 , the fluid passing through the channel is heated effectively, and a large amount of mass flow is transformed to the observation line Γ MF . It should be noted that the biases on the temperature are observed in Figs. 14 and 15, since the distribution of the   Table 4 shows the crosscheck for Fig. 12, and it can be confirmed that the optimized configurations for the periodic problem perform better than other for its particular values of . Figure 16 shows the comparison result of the optimized configurations for three cases of time of the one cycle, i.e., t p = 1.0 × 10 5 , 2.5 × 10 4 , 5.0 × 10 3 . Here, the weight parameter is set to = 0.5 . As shown in Fig. 16, the obtained configuration becomes symmetry about the vertical midplane as t p decreases, which is close to the configuration for the higher case of as seen in Fig. 12c. This result indicates that t p would affect the ratio of the standard deviation of the mass flow j t (t) to the time-average of j t (t) . In fact, the ratio of the second term to the first term in (75) for the three cases of t p , i.e., 1.0 × 10 5 , 2.5 × 10 4 and 5.0 × 10 3 are 0.238, 0.0867 and 0.0135, respectively. Such a response can be understood from the physical point of view as follows. As t p decreases, the boundary condition applied on the lower left side is close to the constant temperature condition, which corresponds to a steady-state problem. As a result, the standard deviation of the mass flow decreases. On the other hand, as t p increases, much time is given for developing convection, so the standard deviation of the mass flow increases. Table 5 shows the crosscheck for Fig. 16, and it can be confirmed that the optimized configurations for the periodic problem perform better than other for its particular values of t p .

3D transient problem
As the last example, we investigate the applicability of the proposed method for a 3D heat sink design problem concerning transient natural convention, while the 2D case has already been discussed in Sect. 4.2. In addition, we   demonstrate that our transient formulation is effective for rapid cooling concerning natural convection, as discussed in the previous work on forced convection (Zeng et al. 2020). The objective functional J and the related parameter settings are similar to those of the 2D problem. The volume constraint in (70) is imposed using v s = 0.05 , and the initial configuration is set to = 0 that corresponds to the fulfilled solid in D and is effective in quickly spreading the temperature distribution to the entire analysis domain. Herein, the heat flux applied to Γ SC is set to q 0 = 1.5 × 10 −2 . Figure 17 shows the problem layout with dimensions and boundary conditions for the analysis domain. Note that a similar problem has already been investigated by Alexandersen et al. (2016) for steady-state problems. The different setting of our study compared with that of the previous work is the design domain in contact with the bottom boundary. This is because we suppose that the 3D problem is a natural expansion from the 2D problem setting shown in Fig. 2. For effectively solving the 3D problem concerning unsteady natural convection, we use a supercomputer, SQUID (Supercomputer for Quest to Unsolved Interdisciplinary Datascience), owned by Osaka University. This study uses 128 CPU nodes with 64 cores per 1 node, and the LBM and the ALBM are implemented using the message passing interface (MPI) for parallel computing. For reference regarding the computational time in the following numerical examples, one optimization step takes about 35 s for t 1 = 5.0 × 10 4 case and 90 s for t 1 = 1.0 × 10 5 case. In this 3D problem, the maximum memory requirement for deriving the sensitivity is approximately 16TB, under about 1.1 × 10 6 lattice points. Due to the usable computational resources, the maximum number of optimization iterations is set to 500. In addition, the calculation domain is given as  Fig. 17, using the symmetrical boundary condition. Figure 18 shows the optimized configurations with the temperature distribution on a x-y plane for different Ra and t 1 settings. Besides, Fig. 19 shows the temperature distribution along with the streamlines developed around the optimized heat sink, and Fig. 20 is top views of the temperature distribution on each optimized configuration. From these results, it can be confirmed that the optimized shape changes significantly with increasing Ra and t 1 . In fact, when the observation time t 1 is large or the number Ra increases, the number of branches decreases, which can be found by comparing Figs. 20a and b, and the optimized shape forms higher vertical plates with a larger surface area standing along with the upwind flow. In other words, a leaf-like shape can be obtained as the optimized design under the larger Ra and t 1 settings. Figure 21 shows the relationship between t 1 and U cn determined by (72) for two Ra settings, including the optimized shapes at each t 1 setting. Similar to the 2D problem discussed in Sect. 4.2, the optimized designs depend on the magnitude of the heat convection effect. That is, the branched tree-like structure is obtained for the smaller settings of Ra and t 1 that correspond to a situation on the heat conduction dominant. On the other hand, for the larger settings of them, the leaf-like structure is obtained owing to the heat convection dominant.
Next, let us compare with the result of a steady-state case. Figure 22 shows the result of steady-state case for Ra = 1.0 × 10 5 setting. Herein, the common parameters with the transient case are identical, and the specific ones for the steady-state problem are similar to the 2D problem in Sect. 4.2. The volume of both the optimized designs is almost the same since the volume constraint is active. Figures 23 and 24 show the time-variation of the average temperature on the heated boundary Γ SC located at the centre bottom of the heat sink. Here, "Steady" represents the result for the optimized design obtained on the steady-state condition, and "Unsteady" is the transient problem under the same Ra settings with the steady-state case. Figure 23 and 24 show that the average temperature is converged around 50,000-75,000, namely, 50-75% of the maximum time steps in the transient problem. From these results, it can be found that the average temperature T for the transient case is always smaller than that for the steady-state case up to t = t 1 . On the other hand, these results are close to each other when t > t 1 , and the magnitude relationship is finally reversed. This result implies that if one wants to lessen the temperature on a heated   Average temperature at the heat source boundary of optimized configurations for Ra = 1.0 × 10 5 in the 3D steady-state and unsteady heat sinks. The ending time of the unsteady case is set to t 1 = 5.0 × 10 4 plate quickly, a heat sink composed of slender trees with branches is effective to enhance the heat transfer via conduction; on the other hand, if one wants to lessen the temperature gradually but sufficiently, a heat sink composed of wider leaves is effective to enhance the heat transfer by convection.

Conclusion
In this study, we proposed a topology optimization method using the LBM for unsteady natural convection problems. We introduced the ALBM, in which the discrete velocity Boltzmann equations are employed to derive the design sensitivities. In addition, we presented detailed numerical implementation strategies, including approximate procedures to solve efficiently. As a verification check for our proposed method, we solved a steady-state problem and confirmed characteristics of the result being as same as those in previous works using the FEM. To demonstrate the efficacy of the proposed method, we tackled several examples concerning unsteady natural convection and achieved the following: (1) In a 2D transient problem, we showed that the cumulative norm of flow velocity could be used to explain the dominant heat transfer mechanism on each optimized configuration. That is, the larger the Rayleigh number or the longer the observation time is considered, the more dominant mechanisms change from conduction to convection in optimized configurations. (2) We proposed an objective functional defined as a weighted sum of the time-averaged flow rate and the deviation to solve a 2D periodic problem. Under using the proposed objective functional, we revealed that the weight parameter affects the optimized configurations, whose shapes can be explained from the physical viewpoint.
(3) We demonstrated the applicability of the proposed method in a 3D transient problem, in which the analysis domain is discretized using about one million grid points with up to one hundred thousand time steps. Furthermore, we showed that the proposed method could lead to better results than the steady-state case in terms of the cooling speed.

Fig. 24
Average temperature at the heat source boundary of optimized configurations for Ra = 1.0 × 10 5 in the 3D steady-state and unsteady heat sinks. The ending time of the unsteady case is set to t 1 = 1.0 × 10 5