Multiphysics and Multiscale Simulation

1.18 Multiphysics and Multiscale Simulation J. G. Korvink, A. Greiner, and Z. Liu, Department of Microsystems Engineering-IMTEK, University of Freiburg, Freiburg, Germany ª 2008 Elsevier B.V. All rights reserved.

1.18.1 1.18.1.1 1.18.1.1.1 1.18.1.1.2 1.18.1.1.3 1.18.1.2 1.18.1.2.1 1.18.1.2.2 1.18.1.2.3 1.18.1.2.4 1.18.1.2.5 1.18.1.2.6 1.18.1.2.7 1.18.1.3 1.18.1.4 1.18.1.4.1 1.18.1.4.2 1.18.1.4.3 1.18.2 1.18.2.1 1.18.2.1.1 1.18.2.1.2 1.18.2.1.3 1.18.2.2 1.18.2.2.1 1.18.2.2.2 1.18.3 References

Multiphysics Simulation and Design Optimization Multiphysical Essence of MEMS Simulation Classification of coupled physical fields in MEMS Domain-type multiphysical coupling Boundary-type multiphysical coupling Coupled FEM Simulation PDE expression of physical field Numerical simulation of single PDE using FEM Numerical simulation of coupled PDEs using FEM Sequentially coupled method Fully coupled method Linearization of nonlinear PDEs Time-dependent solvers Mesh Adaptation Multiphysics Design Optimization Structural layout optimization Simultaneous simulation and optimization method Derivation of the optimization equation Multiscale Simulation Particle Dynamics Methods Coarse-graining microscopic methods vs refinement of continuum models Dissipative particle dynamics Applications of DPD Kinetic Models Approximations of the BTE Applications of CDLBA Conclusion

1.18.1 Multiphysics Simulation and Design Optimization The microelectromechanical systems (MEMS) design process is highly iterative in character. This is because design constraints tend to become more precise as each design evolves, and because complexity issues emerge only as a design takes shape. Computational tools aim at alleviating this aspect by providing two key capabilities: A high level of predictability of the physical processes in a given design realization  A method to modify a design in the direction of optimality 

539 540 540 540 540 541 541 541 543 543 543 543 544 544 545 545 546 547 550 550 551 551 553 553 554 554 555 555

For the second capability, the first is of course a prerequisite. In this chapter we describe some recent techniques that contribute to more powerful tools both for multiscale multiphysics simulations and toward the ultimate goal of automatic design optimization. We have grouped the continuum approaches together and described them first, and have followed this up with the particle simulation approaches for soft matter simulations. As dimensions reduce and soft matter becomes more relevant to MEMS, we see a gradual increase in the use of particle-based methods, and their emergence in commercial tools. 539

540 Multiphysics and Multiscale Simulation

1.18.1.1 Multiphysical Essence of MEMS Simulation Modeling denotes the process of developing an appropriate mathematical description of a particular system. Simulation is the process of using this description to actually compute the behavior of the system analytically or numerically. Originally, silicon-based sensors and actuators constitute the basic building blocks of MEMS, and integrated MEMS employ materials and processing techniques originally designated for integrated circuits. Despite the inherent reference to electrical and mechanical phenomena, the field of MEMS has grown to encompass a much broader application range, including microsystems exploiting electrical, mechanical, thermal, magnetic, fluidic, selected chemical, and biological effects and possible coupling among these fields (Romanowicz 1998, Senturia 2000). In addition, a much wider range of material technologies has become established, including metals, polymers, and ceramics. 1.18.1.1.1 Classification of coupled physical fields in MEMS

In MEMS, certain measurements and outputs are detected in the bulk of MEMS devices, and others only on their surface. The same is true of coupled physical phenomena. For example, sensor effects are modeled by equations that are either coupled through domain effects, where a number of physical models are active within the same region and influence each other, or through common boundaries for the case where uncoupled models operate in nonoverlapping adjoining regions. This is relevant to a large number of physical models. We seek to clarify the basic operating principles of MEMS and to characterize the influence of coupling effects in MEMS devices. This classification is not the only possible one but provides an overview of the coupling effects in MEMS. 1.18.1.1.2 coupling

Domain-type multiphysical

The electro-thermo-mechanical interaction is a commonly used coupling effect in MEMS. This coupling effect is also called domain coupling because the physical phenomena of coupling occurs completely in the material domain. For example, the energy of an electrical current passes through a material and is partially dissipated into heat. At the same time, the

gradient of temperature causes thermal stress and results in structural deformation. This is a natural sequential coupling effect, which can be observed in daily life. In fact, there are fully coupled effects between electrical, thermal, and mechanical signals due to the nature of associated materials. If the independent variables, such as stress, electrical field, and temperature, and the dependent variables, such as strain, electrical displacement, and entropy, are chosen, then the coupling effects can occur between strain–stress as the elastic compliance, strain–electric field as the inverse piezoelectric effect, strain–temperature as thermal expansion, electrical displacement–stress as piezoelectric effect, electrical displacement–electrical field as permittivity effect, electrical displacement–temperature as pyroelectric effect, entropy–stress as piezocaloric effect, entropy–electrical field as electrocaloric effect, and entropy–temperature as heat capacity effect (Nathan and Baltes 1999). On the one hand, it is possible to observe how the independent variables are coupled with each other, based on these coupling effects. On the other hand, a certain coupling may have very small amplitude when compared with other couplings. For example, the value of elastic coefficients may change only in a very small interval (<1%) when retrieved under constant entropy or isothermal conditions. In this case, the coupling effect can be decoupled and the accuracy of the numerical solution will still agree well with measurements.

1.18.1.1.3 coupling

Boundary-type multiphysical

The electrostatic microactuation of a parallel plate structure is a widely used technique in MEMS. One electrode is typically a movable beam while the other electrode is fixed and separated by an air gap, which is within the micrometer range. This electromechanical interaction is a typical boundary coupling effect since the electrically induced mechanical loading is significant with regard to the electrical and mechanical interface. It is clear from the Poisson equation that the electrical force is perpendicular to the material boundary. Then the mechanical deformation of the structure alters the distribution of the electrical potential, thus changing the distribution of the surface charge density and hence the electrostatic forces acting on the structure. The self-consistent solution for both the mechanical deformation and the electrical potential distribution is coupled through the

Multiphysics and Multiscale Simulation

interface of two energy subdomains. It is especially tricky to compute self-consistency when the pull-in in electrostatic actuation occurs.

1.18.1.2

Coupled FEM Simulation

The primary goal of using numerical simulation for MEMS development is to enable and to accelerate the design compromise. It is also an important aid in improving the understanding of some physical processes, given that the model is correct. Of course, simulation is not a replacement for experimental evidence, and it cannot eliminate a hardware development cycle. It also does not make sense to pursue modeling without measured material property data. This necessity cannot be overemphasized. 1.18.1.2.1 field

PDE expression of physical

Physical effects can be described by partial differential equations (PDEs). A PDE relates an unknown scalar or a vector-valued function u[xi] to its partial derivatives for space and time: "

# qu qk u F x1 ; x2 ; . . . ; xn ; u; ; . . . ; k1 ;... ¼ 0 qx1 qx1    qxnkn

½1

The order of a PDE is given by the derivative of the highest order k present in eqn [1], and it is referred to as quasi-linear when all the highest order derivatives appear in a linear relationship. For most of the commonly used second-order PDEs, the characteristic at a point xi is determined by the coefficients of the second derivative term. We speak of elliptic equations when all eigenvalues of the coefficient matrix have the same sign, of hyperbolic equations when the eigenvalues have differing signs, and of parabolic equations when the matrix is singular. Only by appropriately specified boundary conditions, such as Dirichlet, Neumann, and Robin boundary conditions, a unique solution of the PDE is defined (Morton and Mayers 1994). 1.18.1.2.2 Numerical simulation of single PDE using FEM

In most cases, analytical solutions to the partial differential questions are possible only for a linear PDE with certain simple geometries and boundary conditions. For arbitrarily shaped geometries there can be no closed-form solutions. Then a PDE has to be solved numerically using suitable approximations. In this chapter we use the standard weighted

541

Galerkin finite element method (FEM) to discretize a PDE and obtain approximate solutions. In the FEM, the choice of weighting and approximation functions in the discretized weak form of the corresponding PDE has an important influence on the solution. In general, the FEM procedure includes deriving a PDE weak form, choosing an approximation function, meshing the computational domain, assembling the stiffness matrix and load vector, adding boundary conditions, and solving a linear algebraic equation. For more detailed discussions about the accuracy and the efficiency of FEM solutions, see Zienkiewicz et al. (2005). Here we merely discuss some key points that are very important to obtain reasonable numerical solutions in various MEMS field simulations. Mechanical analysis involves solving the Lame´ equation with the given displacement and force boundary conditions. In deformation analysis, it is usual to choose a linear finite element to discretize the computation domain in order to balance the computational cost and accuracy. This is especially important for 3D simulations. If only the numerical solution of the displacement is needed, then the linear finite element is sufficient in most cases, whereas this is not true anymore if the strain or stress distribution needs to be checked. For a 2D beam that is rigidly fixed at one end and is loaded by a force at the other end, the numerical solution of stress with linear finite element shows a major difference (generally >15%) in comparison to the analytical solution. The discrepancy cannot be compensated using more refined meshes. This phenomenon is called self-locking in beam analysis. The reason for this excessively stiff behavior can be attributed to the missing approximate basis function describing the bending terms using a linear finite element. Similar problems can also be met when solving thin plate and shell problems that have large rotations and small strain. One of the simplest remedies is to use a high-order finite element. However, the computational cost of the high-order finite element is much more than that of a linear finite element. An efficient method to obtain an accurate stress solution is to use a nonconforming element. The nonconforming quadratic complement element approximates the solution by the standard linear interpolation function and a nonconforming quadratic function. The nonconforming degrees of freedom (dofs) are condensed at the element level so that dofs of the element stiffness matrix are the same as for the linear element after static condensing. The nonconforming

542 Multiphysics and Multiscale Simulation

element improves the accuracy of stress distribution significantly, and there is no significant increase in the computational cost when compared with the linear element. Electrostatic analysis involves solving the Poisson equation to compute the distribution of the electrostatic potential and electric field in a given domain. The accurate calculation of electrostatic fields is very important in order to simulate microdevice behavior, such as capacitive sensors and electrostatic microactuators. Generally, the Lagrange element, which has potential dofs only on the element nodes, is used. The electric field is calculated through the numerical solution of the potential. In this case, the accuracy of electric field is not as good as the potential. When the solution of the Poisson equation has a singular character caused by a nonsmooth boundary or an inhomogeneous material distribution, the accuracy of the electric field is not good. Hermite-type elements can be used in which both the potential and the electric field are dofs. However, the Hermite type of high-order element is sometimes awkward to use when solving the electrostatic problem, which includes material interfaces. On the nonsmooth points, the distribution of the potential is continuous but the electric field is not continuous. So the numerical solution with standard Hermite-type element is oversmoothed compared with the exact solution. In order to find a remedy for the oversmooth problem of Hermite-type high-order elements, releasing of derivative dof is required. It is well known that on the material interface the potential and electric displacement is continuous, but the electric field is not continuous. The jump of the electric field depends on the ratio of material dielectric coefficients. So only one unknown of electric on the interface should be used when using the standard discretization method. One simple remedy is to rescale the element stiffness matrix before assembling into the global stiffness matrix, by the relationship of the electric displacement. The transformation is completed on the element level and the sum of dofs does not increase when compared with the standard method. Like the electrostatic problem, there is also a material interface problem when solving the electromagnetic problem based on the Maxwell equations. Based on a standard Lagrange element, the electric and magnetic fields and fluxes are determined indirectly by solving a scalar and vector potential. At the same time, the curl of electric field and the divergence of magnetic flux should be zero. It is well known that the tangential components of electric

field and magnetic field are continuous across the material interface, and the normal components of magnetic flux and electric displacement are continuous across the material interface. The difficulty in the modeling of field strengths caused by discontinuities in the properties of the media can be overcome by using edge elements. An edge element has the property of ensuring the continuity of tangential field components across an interface between different media, while leaving the normal field components free to jump across such interfaces. Hence, edge elements can be used for modeling an electromagnetic field along the interfaces between two different media. It can provide sufficient continuity for the electromagnetic field with a conforming-type element. In acoustics, sound pressure waves are often sinusoidal. Associated with a sound pressure wave is a flow of energy. The acoustic wave is extremely sensitive to changes in material properties, such as mass, elasticity, and viscosity of the propagation medium. Similar numerical discretization ideas as used for the electrostatic and electromagnetic problem can be used to deal with the material problem. Meanwhile, there is another problem regarding acoustic simulation with the standard Lagrange element, which uses polynomials as approximate basis. Since the basis solution of the acoustic problem can be expressed using superimposed sine functions, fine enough number of finite elements should be used to interpolate the high-frequency wave solution. Coarse meshes cannot reflect a sound wave. The basic reason for this problem is the limited approximation accuracy of the polynomial function for the high-frequency waves. Bubble elements have shape functions, which are zero at the boundaries of the mesh element, and have a sine or a cosine function in the middle of the mesh element. Combined with a linear basis, the bases used in the bubble elements are more suitable to simulate wave phenomena with coarse meshes when compared with the standard Lagrange element. Fluid problems are also a challenge for the FEM. The challenges arise from two aspects. The first is the consistent interpolation of the velocity and the pressure, and the second is the stable discretization of a pure convection operator. For the first problem, a general rule to obtain a consistent discretization is that the number of pressure dofs should never exceed the number of velocity dofs. When solving incompressible Navier–Stokes equations (NSEs), this demand should be valid independently of the number of elements used to discretize the computational

Multiphysics and Multiscale Simulation

domain. In order to satisfy this criterion, a generally accepted rule is that the order of approximation of the pressure must be one lower than the order of approximation of the velocity. So if the velocity is approximated by a linear polynomial, the pressure is approximated by a constant in each element and so on. For the second problem, nonstandard Galerkin FEM, such as streamline upwind Petrov Galerkin (SUPG) method, can be used to obtain a solution without nonphysical oscillations. 1.18.1.2.3 Numerical simulation of coupled PDEs using FEM

simulation. However, the coupling between the different programs has to be performed via external files. This generates a dependency on the data formats of the different programs. The applicable solution strategies for coupled simulations are limited to decoupled iteration schemes, because the coupling Jacobian matrix is not available. At the same time, the sequentially coupled method is not straightforward when extended to more than two fields, because no general staggering iteration can guarantee the convergence for a variety of physical coupling characteristics in MEMS.

Until now we have merely discussed the numerical simulation of isolated physical effects that are described by a single PDE, but neglecting all coupling effects with other state variables. However, MEMS devices achieve their functionality through coupling of different physical effects. In the following, we introduce coupling simulation based on boundary and domain coupling effects. The treatment is far from being complete but intends to illustrate some possibilities to accurately simulate coupling mechanisms in MEMS.

1.18.1.2.5

1.18.1.2.4

1.18.1.2.6

Sequentially coupled method The sequentially coupled method is typically used in boundary-coupled problems. It is also called partitioned analysis of coupled problems (Felippa 1999), and it was first used to solve the fluid–structure interaction problem. The numerical solution of the structural deformation and fluidic field is implemented separately, and the coupling effect is transmitted through the fluid–structure interface. This means that different programs can be called as black boxes separately, and the finite element meshes need not coincide at the interfaces (http://www.ansys.com/). The decomposition is driven by physical or computational considerations for sequentially coupled methods. Another case using the sequentially coupled method is the loosely coupled problem. For the electrothermal problem, a loosely coupled effect exists if the electric current is not influenced by the change in the temperature. The electric current is first calculated and then solved for the temperature distribution. The advantage of using several programs lies in exploiting the know-how contained in each program, enabling research effort to concentrate on the development of coupling mechanisms and solution strategies, especially since there are many commercially available general-purpose modeling tools, which are extensively used for MEMS

543

Fully coupled method In practical MEMS simulation, many coupled problems are fully coupled. It is natural that one discretizes the coupled PDEs directly. In fully coupled discretization, a unified mesh is used in the whole computational domain and each nodal point has the same dofs, which includes unknowns from all the coupled PDEs. With the standard FEM, the offdiagonal submatrices of the global discretization matrices describe the coupling effect and associated relation (http://www.femlab.com). Linearization of nonlinear PDEs For most of the coupled problems in MEMS, the nonlinearity comes from coefficients or material properties that contain a dependent variable or from certain terms in PDEs that have a dependent variable of degree 2 or higher order. Nonlinear PDE sets should contain at least one nonlinear PDE. In order to solve a system of nonlinear equations, an iterative procedure is necessary. In general a suitable initial estimation is given, and then linearized nonlinear equations are solved based on a previous solution until a converged solution is obtained. Generally there are two ways to derive the iterative method. First, one can discretize the nonlinear PDEs and then solve the nonlinear equation system. An alternative is to linearize the nonlinear PDEs and then to discretize the resulting linear equations (Belytschko et al. 2000, Gockenbach 2002). Sometimes both approaches are identical. The second approach is conceptually easier than the first one, and therefore it is widely used in MEMS simulation. When linearizing the nonlinear PDEs, the nonlinear terms are treated using the Taylor expansion, leaving only the linear terms. This is called the Newton method. Alternative linearization is constructed by the so-called Picard iteration method.

544 Multiphysics and Multiscale Simulation

The nonlinear terms are linearized in which all the nonlinear unknown terms are replaced with an iterative solution from the last time step. An important question with respect to these two iterative methods is how to give a good initial estimate. It is well known that the Newton method converges fast (quadratically) when the initial value is in the neighborhood of the converged solution. However, if the distance between the initial value and the converged solution is too large, the Newton method may converge slowly or may even diverge. The Picard iteration seems to have a larger convergence region, which means that it does not need as accurate an initial estimate. However, Picard iterations converge much slower than the Newton method (linearly). A possible strategy to balance the stable and fast convergence is to use Picard iterations first with an initial value and then use Newton iterations when the iteration residual is smaller than certain value. In strong nonlinear problems, the convergence of nonlinear iterations may be very sensitive to the quality of the initial value even though one is using combined Picard and Newton linearization. In this case, the robustness of the iterative property can be improved using the continuation method (Allgower and Georg 2003). The continuation method defines an easy problem and the original problem that one actually wishes to solve. The solution of the easy problem is gradually transformed to the solution of the hard problem by tracing this path. The path can be defined by introducing an additional scalar parameter that can weaken the nonlinearity of the original problem.

is also taken at the new time level k þ 1, the method is called implicit. In the explicit method, one still has to solve a system of equations if the matrix M is not lumped. In the case of a lumped M matrix, the solution implies only the inversion of a diagonal matrix and the explicit method is relatively cheap. A disadvantage of the explicit method is that the time-step length must be restricted based on the Courant– Friedrichs–Levy (CFL) condition in order to obtain a stable solution. For a stiff ODE, the size of the time step in the explicit method can be so small that the overall computation is not efficient anymore. For the implicit method, it may be unconditionally stable for well-defined problems. The accuracy requirement is the main reason to restrict the size of the time step. In coupled problems, which include both timedependent PDEs and steady PDEs, the matrix M becomes singular. Differential algebraic equations (DAE) have to be solved after the space discretization of the coupled PDEs. This is also the case where there is a big difference in magnitude of the solutions of the coupled time-dependent problems so that the matrix M is almost singular. One has to note that DAEs are different from ODEs, although a general DAE can be transformed into an ODE theoretically by variable substitution. The variable-order, variable-step size backward differentiation formula (BDF) is an efficient method to solve a DAE problem. Solving DAEs is also an important step in the coupled optimization method (COM) for fully coupled physical fields discussed later on in this text (Ascher and Petzold 1998). 1.18.1.3

1.18.1.2.7

Time-dependent solvers

The discretization of the time-dependent PDE results in a system of ordinary differential equations (ODEs). In order to solve this system of equations any classical method for solving ODEs can be used. In general one can distinguish between explicit and implicit methods, and between one-step and multistep methods. In the FEM, the time derivative can be replaced by a forward difference discretization (or method of lines): M

  Ukþ1 – Uk þ KU ¼ F t

½2

where k denotes the present time-level and k þ 1 the next time-level. A method is called explicit if the term KU is taken only at the time level k. When KU

Mesh Adaptation

In MEMS, an efficient coupled field simulation procedure depends on many detailed numerical strategies. In the FEM, a good mesh is a very important issue, which influences the efficiency and the accuracy of the numerical solution (Thompson et al. 1985). When users have limited experience in the properties of the coupled physical fields, it is nearly impossible to decide upon a suitable mesh to discretize the computational domain in advance. In this case, mesh adaptation can bridge the gap. Basically, dynamic mesh adaption can be implemented in two general ways, i.e., adaptive mesh refinement (AMR) (Costa and Alves 2003) and adaptive moving mesh (AMM) (also called the r-adaptive method) (Cao et al. 2003). The mesh is locally remeshed in the AMR method, thereby capturing the important part of the

Multiphysics and Multiscale Simulation

solution. This method is usually accompanied by a relatively complicated data structure updating procedure. For the AMM method, numerical solution of the moving mesh is based on a moving mesh PDE. The nodal points are moved throughout the region to best approximate the solution globally. The AMR method is suitable for domain-type coupled problems in MEMS. When using the FEM to solve a PDE, there are two major AMR methods available, the h-method and the p-method. For the h-method, the mesh is refined or coarsened by adding or deleting mesh grids according to a local error analysis. For the p-method, the mesh is fixed and the order of a finite element is varied to implement a highorder approximation locally. There have been extensive studies on the h-method and p-method, both of which have shown good efficiency for both steadystate and time-dependent problems. The key point of the AMR method is to define a remeshing strategy. This guarantees the quality of the finite element solution and smoothens the expressions of the material domain, while keeping a reasonable number of design variables. Generally, an error estimator is used to determine the adaptation strategy, since the exact solution of the direct problem is unknown. For a single physical field problem, such as a linear elastic mechanical problem, the error estimator by Zienkiewicz and Zhu can be expressed as a least square approximation between the approximated nodal stress and the smoothened stress in the Gaussian points for all finite elements. Since the solution range for different physical fields can have difference in magnitude of several orders, the error estimators for coupled problem cannot be based on the numerical solution directly. One strategy is scaling the numerical solutions based on the first calculation. User can also provide an initial guess based on the physical meaning of the unknowns of the coupled PDEs. In the AMM method, the mesh topology is kept unchanged but the mesh grid points are concentrated inside the physical domain where the solution of the physical problem has large gradient or large solution variations (Beckett et al. 2001). Compared with the AMR method, the AMM method does not require complicated data structures for mesh refinement. Furthermore, the AMM method does not need to interpolate the numerical solution between different levels of mesh refinement. For the boundary coupled problem, the AMM method coincides with the arbitrary Lagrangian–Eulerian (ALE) method, which is popularly used in the fluid–structure interface problem. In a Lagrangian formulation, the update of a

545

mesh follows the change of the computational domain. The advantage of a Lagrangian expression for a moving mesh is that it can retain a sharp material interface when the mesh is deformed properly. However, one of the challenges of using a Lagrangian-type formulation is that the mesh can become severely distorted or can overlap itself. So after the computational domain is updated, the mesh points inside the domain should be smoothly adjusted using the AMM method. Of course, it can also be possible to use both the AMR and the AMM methods during simulation, e.g., when the domain and boundary couplings exist on different subdomains in one simulation.

1.18.1.4

Multiphysics Design Optimization

Microsystem components are hard to design, because manufacturing technology significantly constrains our creative possibilities, and the sensitive regions where sensor and actuator effects are active are usually not neatly lumped but extend over large portions of a device, making system design approaches hard to apply. Currently, very few automated design methods are available for MEMS (Sigmund 2001a, b, Yin and Ananthasuresh 2002). To establish such methods, a number of numerical components must be in place so that together they can form an integrated design and optimization tool. Here we introduce an integrated structural topology optimization method for multiphysics in MEMS. 1.18.1.4.1

Structural layout optimization The optimization of structural layout potentially has a great effect on the performance of a structure through parameters such as weight, volume, stress distribution and dynamic response, etc. Generally, structural layout optimization can be classified into three categories: size, shape, and topology optimization. For size optimization, only simple geometrical parameters, such as the thickness of a plate or the length of a bar, are chosen as the design variable. The shape of the structure is fixed throughout the optimization procedure. For shape optimization, the design strategy to improve the performance of the structure is to adjust the shapes of structural boundaries. The default hypothesis is that the topology of the structure or the connections of the structure remain unchanged. Therefore, shape optimization

546 Multiphysics and Multiscale Simulation

cannot answer topological questions. The topology optimization method is a suitable and an efficient methodology in this case. When compared with size and shape optimization, the topology optimization method is an advanced method, which can obtain an optimal topology, shape, and size at the same time (Rozvany 2001). The continuous topology optimization method has been used to design the layout of structures since the development of the homogenization method (Bendsoe and Kikuchi 1988, Bendsoe and Sigmund 2003). Based on the ground structure method, an initial design domain is discretized by finite element meshing. According to this fixed meshing, an optimal topology can be obtained without guessing an initial topology of the structure. It should be noted that the topology optimization method is equivalent to an inverse problem (first kind is Fredholm integral equation), which considers the material parameter distribution for coupled problems. There is no general theory to obtain a stable solution for this kind of inverse problem. Numerical regularization is generally used in most cases. 1.18.1.4.2 Simultaneous simulation and optimization method

To date, topology optimization has covered several applications, such as functional materials, micro sensors, and compliant actuators, in MEMS design. As an example, let us consider electrothermal microactuator, which is actuated through resistive heating and thereby thermal expansion effects, modeling of electric, thermal, and elastic fields, and then calculates the sensitivity from these three fields to implement the structural optimization iteratively. Numerical instabilities during the optimization procedure, checkerboard patterns and mesh dependency, have been overcome based on the numerical filters or the homogenization method. Geometrically nonlinear deformation, which generally occurs as a snap through or buckling phenomenon in an electrostatic actuated switch, is also a research topic in nonlinear topology optimization. However, further application of structural topology optimization in MEMS is impeded by the complexity of the optimization algorithm itself and the difficulty in finding a unified expression for the optimization procedure for diverse coupled problems in MEMS. For example, there are many surface-type coupling effects in MEMS, such as the structure actuated using electrostatic force. It is known that the deformation of the structure also influences the magnitude of the electrostatic force

and vice versa. During structural optimization the electrostatic force is design-dependent. It means that the force is simultaneously changed following the change of the structural boundary in the material frame during the procedure of optimization. Here an efficient and accurate expression of these forces is a key point and has also been a challenging problem to date. Traditionally, topology optimization is based on the sequentially COM. Specified physical problems can be discretized and solved first, and then sensitivity based on the discretized numerical solutions can be calculated. This is a discretize-then-optimize approach. Another optimization strategy is the COM. One can establish coupled physical equations and an optimization equation at the original PDE level, and then discretize these coupled PDEs together by a numerical method. This is an optimize-then-discretize approach. Here we describe the first as the sequential optimization method (SOM) and the second as the COM (Arora and Wang 2005). The COM has several advantages when compared with the SOM. First, in the SOM, it is obvious that the effect of an optimization iteration is strongly dependent on the accuracy of their discretized physical problems, because the discretized solution of physical problems is often a noisy or a discontinuous approximation of the initial problem. Numerical errors pollute the quality of sensitivity and influence the convergence. For the COM, the sensitivity, which connects the physical equations and the optimization equation, is calculated by the true physical relationship without any numerical discretization error. So it is often less noisy and more accurate than the discretized optimization method. Second, in the SOM, the state variables are treated as implicit functions of the design variable. Then the design sensitivity analysis should include the derivative terms of the state variables and the design variables. In the COM, the state variables are set independent of the design variables, i.e., the state and design variables are defined in a composite vector. One of the distinguishing advantages of the method is that the derivation of the sensitivity is relatively straightforward. This is especially true when one wants to implement structural optimization in which the state variables are controlled by multiphysical problems. For example, equilibrium equations could come from mechanics, electrostatics, and thermoelectric effects. The state variables include displacement, electrostatic field, thermal stress, etc. The coupling effects include surface coupling of

Multiphysics and Multiscale Simulation

mechanism and electrostatics, and domain coupling of mechanics and thermoelectric effects. In this case, the sensitivity analysis of the SOM is extremely complicated, and the expression of the COM is much simpler for coupled problems. Then one question is how to establish a relationship between the decoupled state and design variables. One of the possibilities is to derive a corresponding optimization equation based on the specified objective and to solve the fully coupled equilibrium equations and optimization equations simultaneously. Typically, the element density-type method, such as the solid isotropic microstructures with penalization (SIMP) method, is used to implement structural topology optimization (Allaire et al. 2004). However, the expression of the structural layout is nonsmooth because piecewise constant density values are used to express the material domain. Even though this does not cause serious numerical problems for the structural optimization problems, nonsmooth expression of material leads to unexpected results for other types of structural optimization, such as thermal problems (Yoon and Kim 2006). One of the simplest remedies to overcome the problem of nonsmooth expression of a material is the nodal density-type method, where material density is interpolated smoothly with nodal density values. In MEMS simulation, the actuation force is normally either a body force or a surface force. During structural topology optimization, these two kinds of forces are referred to as designdependent loads, which mean that the load vectors change with the evolution of the structural layout. With either the element- or the nodal density-type method, there is no explicitly defined material boundary. So the implementation of a design-dependent load is not straightforward (Du and Olhoff 2004, Hammer and Olhoff 2000). Basically, the level set method is shape optimization method. Because of its ability to merge holes during the optimization, it has recently become popular in structural topology optimization. When compared with the density-type method, one of the distinguishing advantages of the level set method is the material boundary with zerolevel set contours; therefore the design-dependent load is easily implemented during the optimization procedure. 1.18.1.4.3 equation

Derivation of the optimization

In this section we introduce the use of the level set method to transform a typical optimization procedure to a coupled PDE expression. We first introduce how to

547

establish a coupled simulation and optimization for a single mechanical and electrostatic field problem, so that the structural optimization procedure can be solved as a pseudo-time-dependent coupled PDE problem. Then we extend this procedure to the electromechanical-coupled optimization problem in MEMS. A typical structural topology optimization problem for 2D compliance minimization can be written as: Min :

c½; " ¼

Z

1 E½"T D" d 2



s:t:

Z

r?ðE½Þ ¼ f

½3

H ½d ¼ Vol



where  is the design domain, u is the displacement field calculated using the linear elastic equilibrium equation, " is the strain tensor,  is the stress tensor, E is the design variable and is defined by the level set surface  as E[] ¼ H[] þ (1  H[])Emin, where Emin is a small number, such as 0.001, to avoid singularity of the stiffness matrix, D is the elasticity matrix in which the material Young’s modulus is E0, and Vol is the material volume constraint in the deign domain. The level set function [x] is an implicit function for a given domain  with smooth boundary, which satisfies the following (Osher and Fedkiw 2002): ½x > 0;

xPþ ðmaterialÞ

½x ¼ 0;

xPq ðboundaryÞ

½x < 0;

½4

–

xP ðholeÞ

The area constraint condition, which is typically used in compliance minimization design, is presented. To derive the level set equation, which moves the material boundaries, we combine the objective function and constraints together by using a Lagrangian formulation, and then derive a corresponding Euler–Lagrangian equation (Nocedal and Wright 1999). Using the Lagrangian multiplier , we can rewrite problem [3] as J ½"; ;  ¼

Z  



  1 Vol T d ½5 E½" D" þ  H ½ – 2 Vol

where Vol is the area of the entire design domain. From classical calculus, we know that the extrema of functional [5] is attained at the position where J9 ¼ 0. To obtain the level set equation, which uses a level set surface to express the structural topology implicitly, we need to calculate the variation of the level set surface. Only the normal velocity of the material boundary influences the change of the shape. The

548 Multiphysics and Multiscale Simulation

tangential velocity will not influence the deformation of the geometry, merely its parameterization (Sapiro 2001). We assume that the zero level set contour  ¼ 0 moves only in the normal direction. For an infinitesimal variation l, which is along the normal direction n* , a new zero-level set is 9 and we have j¼0 ¼ 9 ½x – ½x ¼ ½x þ l ? n*  – ½x ¼ r ? ðl ? n* Þ ¼ jrj

½6

where the normal vector n* ¼ r=jrj: For problem [5], the variation of the level set surface on the material boundary can be expressed as  Z  1 ð1 – Emin Þ"T D" þ  ½jrjl d  J ¼ 2

½7



The corresponding Euler–Lagrangian equation at the extreme value point is   1 ð1 – Emin Þ"T D" þ  ½jrj ¼ 0 2

½8

In most cases, it is impossible to solve eqn [8] directly. One general technique is to solve set equation numerically   q 1 – ð1 – Emin Þ"T D" þ  ½jrj ¼ 0 qt 2

½9

The auxiliary variable t is used as pseudo-time. The method is also referred to as gradient descent flow. In the case of a well-posed optimization problem, the solution converges to a local minimum based on different initial values. In this way the 2D compliance minimization topology optimization problem can be transformed into the following coupled plane stress equation and time-dependent reaction diffusion equation r?ðE½Þ ¼ f   q 1 – ð1 – Emin Þ"T D" þ  ½jrj ¼   qt 2



½ ¼ mat H ½ þ ð1 – H ½Þair

where air is the dielectric permittivity of the air and mat is the dielectric permittivity of the design material. The equilibrium equation, which controls the electrostatic problem, is the Poisson equation. Following a similar procedure as for the mechanical problem, a maximization of the electrostatic energy problem can be expressed as C½; rV  ¼

Max :

1 ½rV T rV d 2



½12



A fully coupled PDE expression of topology optimization for the electric potential energy problem with area constraint can be similarly expressed as r?ð½rV Þ ¼ 0   q 1 – ðmat – air ÞrV T rV þ ½jrj ¼   qt 2

½13

where

½10 ¼

½11

Z

s:t: r?ð½rV Þ ¼ 0 Z H ½d ¼ Vol

Z 

with a suitable initial condition 0. Here an artificial diffusion term   is added to overcome numerical oscillations during the evolution of the level set surface. The coefficient  is proportional to the mesh size and the reaction term. The Lagrangian multiplier is expressed as  Z  1 ð1 – Emin Þ"T D" 2 ½jrjd 2 Z ¼  2 ½jrjd

Here we call the first equation in eqn [10] as the state equation for the direct problem and the second equation as the optimization equation for the structural optimization problem. For the electrostatic problem, one of the most commonly used objective functions is the stored potential energy inside a design domain. The electrical energy can be expressed as []rVTrV/2 where V is the electrical potential, rV is the electric field vector,  is the dielectric permittivity, which is defined by a level set surface  as



 1 ðmat – air ÞrV T rV 2 ½jrjd 2 Z 2 ½jrjd

½14



Hence as for the mechanical problem, an artificial diffusion term is added to the right-hand side of the second equation in eqn [13]. During the procedure, the amplitude and the direction of the electrostatic force will change with the update of the topology of the structure and its deformation. Hence we have a coupled optimization problem. There are two different methods to solve this surface coupling problem: the sequentially coupled solver and the directly coupled solver. To

Multiphysics and Multiscale Simulation

implement structural optimization with the coupled method, the directly coupled solver is used to calculate the structural deformation. The design domain is shown in Figure 1. There are two subdomains. Subdomain 1, the upper strip of the beam marked in black, has specified constant material properties. Subdomain 2 is the optimization domain and is white in Figure 1(a), 1(c), and 1(e). Here we discuss three different load cases. In the first case, a fixed load vector is added on the middle point of the upper surface of the whole domain. In the second case, the electrostatic force acts on the structure. In the last case, both the fixed load and the electrostatic force act together on the whole structure. The objective of structural optimization is still compliance minimization with area constraint. The material inside both subdomains is aluminum with Young’s modulus of 70 GPa and a Poisson ratio of 0.33. In the first case, we use the compliance minimization model in eqn [3]. In addition, we consider the geometrically large deformation for the elastic equilibrium equation. So the strain tensor " has a nonlinear relationship with the displacement vector. Since we use the COM (Liu et al. 2005) in which the design variable  and the displacement vector are independent of each other, the optimization equation is still the same as the second equation in eqn [10]. The first equation in eqn [10] is solved using a nonlinear solver for large deformation analysis. The optimized structural topology is shown in Figure 1(b).

549

In the second load case, a consistent solution for the electromechanical problem is needed. To coincide with the fully coupled optimization method, we choose a direct coupled solver where the deformed mesh caused by the structural deformation is simulated by the ALE method ( Jakobsen et al. 2001). The couple ALE method can be implemented using the following coupled equations: r ? ð½rVX Þ ¼ 0 r ? ðE½Þ ¼ Fele r ? ðGrX Þ ¼ 0

½15

where X represents the deformed nodal position because of structural deformation, G is the mesh monitor function, which controls a local property of the deformed mesh, VX means that the value of the electrical potential is calculated on the deformed mesh, and ½ ¼ max H ½ þ ð1 – H ½Þair

where max is a large number to approximate the distribution of the electrical potential inside the conductor. In this example max ¼ 100. Using the moving mesh equation (the third equation in eqn [15]), all of the displacement, force, and electrical potential boundary conditions can be implemented in the original coordinates. This method can deal with the electromechanical coupling directly before pull-in occurs. For the optimization examples we presented here, the fixed–fixed end beam in subdomain 1 is used (b)

(a) 10 µm 1 µm

Ffixed

60 µm

(c)

V = 200 V

(d)

V=0 V

(e)

V = 200 V

Ffixed

(f)

V=0 V Figure 1 Compliance-minimized topology optimization of a fixed–fixed end beam and three load cases. (a) and (b) Initial design domain and optimized topology with fixed load vector; (c) and (d) initial design domain and optimized topology with electrostatic force; (e) and (f) initial design domain and optimized topology with both the fixed load vector and the electrostatic force. The beam in (a) has a stiffness of 0.019 N m1 at its center point, and for the situation in (c), it has a theoretical pull-in voltage of 1700 V. In all simulations F was set to 10 N. (Source: Lui Z, Korvink J G, Reed M L, Multiphysics for structural topology optimization. Sensor Letters. 4, 191–199.)

550 Multiphysics and Multiscale Simulation

to avoid pull-in occurring during the optimization procedure with reasonable electrical potential conditions. The coupled optimization equations are: r?ð½rVX Þ ¼ 0 r?ðE½Þ ¼ Fele r?ðrX Þ ¼ 0  q 1 T – ð1 – Emin Þ" D" þ  ½jrj ¼   qt 2

½16



Here the mesh monitor function is the unit function so that the mesh is smoothly modified on the whole design domain. On each mesh node, there are six dofs, two displacement dofs, one electrical potential dof, two moving mesh position dofs, and one level set surface dof. However, the optimization equations remain unchanged despite the complicated coupled forward problems. The optimized result is shown in Figure 1(d). In the last case, both the fixed load vector and the electrostatic force act on the structure. So this is a multiple load optimization problem. Based on the analysis of the first two cases, the optimization equations for these two loads are the following: r?ðE½1 Þ ¼ Ffixed r?ð½rVX Þ ¼ 0 r?ðE½2 Þ ¼ Fele r?ðrX Þ ¼ 0    T  q 1 T – ð1 – Emin Þ "1 D"1 þ "2 D"2 þ  qt 2

½17

 ½jrj ¼  

where "1 and "2 are the strain vectors corresponding to two load cases. The optimized result is shown in Figure 1(f).

1.18.2 Multiscale Simulation Modeling and simulation of the behavior of materials in various applications of modern microsystems is a challenging task. On the other hand, high costs or restricted availability of samples for analysis requires modern microfluidic systems to operate in the picoliter regime and below. The ultimate goal is to make available simulation tools that emerge from a theoretical modeling that ranges from quantum physics of atomic-scale phenomena, covering up to continuum descriptions of macroscopic behavior, and going even further to system behavior. Such engineering tools, supporting design processes in microsystem technology, are still missing.

In microsystems, the presence of length scales of several orders of magnitude requires the coupling of different tools appropriate for simulation at the respective spatial scale. The use of advanced computer-aided design (CAD) tools promises reduction in the extent of physical testing necessary to prototype a device. Activities to incorporate various physical models at different length scales (e.g., Nielaba et al. 2002) exist. The simulation of complex processes in liquids and the modeling of the material properties of liquids are tasks that require resolutions on the molecular scale and below. The computational effort at the atomistic level is demanding. Consequently, atomistic approaches should undergo a coarse graining procedure. Simulating the behavior of a reasonable amount of liquid with various components in a common framework with a reasonable computational effort – as it is mandatory for the application of microfluidic systems in a design process – can be realized only by superimposing microscopic molecular dynamics (MD) approaches on coarse-grained particle dynamics methods. Typically, the microscopic MD simulation is computationally more expensive than coarse-grained particle dynamics solvers and therefore MD must be restricted to very specific parts of the whole system where an atomistic resolution is desirable. It is even conceivable that a successive coupling to continuum Navier–Stokes models is necessary. To overcome the problem of the costly particle dynamics calculation being repeated in every iterative loop of the continuum solver, a microscale simulation featuring noniterative and explicit time marching must be coupled to an equivalent model acting at a much coarser length scale. This eventually leads to a feasible simulation process for engineering CAD applications in multiscale microfluidic systems. A promising technique to fulfill the above-mentioned requirements is dissipative particle dynamics (DPD), a particle dynamic simulation technique for simulations above the molecular length scale (Hoogerbrugge and Koelman 1992), which is described in Section 1.18.2.1. A second approach to fluid dynamics that couple the DPD model to more macroscopical environments is given by the lattice Boltzmann method (Wolf-Gladrow 2000). 1.18.2.1

Particle Dynamics Methods

Molecular dynamics is the most popular particle dynamics method, and it has a long-standing tradition. Its main ingredients are the interaction potentials of the atomistically resolved particles together with an algorithm to integrate the set of differential equations for this many-body problem. In essence it is an

Multiphysics and Multiscale Simulation

integration of Newton’s equation of motion. Much effort has been put into the development of efficient algorithms. On the other hand the interaction potentials can be derived from first principles and their form is crucial for the properties of the simulated material. Coarse-grained particle dynamics methods loose this direct connection to the first principles derivation of the interaction potentials. Instead, they must rely on either MD simulations or macroscopic measurements to adjust the respective parameters. Nevertheless, they benefit from the development of the efficient integration techniques for MD systems. 1.18.2.1.1 Coarse-graining microscopic methods vs refinement of continuum models

At a microscopic length scale on the order of interatomic distances in fluids or solids, MD represents an attractive simulation technique in computational fluid dynamics. However, at the length scale from 10 to 1000 nm a coarse graining of the molecular model due to rising computational costs is desirable. There are two ways to arrive at a particle description of microfluidic phenomena: coarse graining of MD and a direct Lagrangian discretization of the equations governing hydrodynamics. The former leads to the DPD variants discussed in Section 1.18.2.1.2, while the latter is exemplarily discussed here by the method of smoothed particle hydrodynamics (SPH) (Koumoutsakos 2005). SPH is a method for the simulation of hydrodynamic phenomena. Originally, SPH was developed to overcome hydrodynamic problems in astrophysics. Eventually it was adopted for the investigation of laboratory scale problems as well. It is a direct discretization of the macroscopic hydrodynamic equations on Lagrangian grid. The spatial derivatives are found by means of an interpolation formula. One of the biggest drawbacks of SPH is the fact that thermal fluctuations are not included in the equations of motion. For modeling flows at astronomical length scales, fluctuations most certainly do not play a role. In applying SPH to microfluidics where the dimensions are orders of magnitude smaller, fluctuations definitely do come into play. In a recent effort a generalization of SPH including thermal fluctuations has been proposed (Espanol and Revenga 2003) and termed smoothed dissipative particle dynamics (SDPD). 1.18.2.1.2

Dissipative particle dynamics

DPD, first introduced by Hoogerbrugge and Koelman (1992), is expected to be suitable for the

551

intermediate range of length scale. Its relation to MD as well as to continuum methods has been investigated (Flekkoy and Coveney 1999, Flekkoy et al. 2000). DPD is a mesoscopic simulation method capable of bridging the gap between atomistic and mesoscopic simulations (Groot and Warren 1997). In essence, DPD is a treatment of the dynamics of quasi-particles representing small sets of the liquid’s molecules by stochastic differential equations similar to a Langevin approach. It combines features from MD and lattice gas methods. Figure 2 provides a schematic representation of the context in which DPD has to be seen with respect to other simulation methods. Since its first introduction this method has been applied for the simulation of a wide range of phenomena, especially in material science. Thorough investigations have been carried out to understand the capabilities of DPD when applied to computational fluid dynamics. Recently, DPD has experienced an increasing interest of different research groups for its improvement in several aspects. Its algorithmic optimization was one focus, making it an appropriate method for application in engineering fields (Kauzlaric et al. 2005). The inclusion of energy and momentum conservation in the particle–particle interaction for the set of stochastic differential equations describing a DPD models is an important extension for heat flow applications. Recently, phase change models have been included, which is built on the energy-conserving DPD (DPDE) models. A solid– liquid phase transition and liquid–vapor coexistence can be modeled with DPD without turning down all the advantages of the method that come from using larger time steps than in MD. We will briefly describe the basic ingredients of DPD. Every dissipative particle i with its position ri and its momentum pi or its velocity vi is propagated according to Newton’s law with the force Fi acting on a particle. Thus the equations of motion is given as follows: r_ i ¼ vi  X R FCij þ FD v_ i ¼ Fi ¼ ij þ Fij

½18 ½19

j 6¼i

where the force has three main contributions due to other particles in the flow, namely a conservative part FC, a dissipative contribution FD, and a fluctuating force FR. The latter two are defined as FijD ¼ – ij w½rij ðeij ? vij Þeij

½20

FijR ¼ ij w1=2 ½rij eij ij ½t 

½21

552 Multiphysics and Multiscale Simulation

m Micro fluidics Lumped network models

Navier–Stokes equations (CFD, SPH ...)

Local balance Length scale

Ense m aver ble age

Spatial average

Global balance

Mo

me

Qunatum dynamics (Car-Parinello...) Weighted average

nts

Boltzmann transport equation

Qunatum dynamics (Car-Parinello...)

Wigner transport equation

Locality

Locality

Classical molecular dynamics

Mesoscopic models

Hamiltonian mechanics

Quantum nm mechanics

Figure 2 DPD model in the context of other simulation methods. BD, Brownian dynamics; CFD, computational fluid dynamics; DPD, dissipative particle dynamics; SPH, smoothed particle hydrodynamics. Deriving a mesoscopic technique from a microscopic one or a macroscopic technique from a mesoscopic one usually involves some sort of projection (averaging).

where  and  are the dissipation strength and the fluctuation amplitude, respectively. w½rij  is a weighting function of the distance rij of the two dissipative particles i and j, while eij is the unit vector pointing from particle i to particle j, and vij is the relative velocity of the two particles. ij ½t  describes the Gaussian white noise according to the following equation  ij ½t  ¼ 0

½22

 ij i½t  i9j 9 ðt 9Þ ¼ ðii9 jj 9 þ ij 9 ji9 Þ½t – t 9

½23

The modeling of material properties is done in DPD by adjusting the parameters entering the forces given in eqn [19]. These are  and  and form the conservative potential. An advantage in computational effort results only if FC is a soft repulsive potential. This allows for much larger time steps than in MD. On the other hand, the lack of an

attractive part in FC does not allow for free surface flow. Hence a density-dependent part must be added to the conservative potential accounting for the excess surface energy in order to describe liquid– gas interfaces. Since this is an extension of the simple DPD model to many-body behavior, it is termed as MDPD in the literature (e.g., Paganobarraga and Frenkel 2001). So far we treated only isothermal DPD even with its extension to MDPD. To include thermal properties of the material the internal energy of the individual dissipative particle and its interactions with the translational dofs should be taken into account. The internal energy as an example for the extension of the DPD to DPDE (e.g., Pastewka et al. 2006) devises a way for the introduction of various other dofs. Clearly, their equation of motion should be coupled to those of the existing dofs, and is a task for material modeling. With the extension to

Multiphysics and Multiscale Simulation

553

MDPDE, the phase description of phase transitions has become possible, which makes this method also very attractive in engineering applications when compared with direct numerical treatment of the hydrodynamic equations. 1.18.2.1.3

volume fraction

Applications of DPD

The area of applications for DPD is growing very rapidly. Its range for different simulation tasks is listed as follows. A discussion on mesoscopic dynamics of colloids can be found in Dwinzel et al. (2002), the treatment of binary fluids and immiscible fluids is reported in Coveney and Novik (1996) and Dwinzel and Yuen (2000), the rheology of dense colloidal suspensions can be found in Boek et al. (1997), and many more. The important step was the derivation of a thermodynamically conservative interacting force from a given free energy density, which allowed for the implementation of phase coexistence and the description of phase boundaries (Paganobarraga and Frenkel 2001). The implementation of DPDE was a very important development to make the model applicable for phase transitions (Boneyt Avalos and Mackie 1997). Newer developments overcome the limitations of DPDE with respect to the time steps in simulations (Pastewka et al. 2006), and thus allow for its application to the study of complex fluids. There are recent investigations on the applicability of DPD in the framework of micelle formation (Schulz et al. 2004), which makes it a promising candidate for the task of bridging the gap of length and time scales for the simulation of ionic liquids as well, even if there are still strong implications for the applicability of the method. In Pool and Bolhuis (2006), the failure in the prediction of a critical micelle concentration by an error of thirteen orders of magnitude was exposed impressively. We noted that there is a need for further developments to make DPD applicable in specific materials systems, especially when it comes to difficult rheological properties (Kauzlaric et al. 2003a, b, 2004). In Figure 3 the filling of a bending specimen mold was simulated. Experimental findings reported on the jetting behavior around the obstacle in the flow field. This was exactly reproduced by the DPD simulation,

0.450 0.362 0.275 0.188 0.100

Figure 4 Characteristic powder distribution in micropowder injection molding.

while a conventional simulation tool failed in the prediction of jetting. A very important information in powder injection molding is the ceramic or metal powder distribution after the molding process. A homogeneous mold mass is subject to demixing due to the shear rates occurring in the molding process. Figure 4 shows that by the extension of the DPD model a direct prediction of the powder concentration is possible.

1.18.2.2

Kinetic Models

Kinetic theory, as opposed to particle dynamics, starts from a different point of view for a system of many dofs. A state of N particles is given by its 3N coordinates ri and its 3N momenta pi to give 6N dofs. Instead of solving the coupled 6N equations of motion, as done in a particle dynamics approach, we rather regard the system as represented by a probability density fN ½r1 ; . . . ; r3N ; p1 ; . . . ; p3N . The quantity fN has the meaning of a correlation function for the 6N dofs. By knowing the interaction potentials of the N particles we can form the Hamilton function H ½r1 ; :::; r3N ; p1 ; :::; p3N  and write down the 6N Hamiltonian equations of motion as follows r_i ¼

H qpi

p_i ¼ –

Figure 3 Filling simulation of a bending specimen mold performed by dissipative particle dynamics (DPD).

H qri

½24 ½25

This is the starting point of particle dynamics methods, for example, MD, which proceed by solving eqns [24] and [25]. Instead of doing so we can write

554 Multiphysics and Multiscale Simulation

down the equation of motion for fN, which in this case in treated as an observable quantity, to give qfN ½ri ; pi ; t  qrj qfN ½ri ; pi ; t  qpj qfN ½ri ; pi ; t  þ þ ¼ 0 ½26 qt qrj qpj qt qt

For solving eqn [26] the information from eqns [24] and [25] is needed and thus there is no advantage. In kinetic models the interest is not in the full correlation function fN, but instead in the probability of finding one particle in phase space at a certain point f1 [r, p]. This function can be obtained by integrating fN over all the 6(N  1) variables that are not of interest. Together with Boltzmann’s assumption of molecular chaos, the equation of motion for a single particle distribution function then is given as follows: qf ½r; p; t  qf ½r; p; t  qf ½r; p; t  þ r_ þ p_ ¼ C½r; p; t coll ½27 qt qr qp

where C[r, p, t]coll is a collision term that derives from the fact that we averaged over the 6(N  1) dofs. Eqn [27] is called the Boltzmann transport equation (BTE). Solution techniques for eqn [27] that make use of cellular automata is briefly discussed in Section 1.18.2.2.1. 1.18.2.2.1

Approximations of the BTE Although it is a discretization method of the BTE, the lattice Boltzmann automaton (Chen and Doolen 1998) was introduced for the purpose of solving the NSE, which in turn can be derived from the BTE via a Chapman–Enskog Expansion (CEE). The NSE is the equation of motion of the first moment of the single-particle distribution function. An approximation of the BTE up to a finite order satisfies the NSE by design, even though it could be more general. At first, the most simplified lattice Boltzmann model still consistent with the NSE has been sought. Examples for simplifications that do not necessarily improve efficiency are the single relaxation time approximation (Bhatnagar et al. 1954) and models with a minimal set of dofs like the D3Q13 model (13 speeds in three dimensions) (d’Humieres et al. 2001). In these models, efficiency is often decreased by the introduction of instabilities. However, there is a very simple way to make even the worst model stable, e.g., through increased viscosity. By applying the scaling laws of fluid dynamics it is still possible to simulate any flow condition by choosing appropriate time steps and grid spacings for the model-dependent lowest attainable viscosity.

The performance penalty is severe and can be calculated from the same scaling laws. Two flow conditions are said to be similar if their Reynolds number is identical. The Reynolds number is defined as Re ¼ Lv/ . Here L is a characteristic flow length, v is the flow’s average velocity, and is the kinematic viscosity. Scaling laws are useful to derive numerically tractable dimensionless equations. It is easily seen that, for a constant Reynolds number, increasing viscosity can only be compensated for by increasing either v or L. Higher velocities will typically not result in more stable solutions since the lattice Boltzmann method is an explicit finite difference scheme and has to satisfy the CFL condition (Courant et al. 1967). The only remaining option is to take the expensive road and decrease L. This translates to finer grid spacings and shorter time steps. Since we are interested in time-varying 3D simulations we see that lowering the viscosity by one order of magnitude leads to four orders of magnitude decline in cell update operations and three orders of magnitude decline in the number of cells. In contrast, the above-mentioned simplifications result only in small gains per update step, typically smaller than a factor of 2 in memory and computational time. This cannot be compared at all with what we gain by finding a way of lowering the viscosity by many orders of magnitude without sacrificing stability. This is done, not by introducing new concepts or models, by merely removing other types of simplification, which crept into the scheme due to the effort to relate the lattice Boltzmann automata to the NSE rather than to see it as a discretization of the BTE. A recently developed new scheme called cascaded digital lattice Boltzmann automata (CDLBA) is able to overcome all the instability issues of conventional LBA (Geier et al. 2006a). 1.18.2.2.2

Applications of CDLBA The stability of CDLBA is impressively demonstrated by the turbulent flows in Figures 5 and 6.

Figure 5 Top and side view of the wake field of a rectangular object in a gas flow, Re 500 000.

Multiphysics and Multiscale Simulation

CDLB

y

x

speed magnitude 0.00 0.0375 0.0750 0.112 0.150

z

Imtek

Figure 6 Arbitrarily shaped object in a gas flow from an orifice, Re  1 000 000.

The flow over a rectangular plate has a Reynolds number of Re  500.000 while the orifice flow has Re  1.000.000. CDLBA is undoubtedly a major improvement for the calculation of flow fields in subsonic flow. In micrometer dimensions one may argue that turbulence does not play an important role. This is not true for all applications, especially when fast mixing devices are considered. In this application the onset of turbulence on short length scales is desired and the design of the respective structure can be supported by CDLBA, as it has already been done for the topology optimization of microflows (Geier et al. 2006b).

1.18.3 Conclusion Whenever there is talk of automatic simulation tools, meaning that we wish to eliminate the need for a human operator, we deal either with highly linear and therefore deterministic decisions and computations or with a highly restricted application area and need to invest enormous effort to make the implementations robust. Nevertheless, engineering design practice yearns for ever more automation as a means to become more efficient, especially in avoiding the tedious work of scanning for optimal design layouts for those cases where at least the target performance quantities are evaluable. In our search for better tools to achieve these goals, we and others are working at cointegrating continuum and mesoscopic methods, and on extending algorithms to become more efficient and robust, and general enough to cope with the insatiable demand of MEMS engineers for design support tools in an ever-broadening application field.

555

Our choice for the immediate future has fallen on the level set method, treated as a coupled scalar field and tied to the rest of the physics, as a way to capture the engineering design topology. It provides a natural setting for a range of complex phenomena when performing topology optimization together with the finite element and related methods. For soft matter physics, two particle approaches are followed: the DPD method and the digital lattice Boltzmann method. We are exploring ways to tie them together with conventional (continuum) methods, so as to achieve efficient multiscaling simulators, and we are embedding them in a topology optimization cycle, so as to bring these methods also into the design tool domain.

Acknowledgments The authors wish to thank M. Geier, D. Kauzlaric, and L. Pastewka for their support in preparing the figures.

References Allaire G, Jouv F, Toader A 2004 Structure optimization using sensitivity analysis and a level-set method. J. Comput. Phys. 194, 363–93 Allgower E, Georg K 2003 Introduction to Numerical Continuation Methods. SIAM, Philadelphia, PA Arora J S, Wang Q 2005 Review of formulations for structural and mechanical system optimization. Struct. Multidisc. Optim. 30, 251–72 Ascher U M, Petzold L R 1998 Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equation. SIAM, Philadelphia, PA Beckett G, Mackenzie J A, Robertson M L 2001 A moving mesh finite element method for the solution of two-dimensional Stefan problems. J. Comput. Phys. 182, 478–95 Belytschko T, Liu W K, Moran B 2000 Nonlinear Finite Elements for Continua and Structures. Wiley, New York Bendsoe M P, Kikuchi N 1988 Generating optimal topologies in optimal design using a homogenization method. Comp. Meth. Appl. Mech. Eng. 71, 197–224 Bendsoe M P, Sigmund O 2003 Topology Optimization Theory, Methods and Applications. Springer, Berlin Bhatnagar P, Gross E, Krook M 1954 A model for collisional processes in gases I: Small amplitude processes in charged and neutral one-component system. Phys. Rev. 94, 511–25 Boek E S, Coveney P V, Lekkerkerker H N W, van der Schoot P 1997 Simulating the rheology of dense colloidal suspensions using dissipative particle dynamics. Phys. Rev. E 55, 3124–33 Boneyt Avalos J, Mackie A D 1997 Dissipative particle dynamics with energy conservation. Europhys. Lett. 40, 141–6 Cao W, Huang W, Russell R D 2003 Approaches for generating moving adaptive meshes: Location versus velocity. Appl. Numer. Math. 47, 121–38 Chen S, Doolen G 1998 Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech. 30, 329–64

556 Multiphysics and Multiscale Simulation Costa J, Alves M K 2003 Layout optimization with H-adaptivity of structures. Int. J. Numer. Meth. Eng. 58, 83–102 Courant R, Friedrichs K, Lewy H 1967 On the partial difference equations of mathematical physics. IBM J. Res. Dev. 11(2), 215–324 Coveney P V, Novik K E 1996 Computer simulations of domain growth and phase separation in two-dimensional immiscible fluids using dissipative particle dynamics. Phys. Rev. E 54, 5134–41 d’Humieres D, Bouzidi M, Lallemand P 2001 Thirteen-velocity three-dimensional lattice Boltzmann model. Phys. Rev. E 63, 066702 Du J, Olhoff N 2004 Topological optimization of continuum structures with design-dependent surface loading – Part 1: New computational approach for 2D problems. Struct. Multidisc. Optim. 27, 151–65 Dwinzel W, Yuen D A 2000 Matching macroscopic properties of binary fluids to the interactions of dissipative particle dynamics. Int. J. Mod. Phys. C 11, 1–25 Dwinzel W, Yuen D A, Boryczko K 2002 Mesoscopic modeling of colloids simulated with dissipative particle dynamics. J. Model. 8, 33–43 Espanol P, Revenga M 2003 Smoothed dissipative particle dynamics. Phys. Rev. E 67, 026705 Felippa C A 1999 Partitioned analysis of coupled mechanical systems. Report Cu-CAS-99-06 Flekkoy E G, Coveney P V 1999 From molecular dynamics to dissipative particle dynamics. Phys. Rev. Lett. 83, 1775–8 Flekkoy E G, Wagner G, Feder J 2000 Hybrid the large model for combined particle and continuum dynamics. Europhys. Lett. 52, 271–6 Geier M, Greiner A, Korvink J G 2006a Cascaded digital lattice Boltzmann automata for high Reynolds number flow. Phys. Rev. E 73(6), 066705 Geier M, Liu Z, Greiner A, Korvink J G 2006b Lattice Boltzmann based structural topology optimization of micro channels. Second Int. Conf. Transport Phenomena in Micro and Nano devices, Barga, Italy Gockenbach M S 2002 Partial Differential Equations: Analytical and Numerical Methods. SIAM, Philadelphia, PA Groot R D, Warren P B 1997 Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J. Chem. Phys. 107, 4423–35 Hammer V B, Olhoff N 2000 Topology optimization of continuum structures subjected to pressure loading. Struct. Multidisc. Optim. 19, 85–92 Hoogerbrugge P J, Koelman J M V A 1992 Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys. Lett. 19, 155–60 http://www.ansys.com/ http://www.femlab.com Jakobsen L A, Lund E, Møller H 2001 Shape sensitivity analysis of time dependent fluid-structure interaction problems using the ALE method. Proc. 4th World Congr. Structural and Multidisciplinary Optimization, Dalian, China Kauzlaric D, Greiner A, Korvink J G 2003a An alternative simulation method for micro injection moulding. Conference, Proceedings Polytronic 2003, Montreux, Switzerland, pp. 79–85 Kauzlaric D, Greiner A, Korvink J G 2003b Modelling of microrheological behaviour of fluids with dissipative particle dynamics. VDI Berichte 1803, pp. 179–82 Kauzlaric D, Greiner A, Korivink J G 2004 Modelling microrheological effects in micro injection moulding with dissipative particle dynamics. Tech. Proc. 2004 Nanotechnol. Conf. Trade Show, Boston, MA, USA, vol. 2, pp. 454–7

Kauzlaric D, Greiner A, Korvink J G, Schulz M, Heldele R 2005 Modelling micro-PIM. In: Balites H, Brand O, Fedder G K, Hierold C, Korvink J G, and Tabata O (eds.) Microengineering of Metals and Ceramics. Wiley-VCHVerlag, Weinheim, Germany, ISBN 3–527–31208–0 Koumoutsakos P 2005 Multiscale flow simulations using particles. Annu. Rev. Fluid Mech. 37, 457–87 Liu Z, Korvink J G, Huang R 2005 Structure topology optimization: Fully coupled level set method via FEMLAB. Struct. Multidisc. Optim. 29, 407–17 Morton K W, Mayers D F 1994 Numerical Solution of Partial Differential Equations: An Introduction. Cambridge University Press, Cambridge Nathan A, Baltes H 1999 Microtransducer CAD, Physical and Computational Aspects. Springer, Berlin Nielaba P, Mareschal M, and Ciccotti G (eds.) (2002) Bridging Time Scales: Molecular Simulations for the Next Decade. Lecture Notes in Physics. Springer, Berlin, Vol. 605 Nocedal J, Wright S J 1999 Numerical Optimization. Springer, Berlin Osher S, Fedkiw R 2002 Level Set Methods and Dynamic Implicit Surfaces. Springer, Berlin Paganobarraga I, Frenkel D 2001 Dissipative particle dynamics for interacting systems. J. Chem. Phys. 115, 5015–26 Pastewka L, Kauzlaric D, Greiner A, Korvink J G 2006 Thermostat with local heat bath coupling for exact energy conservation in dissipative particle dynamics. Phys. Rev. E 73(3), 037701 Pool R, Bolhuis P G 2006 Can purely repulsive soft potentials predict micelle formation correctly. Phys. Chem. Chem. Phys. 8, 941–8 Romanowicz B F 1998 Methodology for the Modeling and Simulation of Microsystems. Kluwer Academic Publishers, Dordrecht, The Netherlands Rozvany G I N 2001 Aims, scope, methods, history and unified terminology of computer-aided topology optimization in structural mechanics. Struct. Multidisc. Optim. 21, 90–108 Sapiro G 2001 Geometric Partial Differential Equations and Image Analysis. Cambridge University Press, Cambridge Schulz S G, Kuhn H, Schmid G, Mund C, Venzmer J 2004 Phase behaviour of amphiphilic polymers: A dissipative particle dynamics study. Colloid. Polym. Sci. 283, 284–90 Senturia S D 2000 Microsystem Design. Kluwer Academic Publishers, Dordrecht, The Netherlands Sigmund O 2001a Design of multiphysics actuators using topology optimization – Part I: One-material structures. Comput. Meth. Appl. Mech. Eng. 190, 6577–604 Sigmund O 2001b Design of multiphysics actuators using topology optimization – Part II: Two-material structures. Comput. Meth. Appl. Mech. Eng. 190, 6605–27 Thompson J F, Warsi Z U A, Mastin C W 1985 Numerical Grid Generation: Foundations and Applications. North-Holland, New York Wolf-Gladrow D A 2000 Lattice Gas Cellular Automata and Lattice Boltzmann Models: An Introduction. Springer, Berlin Yin L, Ananthasuresh G K 2002 A novel topology design scheme for the multi-physics problems of electro-thermally actuated compliant micromechanisms. Sens. Actuators A 97, 599–609 Yoon G H, Kim Y Y 2006 The element connectivity parameterization formulation for the topology design optimization of multiphysics systems. Int. J. Numer. Meth. Eng. 64, 1649–77 Zienkiewicz O C, Taylor R, Zhu J 2005 The Finite Element Method. Elsevier, Amsterdam

Multiphysics and Multiscale Simulation

Biographies Dr. Jan G. Korvink obtained his M.Sc. in computational mechanics from the University of Cape Town in 1987, and his Ph.D. in applied computer science from the ETH Zurich in 1993. After his graduate studies, he joined the Physical Electronics Laboratory of the ETH Zurich, where he established and led the MEMS Modelling Group. This was followed by a move in 1997 to the Albert Ludwig University in Freiburg, Germany, where he holds a Chair position in microsystem technology and runs the Laboratory for Microsystem Simulation. He has written more that 180 journal and conference papers in the area of microsystem technology, and co-edits the review journal Applied Micro and Nanosystems, see http:// www.wiley-vch.de/books/info/amn. His research interests cover the modelling and simulation and low cost fabrication of microsystems.

557

Andreas Greiner is a senior scientist and head of the group of physical modeling at the laboratory for Simulation of the Department of Microsystems Engineering (IMTEK), University of Freiburg, Germany. He graduated in physics from the University of Stuttgart, Germany, where he received his PhD in Theoretical Physics in 1992 and where he was concerned with the simulation of carrier transport in semiconductor lasers until 1995. From June 1995 on he held several grants as a research fellow at the University of Bologna, the University of Lecce and the CNRS research center for Microoptoelectronics in Montpellier. His research interest is the multiscale modeling and simulation of micro- and nanosystems as well as process simulation.

Zhenyu Liu obtained his PhD in mechanical engineering from the Dalian University of Technology, China in 2000. From 2001 to 2002, he worked as postdoc in the Department of Microsystems Engineering (IMTEK), University of Freiburg, Germany. From 2003 he joined IMTEK as a senior scientist and group leader at the laboratory for Simulation. His research interests cover the simulation and optimization of devices in MEMS.