key: cord-0055696-238483s3
authors: Yoon, Gil Ho; So, Hongyun
title: Development of topological optimization schemes controlling the trajectories of multiple particles in fluid
date: 2021-01-27
journal: Struct Multidiscipl Optim
DOI: 10.1007/s00158-020-02817-8
sha: 463be9845c8d673c0ec23f38ebcda6cfc9202052
doc_id: 55696
cord_uid: 238483s3

This paper describes the development of a new topology optimization framework that controls, captures, isolates, switches, or separates particles depending on their material properties and initial locations. Controlling the trajectories of particles in laminar fluid has several potential applications. The fluid drag force, which depends on the fluid and particle velocities and the material properties of particles, acts on the surfaces of the particles, thereby affecting the trajectories of the particles whose deformability can be neglected. By changing the drag or inertia force, particles can be controlled and sorted depending on their properties and initial locations. In several engineering applications, the transient motion of particles can be controlled and optimized by changing the velocity of the fluid. This paper presents topology optimization schemes to determine optimal pseudo rigid domains in fluid to control the motion of particles depending on their properties, locations, and geometric constraints. The transient sensitivity analysis of the positions of particles can be derived with respect to the spatial distributed design variables in topology optimization. The current optimization formulations are evaluated for effectiveness based on different conditions. The experimental results indicate that the formulations can determine optimal fluid layouts to control the trajectories of multiple particles.

This paper describes the development of a new topology optimization scheme that manipulates and controls the trajectories of multiple particles in laminar flow. Controlling the motion of particles suspended in a carrier fluid has been an important subject in engineering. Several relevant studies can be found pertaining to microfluid, micromanufacturing, biology, and medicine (Qiu et al. 2015 (Qiu et al. , 2018 Wylie and Koch 2000; Mo and Sangani 1994; Lee and Balachandar 2010; Dinler and Okumus 2018; Wu and Yang 2019; Makhija et al. 2012) . Viscosity is widely employed to control the motion of particles. Depending on the size and Responsible Editor: Emilio Carlos Nelli Silva Gil Ho Yoon ghy@hanyang.ac.kr 1 School of Mechanical Engineering, Hanyang University, Seoul, South Korea mass of particles, the viscosity and inertial forces act differently. These differences can be employed to control the trajectories of particles in fluid. This type of technique has become a major standard to manipulate the motion of particles with additional forces (Yoon 2020) . To the best of our knowledge, there have been few studies on the topology optimization of structures or channels to manipulate the motion of particles (Yoon 2020; Andreasen 2020 ). Therefore, the aim of this study is to conduct transient sensitivity analysis for the location of particles and develop new topology optimization formulations to analyze the motions of particles in steady-state laminar flow (Fig. 1) . There have been a few studies regarding the topology optimization of the trajectories of particles in fluids (Yoon 2020; Andreasen 2020; Damiri and Bardaweel 2015; Prohm et al. 2013; Pagano et al. 2014) . After the introduction of structural topology, the aforementioned optimization scheme has been applied to several engineering applications (Bendsøe and Kikuchi 1988; Wang et al. 2003; Zhang et al. 2018; Maute and Reich 2006) . Without initial optimal topologies, the topology optimization scheme can provide the optimal layouts of various engineering structures (Zhang et al. 2018 ; Evgrafov et al. 2008; Dede et al. 2014; Yaji et al. 2016; Deng et al. 2011) . The application of this optimization to fluid problems has been extensively studied. Moreover, the topology optimization of turbulent fluids was studied previously (Papoutsis-Kiachagias and Giannakoglou 2016; Yoon 2016) . The sensitivity analysis and topology optimization for unsteady flow problems were formulated and presented (Deng et al. 2011 ). The multiphysics system of the conjugate heat transfer problem was considered (Yoon 2010; Dede et al. 2014; Yaji et al. 2016) . A review of topology optimization for the conjugate heat transfer problem can be found (Dbouk 2017) . The present study extends the work presented to consider the trajectories of multiple particles (Yoon 2020) . In this study, we focused on specific issues that consider multiple particles.

Subjects related to controlling particles are considered to be important and complex in several science and engineering fields, e.g., microfluidics systems, micromanufacturing, biology, lab-on-a-chip, and medicine. The motions of particles are mainly influenced by the drag force in laminar flow or external forces such as electrostatic, electromagnetic forces, and micro-scale hydrodynamic effects (Bockelmann et al. 2012; Hu et al. 1992; Wu and Yang 2019; Yoon and Park 2010) . Several numerical approaches have been developed to analyze the motions of particles (Wylie and Koch 2000; Mo and Sangani 1994; Lee and Balachandar 2010; Wu and Yang 2019) . The analysis of this multiphysics system has challenges; e.g., the positions of multiple particles and the effect of the particles on fluid should be mutually coupled. Owing to difficulties in mutual coupling, assuming small particles, it is possible to adopt a one-way coupling approach, i.e., fluid to particle only. Not fully accounting for the influence of particles on the surrounding fluid, Eulerian-Lagrangian approaches for fluid and particle, respectively, have been employed for the simulation purposes. This one-way coupling approach is beneficial for the topology optimization for forward and sensitivity analyses (see Andreasen 2020 for the rigorous sensitivity derivation). This one-way coupled multiphysics equation (Eulerian-Lagrangian approach) can provide details of the flow and enable interactions with small particles.

Based on the one-way coupling multiphysics analysis, this study develops new topology optimization schemes to split, capture, and control the trajectories of multiple particles. The topology optimization of particles considering the fluid velocity as the objective function is reported (Yoon 2020) . A similar topology optimization method for the design synthesis of microfluidic particle manipulators was independently presented with rigorous sensitivity derivation (Andreasen 2020) . The controlling, separating, and gathering of multiple particles via the design of topologically optimized pseudo rigid bodies are considered in the current study. In the Navier-Stokes equation, Darcy's force is applied to model fluid in porous media and is parameterized with respect to spatial design variables using the Solid Isotropic Material with Penalization (SIMP) (Yoon 2020) . By modifying the multiphsyics equations and the transient sensitivity analysis for the trajectories of multiple particles, it is possible to determine the optimal distributions of porous media.

The paper is organized as follows: Section 2 provides a background pertaining to the coupled analysis of particles and steady-state fluid motions and the development of the sensitivity analysis of the objective function with particle motion. Section 3 describes several topology optimization formulations and studies. Section 4 presents the conclusions of the research and provides suggestions for future studies.

This section describes the development of the topology optimization scheme and the multiphysics analysis coupling between the Navier-Stokes equation and Newton's second equation of particles. The finite element procedure is implemented to model the steady-state laminar flow, and the motions of particles are solved using Newton's second equation. The objective of the study is set as a function of the final positions of the particles. Compared with the relevant studies that consider the velocity of a particle, this study considers the positions of multiple particles focusing on the control of the particle locations (Yoon 2020 ). Newmark's scheme is employed as a numerical solver for the transient analysis of the particles. To model the rigid domain in fluid, the pseudo force for the porous material model is applied to the Navier-Stokes equation. Transient sensitivity analysis of the positions of multiple particles is performed.

Several innovative numerical theories have been developed to track the motions of particles in fluid. This study assumes that the particle size is small and does not affect the viscosity and the motion of the fluid. With larger particles or external forces, these effects cannot be negligible. The present study advances the effort made to consider and track the positions of the particles.

The steady-state laminar flow can be computed by solving the Navier-Stokes equation with Darcy's force, αu to set the velocities of pseudo-rigid domain close to zero.

with the following boundary conditions:

The fluid domain is denoted by Ω. The velocity, pressure, mass density, and dynamic viscosity are denoted by u, p, ρ, and μ, respectively. The boundary conditions are described either by the fluid velocity condition (no-slip boundary condition along Γ u 0 , the inflow/outflow boundary condition along u * ) or the pressure boundary condition along Γ p * .

To model the pseudo-rigid domain, it is possible to include Darcy's force with the SIMP method.

In the Navier-Stokes equation (1), the Darcy force αu is included with the maximum Darcy force α max . For topology optimization, this fictitious force is interpolated with respect to the spatial design variables γ with the SIMP n. Domains discretized with finite elements with small values for γ correspond to fluid domains where domains with ones for γ correspond to pseudo-rigid domains (Yoon 2020; Andreasen 2020; Dinler and Okumus 2018; Qiu et al. 2015; Damiri and Bardaweel 2015; Pagano et al. 2014; Deng et al. 2011 ).

The domain with ones for the design variables is defined as the pseudo rigid domain and the deformation of the solid domain is not considered. In the related applications, i.e., microfluid, the streamline of fluid mainly moves particles and the deformation of structure is neglected. Therefore, this research also ignores the structural deformation. Newton's second law can be analyzed with respect to the fluid force by neglecting the rotational motions of particles, as given below:

where the drag force coefficient is denoted by F D . The mass, diameter, density, and velocity of a particle are denoted by m p , d p , ρ p , and v, respectively. The absolute coordinate of the particle is denoted by X. The drag force inside the laminar fluid is formulated as the function of the difference between the velocities of particle u and fluid v. The proportional constant of the force in (4) is m p F D .

Based on the product of the difference and the constant, the direction of the drag force is opposite the direction of the relative motion of the particle to fluid. Several experimental and theoretical models have been developed (Hu et al. 1992; Kulkarni and Morris 2008; Walsh 1976; Bagheri and Bonadonna 2016) . For the stability of the analysis of Newton's equation, a smaller drag force is used. The velocity vector is numerically integrated in order to evaluate the loci of the particles.

The objective of the optimization problem is to manipulate the loci of particles depending on the size and mass. Therefore, it is intended to maximize the control capability of a device by structural optimization. To achieve this objective, several optimization formulations are considered. To achieve the objective, several simulation conditions are assumed. The first assumed condition is the steadystate laminar flow. Several particles of different sizes and masses are placed in one location or different locations, and they follow a streamline of laminar flow. Depending on the properties of particles, they show different loci that result from the control of the particles. To determine the locations of the particles, the transient forward analysis and the transient sensitivity analysis are performed and the contact conditions among particles are ignored. In the transient simulation of the particles, the consideration of the contact increases the complexity in the forward and the sensitivity simulations. Thus, the contacts among particles toward the pseudo-rigid body and the boundary condition are not considered.

For the topology optimization considering the motions of particles, the objective function should be formulated with the position vector, which is the integrated velocity. The velocity of the particle is considered, and to consider the positions of multiple particles, a new optimization formulation and its sensitivity analysis should be derived (Yoon 2020) . The positions of the particles can be defined as follows:

The velocity vectors and the initial locations of the particles are denoted by v and X init , respectively. The simulation time is denoted by t f . The Lagrangian L over a time period is defined by the Lagrange multipliers. The initial position of the particle is denoted by X init .

where the Lagrange multipliers are denoted by λ for Newton's equation and ψ for Navier-Stokes equation. The objective function is a vector quantity rather than a scalar quantity. For the sensitivity, the following differentiation of the Lagrangian is obtained.

The differentiation is summarized as follows:

The first term of the second equation above is integrated as follows:

After separately gathering the terms involved for ∂v ∂γ e and ∂u ∂γ e and letting the involved factors be zeros, the following transient sensitivity analysis with adjoint sensitivity equations in (13) and (14) can be derived.

Adjoint system 1:

withλ = 0 and λ = 0 at t = t f based on the reversal time integration scheme.

Adjoint system 2:

The first adjoint variable λ can be solved for the adjoint sensitivity equation in (13), and it is used to solve the second adjoint system equation in (14) for ψ. The first Lagrange multiplier λ in (13) is the solution of the transient system with the time reversal. The Newmark integration scheme has been used for the solutions (see Yoon 2020 for more details). For large particles, the contact conditions between the particles become important. To consider the contact conditions among particles, the governing equatinos should be modified. For an example, Newton's second law should be modified by adding the contact force as follows:

where the contact force is denoted by F contact whose direction and magnitude depend on the displacement, velocity, design parameters, and parameters of the numerical integration scheme. Therefore, to consider the contact condition, the new design sensitivity analysis should be derived considering these parameters and conditions. To simplify the derivations, this research does not consider the contact condition. The finite element procedure with 9-node element is implemented for the fluid motion with the Newton-Raphson solver. The particle motion is solved by the Newmaker scheme (beta = 1/6, gamma = 1/3) considering the current particle position and its associated fluid velocity. The incremental time in the Newmark scheme is an important. A small value for this incremental time is selected to maintain the stability in the numerical simulation of particles. The time reversal scheme is used for the first adjoint system.

The inversion of the tangent stiffness of the second adjoint system is carried out one time owing to the steady-state fluid motion. One crucial computational procedure is the step to find out the particle velocity of particles. In the current implementation, after finding the elements containing particles, the corresponding natural coordinates at the current positions of the particles are calculated. These natural coordinates are used to calculate the derivative of the force with respect to the velocity. As the particles do not interact with each other, Matlab's parallel approach can speed up the procedure. In our simulation, a HPC cluster with 100 GB memory and 72 cores is used with the parallel toolbox in Matlab.

When investigating the dynamics of multiple particles or polyatomic molecules, geometric constraints must be imposed among particles in the simulation. For example, bond-length and bond-angle constraints among particles can be imposed for the particle dynamic simulation. Several approaches have been developed, and a method should be chosen considering the accuracy of the simulation and the computation time. In this study, the fast SHAKE algorithm is implemented to impose constraints among particles (Kräutler et al. 2001) . The purpose of the SHAKE algorithm is to add the constraint force, which imposes geometrical constraints in Newton's equation. The N c distance constraints can be defined as follows:

where the distances between the particle k 1 and k 2 involved in the k-th constraint and the corresponding constraint distance are denoted by r k 1 k 2 and d k 1 k 2 , respectively. Based on the constraint condition, the Newton equation is modified considering the constraint force:

where the Lagrange multipliers for the constraint conditions are denoted by l k . The second term on the right-hand side of (17) represents the constraint force, f c i . To determine the constraint force and the Lagrange multiplier, the SHAKE algorithm is employed (Kräutler et al. 2001 ). In the SHAKE algorithm, the updated displacements without the constraint forces are denoted by X uc i (t + Δt). Once the constraint forces are added, the displacements are re-updated to

The condition that these new displacements should satisfy the constraint equation is used to determine Lagrange multipliers, l k (t):

where X k 1 k 2 is defined as X k 1 − X k 2 Then, the displacement can be updated considering the Lagrange multipliers. Owing to the aforementioned approximation in the determination of X i (t + Δt), the above procedure should be iterated , fluid: density = 1000 kg/m 3 , dynamic viscosity = 352 × 10 −6 Pa×s, particle: mass 1 = 1.0865 × 10 −13 kg (radius = 2.39 μm, density = 1900 kg/m 3 ), mass 2 = 1.0865 × 10 −11 kg (100 times larger than that of the first particle), F D1 = 3.1715 ×10 −12 N·s kg·m , F D2 = 1.5857 × 10 −12 N·s kg·m , mass 0 : 30% of the design domain discretized by 120 by 180) and b the trajectories of the two particles in no-structure An optimization result when separating the particles. a An optimized layout with 30% mass t f = 150 s with 20,000 time steps (n = 4, α max = 10 4 , objective function: −1.58019 ×10 −2 m), and b fluid velocity distribution and the particle trajectories (green: particle 1 (mass 1 = 1.0865 × 10 −13 kg, F D1 = 3.1715 × 10 −12 N·s kg·m ) and red: particle 2 (mass 2 = 1.0865 × 10 −11 kg, F D2 = 1.5857 × 10 −12 N·s kg·m )), b the fluid velocity, and c the optimization history until the convergence of the corrected displacement. In the adjoint equation, the Lagrange multipliers are considered as the stiffness term. Although the Lagrange multipliers are a function of the displacement, this study assumes that they are constants in the adjoint equation as the Lagrange multipliers that converged in the iterative procedure of the SHAKE algorithm are insensitive to the displacements. This assumption simplifies the implementation.

To show the application of the present topology optimization and explore its limitation in designing optimal layouts, which control the motions of particles, this section presents optimization formulations whose objective function and constraints and their derivatives can be obtained from the formulations given in the previous section. To observe the effect of different material properties and maintain numerical stability, the values of mass and drag force are arbitrarily modified. The simulation time and time steps are selected based on the stability of the motion of particle. For the topology optimization, several optimization formulations considering the locations of several particles are presented. Therefore, using the derived sensitivity analysis using the Lagrange multipliers, it is possible to calculate the sensitivity of the objective function and constraint. The present study does not apply sensitivity filter or density filter in order to focus the optimization problem itself. It is possible to apply some filters to control feature of optimum design and to consider the manufacturability. To solve the optimization problem, the method of moving asymptotes (MMA) algorithm is employed (Svanberg 1987 ).

For the first optimization example, the topological optimization problem separating two heterogeneous particles with different masses and drag forces but with identical release positions is considered (Fig. 2) . For the fluid boundary condition, the parabolic distributed fluid (the inlet boundary condition) is defined at the bottom left boundary of the design domain and the fluid outlet boundary condition (the pressure boundary condition) is defined at the upper left boundary. In order to show the effect of the inertia force and the drag force on the optimum design, some artificial material properties are assumed here. It is assumed that the first particle is set lighter than the second particle, which significantly changes the inertia force but the drag coefficient of the first particle is 2 times larger that of the second particle. Although they are released in the same location (0.54 cm, 0.5 cm) just after the fluid inlet, their trajectory responses differ from one another in the fluid. This may be regarded as a limited application of the inertia micro fluid device. The almost identical initial trajectoreis of the particles are shown in Fig. 2b . To efficiently separate them, the optimization formulation in (20) is formulated and solved. The objective function in (20) is the difference in the x positions of the two particles. The mass constraint and the difference in the y positions are imposed as the constraints. The second constraint in (20) is imposed for the efficient movement of the particles. Without this second constraint, the motion of one particle may be frozen to maximize the objective function, which we do not intend to. By distributing spatial materials at the fluid design domain and increasing the objective function, the particles can be separated effectively. (20) where the x positions of the particles are denoted by p 1

x and p 2

x and the y positions are denoted by p 1 y and p 2 y . The N e design variables are γ and the allowable mass is denoted by mass 0 . Note that the mass constraint is added to obtain a local optimum.

To maximize the objective function satisfying the constraints, the optimal distributions of large Darcy's force or the pseudo solid domain should be designed through the topology optimization process. By solving (20), a complex channel utilizing the difference in the inertia forces of the particles can be obtained. Figure 3a shows an optimized layout obtained via the present optimization formulation with a simulation time of 150 s. The loci of the particles in this system smoothly incline toward the outlet. The complex curved channel shape above the initial location can be obtained. Just before the initial release, the expansion chamber shape maximizes the fluid velocity (see Fig. 3b ). Subsequently, the curved expansion chamber makes a streamline deflect toward the upper domain to detour the particles with the different inertia forces. Then, the upper structures maximize the difference in the x positions with similar y positions imposed by the second constraint. Some porous domains with various permeability values appear. Figure 3c shows the optimization history.

To further illustrate the feature of the present formulation, the different initial locations are considered in Fig. 4a   Fig. 4 Optimization results when separating the particles. a An optimized layout with the different locations (initial locations: (1 cm, 0.5 cm)) and b an optimized layout with different locations (initial locations: (2 cm, 0.5 cm)) and b. The spatially optimized fluid motions optimized for each condition can be obtained by solving the optimization problem using the present topology optimization. The optimized layouts use the inertia forces by developing the curved layouts. Some isolated particles appear on the right side of Fig. 4b . This may be because the present formulation optimizes the final positions of the particles without considering the fluid motion. After several iterations with intermediate design variables, materials are distributed in that region. As the optimization process continues, some materials gather, and the optimization algorithm does not remove them as that region is not sensitive to the objective function or the constraints. In addition, the fixed mass constraint is added to penalize the intermediate design variables. By maximizing the velocities or minimizing the head loss, small particles may disappear, but the current study focuses on the motions of the particles. These problems show that the present topology optimization framework can be an effective tool to design a channel that controls or sorts particles.

For the second example, structural topology optimization to gather multiple particles is considered, as shown in Fig. 5 . This optimization intends to find an optimal layout that gathers identical or heterogeneous spatially distributed particles. It is intended to modify the flow to make the particles gather at an area whose location should be determined by the optimization. The design domain and the boundary conditions are shown in Fig. 5 .

For this purpose, the optimization problem in (21) is proposed herein. To control the loci of the four particles, the variations in the coordinates, i.e., σ x for the variation in the x coordinates and σ y for the variation in the y coordinates, are considered in the optimization formulation. The optimization formulation "OP 2X" sets the variations in the x coordinates of the particles as the objective function with the mass constraint and the variations in the y coordinates of the particles as the constraint. The allowable limit of the variation is set to Δ. In addition, it is possible to recombine the above objective function and the constraints for a specific purpose. Using the derived adjoint equations in the previous section, it is possible to obtain the sensitivity values of the above objective and the constraint with respect to the design variables. To demonstrate the above optimization formulation, the four particles of equal mass (mass = 1.0865 ×10 −13 kg and F D = 3.1715 × 10 −12 N·s kg·m ) are initially located at (0.0054 m, 0.006 m), (0.01 m, 0.006 m), (0.02 m, 0.006 m), and (0.03 m, 0.006 m), respectively, in the domain. The simulation time is set to 300 seconds with 40,000 time steps. An optimal layout to gather the particles using the optimization formulation in (21) is shown in Fig. 6a . The trajectories are plotted with different colors. The loci of the particles in this system tend to incline toward the triangular structure. As the particles are gathered after 300 s, an optimal layout with longer trajectories inversely prorated according to the initial position of x is obtained. The left larger triangular structure and the bridgelike structure beneath the triangular structure create the flow channel to optimize the trajectories of the particles. The small arc structure appears along the end of the trajectory of the first particle to minimize the fluid leakage to the outlet. This small arc structure helps the leftmost particle to travel along with the other particles. The bottom bridgelike structure should be post-processed as the physical interpretation of the inner holes is meaningless. In other words, the outer rim of the bridge structure dominantly changes the direction of the flow to optimize the fluid motion for the given optimization formulation. It appears owing to the modeling method of the pseudo rigid body, which includes Darcy's force. This aspect also can be addressed by considering the manufacturing constraint. The trajectories of the second and the third particles are almost the same. Figure 6c shows the optimization history.

Separated domains can appear in the current optimization framework with the modeling of three-dimensional fluid channel with the two-dimensional simulation. Isolated domains can be fabricated using some fabrication techniques such as etching, milling, or CNC machining. To study and explore the effects of the isolated structure on the triangular structure, the post-processing of the structure is performed, as shown in Fig. 7 . After removing the structure, the fluid motion is modified such that the particles are spread, as presented in Fig. 7 . As displayed in Fig. 7c , an investigation of the trajectories reveals that the first particle moves further toward the left side after the structure is removed. This implies that the present optimization scheme finds out the auxiliary structure to gather the particles. In addition, the effect of the mesh refinement is tested with the different mesh refinements in Fig. 8a and b. The similar layouts are obtained with a little subtlety between them. In Fig. 8c , the post-processed designs especially the bridge-like structure and the trajectories of the particles are tested. The similar trajectories can be obtained. Fig. 7 Post-processing of the design in Fig. 6 . a The post-processed design and the trajectories, b a comparison of the trajectories, and c the trajectory of each particle Fig. 8 Effect of fluid mesh refinement and the post-processing. a An optimal design with a 120 × 90 mesh, b an optimal design with a 160 × 120 mesh, and c the post-processing result Note that the optimal layouts with the isolated structure can be manufactured.

The optimization problem with the equivalent formulation and boundary condition without considering the initial locations and mass of the particles is tested, as shown in Fig. 9 . The initial locations of the particles are set to (0.0054 m, 0.002 m), (0.01 m, 0.002m), (0.02 m, 0.0005 m), and (0.02 m, 0.0003 m), respectively. In addition, the first two particles are 10 times heavier than the other two, and this naturally induces different initial forces. With this configuration, an optimized layout in Fig. 9a can be obtained. The feather-like structures appear at the end of the triangular structure. The bill-shaped void part appears at the end of the triangular structure, which appears to control the fluid motion for smooth transition of the third particle.

Moreover, an optimization with 10 particles is tested using the proposed formulation (Fig. 10) . The conditions, except the number of the masses, are the same as those of the first test of the second example. By increasing the number of particles, it appears that the fluid channel between the triangular and bottom structures becomes more elaborate and smoothed. The upper arc structure is attached to the main triangular structure to control fluid motion, thereby preventing the spread of the particles. Figure 11 shows an optimal layout with 100 identical masses to demonstrate the potential, capability, and limitation of the present algorithm and optimization framework. As our numerical code may not be optimized with limited computing resources, a shortage of memory occurs with a long simulation time by increasing the number of involved masses. In particular, the first Lagrange multiplier requires a large amount of physical memory. To address this, the simulation time is reduced to 50 s with 10,000 time steps for the optimized layout in Fig. 11 . The optimized channel is similar to the one in the previous example but with smoother boundaries. As information about the final positions of the particles is only considered in the optimization framework, the greater the number particles and the involved information, the smoother and clearer are the optimized channels. The upper arc structure in Fig. 10 does not appear as the particles do not travel to the upper part. In a cluster computer with sufficient memory and cores, the simulation time is increased again for the optimal layout in Fig. 12 . The optimized shape is similar to that of Fig. 10 , with an upper arc structure. The channel shape becomes smoother in this design. This illustrates that the richness of the trajectories and the dataset allows the optimization algorithm to refine the design where fluid velocities control the particle velocity and the final locations of the particles.

Additionally, the geometric conditions imposed by the SHAKE algorithm are considered in the optimization problem of the second example in Fig. 5 . With the same optimization formulation in (21), the three pairs of two particles are linked with the geometric conditions in that the distances of each pair should be 0.002 m in Fig. 13a . This implies that optimization should be conducted while maintaining these distances. The optimized design of this optimization problem is shown in Fig. 13b . The trajectories of the pairs are indicated in different colors. This example demonstrates the optimized layout in a manner similar to previous examples with sharp edges. Figure 13c shows the optimization history. It demonstrates that the SHAKE algorithm works well with the present optimization algorithm. 

To demonstrate the geometric constraint in the present optimization formulation, a triangular constraint is imposed for three particles in Fig. 14. The addition of the geometric constraints would make it possible to impose the triangular constraint during analysis and design. The dimensions of the design domain, including the left bottom fluid inlet and the right upper outlet, are set as 4 cm × 3 cm in where x and y locations of the i-th particle are denoted by p i x and p i y , respectively. The objective of the above optimization problem is to maximize the location of the x coordinates of NM particles. Figure 14b shows an optimized layout with the optimization formulation. The channel shape appears in front of the pairs, and the expansion chamber is designed based on the present optimization formulation to guide the fluid motion toward the outlet. During the optimization process, the geometric constraints with constant distances are successfully imposed and maintained. In addition, a saw-like structure is formed at the outlet, and the outlet is partially covered to maximize the objective function. Figure 14c shows the optimization history. This demonstrates that the present optimization formulation works well with the SHAKE algorithm.

To further demonstrate the SHAKE algorithm in the present optimization formulation, 10 pairs with 20 particles are considered, as shown in Fig. 15 . As the initial locations of the particles are randomly chosen, the distances of the pairs are different. The optimal layout shown in Fig. 15 is similar to that with the triangular particles. To maximize the initial velocity, the channel collecting the inlet fluid is designed based on the present optimization formulation.

In the last example, the optimization problem related to switching or sorting two particles using different material properties and inertial forces is considered (Fig. 16) . This optimization problem focuses on utilizing the difference in the inertial forces exerted on the two particles for the sake of maximizing the distance in the y direction Fig. 16 Example 4. A channel design to switch particles (Reynolds number: 85.2273 (1000 × 0.03 ×1 × 10 −3 / 352 × 10 −6 ), n = 3, α max = 5 × 10 4 , fluid: density = 1000 kg/m 3 , dynamic viscosity = 352 × 10 −6 Pa×s, mass 0 : 20% of the design domain, the analysis domain discretized by 400 × 120) between the particles. A local optimal design restraining the motion of one particle can be obtained by considering only the distance. Thus, the two constraints with the x displacements of the two particles are included in the optimization formulation (24). The inlet boundary condition is defined at the center of the left side, and the right side includes two side outlets. The design is limited to the middle domain in Fig. 16 . One of the challenges in this optimization problem is to switch the particles first and then maximize the distance. Based on several numerical exercises, it has been found that the initial distance between the particles may not be large in an optimized layout that enables switching the particles with the given material properties. This is regarded as one of the main physical difficulties in this optimization. The initial distance is set to 0.5 cm in the y direction on the left side of the design domain. To obtain the differences among the inertial forces, the mass of the second particle is set to be 100 times greater than that of the first particle, mass 2 = 100 × mass 1 , and the drag force is set as 0.5 times greater than that of the first particle, F D1 = 2×F D2 . The optimization algorithm should challenge the differences of the material properties to find out an optimal layout for the optimization formulation in (24). The locations of the particles in the y direction are denoted by p 1 y and p 2 y , and those in the x direction are denoted by p 1

x and p 2 x . The lower bounds of the x positions of the particles are set to 3 cm. Without any constraints on the x displacements, few local optimal layouts trapping the motion of the heavier particle can be obtained. That is regarded as a solution for the optimization but not practical. In other words, the second and third constraints help the optimization algorithm determine a channel structure that smoothly moves and switches the particles. Figure 17a shows an optimal layout that separates the particles. It is interesting that a small triangular structure is obviously formed between the initial locations of the particles. This triangular structure separates the input fluid. The fluid above the triangular structure moves downward to guide the first particle whose trajectory is rendered in red color. The fluid underneath the triangular structure moves upward to guide the second particle whose trajectory is rendered in green color. Then, the upper and lower structures blocking the flow appear. As per our analysis, the upper and bottom structures near the triangular structure play an important role in controlling fluid motion and drag force. Figure 17b and c show the fluid velocity and the optimization history. The present structure is a local optimum. The number of particles, upper bound of the pseudo force, penalization in topology optimization, and lower bounds of the constraints mainly affect the optimized layout.

Subject to Mass ≤ mass 0 δ ≤ p 1 x (δ = 3 cm) δ ≤ p 2 x (δ = 3 cm) γ = [γ 1 , γ 2 , ..., γ N e ], γ min γ 1, γ min = 0.001 (24) As a next test, by assigning the same mass to the two particles initially separated 0.5 cm, an optimal layout switching and separating the particles is pursued in the example. A symmetric design shown in Fig. 18a can be obtained using a symmetric condition. Above and beneath the triangular structure, the funnel-shaped structure collecting the input fluid before the triangular structure and pouring the fluid after the funneled shape structure appears. As the fluid motion is guided by the optimal structure, the switching of the particles is successfully Fig. 17 An optimization result with two particles switching their trajectories (initial locations: (1 cm, 1 cm), (1 cm, 1.5 cm), particle: mass 1 = 1.0865 × 10 −13 kg, mass 2 = 1.0865 × 10 −11 kg, F D1 = 3.1715×10 −12 N·s kg·m , F D2 = 1.5857 × 10 −12 N·s kg·m ). a An optimized layout with 20% mass t f = 225 s with 30,000 time steps (objective function: −7.5149 ×10 −4 m), b fluid velocity distribution, and c the optimization history achieved. The aforementioned example shows that the Fig. 18 An optimization result with the identical two particles switching their trajectories (initial locations: (1 cm, 1.25 cm), (1 cm, 1.75 cm), particle: mass 1 = mass 2 = 1.0865 × 10 −12 kg, F D1 = F D2 = 1.5857 × 10 −12 N·s kg·m ). a An optimized layout with 20% mass t f = 120 s with 16,000 time steps (objective function: −1.2080686 ×10 −3 m), b fluid velocity distribution, and c the optimization history present optimization framework can effectively control the trajectories of several particles.

This study describes the development of topology optimization formulations for optimal layouts that gather, control, or separate multiple particles by manipulating the fluid drag force. Recently, studies on the motions of particles in fluids have become increasingly important. For example, SARS-CoV-2 can spread via coughing in a corridor, and it is important to control the aerosol cloud traveling in the corridor to prevent infection. To contribute to this subject from a structural optimization viewpoint, this study derives a transient sensitivity analysis for the trajectories of multiple particles and develops several topology optimization formulations considering the trajectories of particles in fluids. Steady-state laminar flow is assumed, and the effect of the particle on the surrounding fluid is neglected. Two Lagrange multipliers are required for the Navier-Stokes equation and Newton equation. As the force terms of the adjoint equation originated from Newton's second law (of particles) are time-varying, transient sensitivity analysis with the time reversal scheme is essential. The final locations of particles are mathematically considered in the objective function and constraints to control the trajectories. To show the validity of the present formulations, several demonstration problems are solved, and the optimal layouts to control the motions of the particles and fluid can be obtained. It has been shown that the richness of the trajectory and the sensitivity analysis allows the optimization algorithm to refine optimal layouts, whose fluid velocities can control the particle velocity and the final locations of the particles. The present topology optimization scheme can help control the trajectories of particles suspended in fluids. As part of future research, the motions of several particles must be optimized, and the contact between the particles and the wall needs to be considered to examine local particle concentration, migration, and segregation. Furthermore, the rudimentary theory of the multiphysics system can be improved. In addition, the mixing of particles in a 3D turbulent flow can be researched. 

Conflict of interest The authors declare that they have no conflict of interest.

The replications of the presented results require the FE-based Navier-Stokes equation and the Newmark scheme to solve the particle trajectory. The Newmark scheme can be found in the standard FE book and the Navier-Stokes equation with rectangular elements should be used to efficiently calculate the fluid velocities of fluid. The sensitivity analysis requires the time reversal scheme to compute the Lagrange multipliers.

A framework for topology optimization of inertial microfluidic particle manipulators

On the drag of freely falling nonspherical particles

Generating optimal topologies in structural design using a homogenization method

Optimization of an electrokinetic mixer for microfluidic applications

Numerical design and optimization of hydraulic resistance and wall shear stress inside pressure-driven microfluidic networks

A review about the engineering design of optimal heat transfer systems using topology optimization

Multiphysics simulation: electromechanical system applications and optimization

Topology optimization of unsteady incompressible navier-stokes flows

Inertial particle separation in curved networks: a numerical study

Topology optimization of fluid domains: kinetic theory approach

Direct simulation of fluid particle motions

A fast shake algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations

Pair-sphere trajectories in finitereynolds-number shear flow

Drag and lift forces on a spherical particle moving on a wall in a shear flow at finite re

Topology optimization of multi-component flows using a multi-relaxation time lattice boltzmann method

Integrated multidisciplinary topology optimization approach to adaptive wing design

A method for computing stokesflow interactions among spherical objects and its application to suspensions of drops and porous particles

Optimizing design and fabrication of microfluidic devices for cell cultures: an effective approach to control cell microenvironment in three dimensions

Development of topological optimization schemes

Continuous adjoint methods for turbulent flows, applied to shape and topology optimization: industrial applications

Optimal control of particle separation in inertial microfluidics

Microfluidic device for mechanical dissociation of cancer cell aggregates into single cells

Microfluidic channel optimization to improve hydrodynamic dissociation of cell aggregates and tissue

The method of moving asymptotes-a new method for structural optimization

Influence of particle drag coefficient on particle motion in high-speed flow with typical laser velocimeter applications

A level set method for structural topology optimization

An overview of numerical methods for incompressible viscous flow with moving particles

Particle clustering due to hydrodynamic interactions

Topology optimization in thermal-fluid flow using the lattice boltzmann method

Topological design of heat dissipating structure with forced convective heat transfer

Topology optimization for turbulent flow with spalart-allmaras model

Transient sensitivity analysis and topology optimization for particle motion in steady state laminar fluid

Topological design of electrode shapes for dielectrophoresis based devices

A moving morphable void (mmv)-based explicit approach for topology optimization considering stress constraints