Suche

Optimierungsprozesse mit partiellen Differentialgleichungen

'Multilevel Based All-At-Once Methods in PDE Constrained Optimization with Applications to Shape Optimization of Active Microfluidic Biochips

Projektstart: 01.02.2006
Projektträger: DFG (Deutsche Forschungsgemeinschaft)
Projektverantwortung vor Ort: Prof. Dr. R. Hoppe, Prof. Dr. K. Siebert, Prof. Dr. A. Wixforth
Beteiligte Wissenschaftler der Universität Augsburg: Prof. Dr. K. Siebert
Prof. Dr. A. Wixforth

Zusammenfassung

This project within the area of PDE constrained optimization focuses on the development, analysis and implementation of optimization algorithms that combine efficient solution techniques from the numerics of PDEs, namely multilevel iterative solvers, and state-of-the-art optimization approaches, the so-called `all-at-once' optimization methods. It is well-known that multilevel techniques provide efficient PDE solvers of optimal algorithmic complexity. On the other hand, optimization methods within the all-at-once approach, such as sequential quadratic programming (SQP) methods and primal-dual Newton interior-point methods, have the appealing feature that in contrast to more traditional approaches, the numerical solution of the state equations is an integral part of the optimization routine. This is realized by incorporating the PDEs as constraints into the optimization routine. These strategies allow to save a considerable amount of computational work compared to methods that treat the PDE solution as an implicit function of the control/design variables. Moreover, the proper combination of multilevel techniques and optimization algorithms makes it possible to extract essential structural information from the originally infinite dimensional optimization problem. This can not be done with respect to a single grid. We aim to develop and analyze multilevel preconditioners for optimization subproblems arising in SQP and primal-dual Newton interior-point methods including strategies to control the level of inexactness allowable in optimization subproblems, when using iterative subproblem solvers. Moreover, we will investigate strategies to use multilevel methods for detection of negative curvature and in path following methods.

Beschreibung

Multilevel Based All-At-Once Methods in PDE Constrained Optimization with Applications to Shape Optimization of Active Microfluidic Biochips


Project leaders:

a) Prof.Dr. Ronald H.W. Hoppe
Universitätsprofessor (C4)
geb. 10.04.1951, Nationalität: deutsch
Adresse: Universität Augsburg
Lehrstuhl für Angewandte Analysis/Numerik
Institut für Mathematik
Universitätsstr. 14, D-86159 Augsburg
Tel.: 0821-598-2194, Fax: 0821-598-2339
E-mail: hoppe@math.uni-augsburg.de

Address: University of Houston
Department of Mathematics
651 Ph. G. Hoffman
Houston, TX 77204-3008, U.S.A.
Tel.: +1-713-743-3452, Fax: +1-713-743-3505
E-mail: rohop@math.uh.edu

b) Prof.Dr. Kunibert Siebert
Universitätsprofessor (C3)
geb. 09.05.1963, Nationalität: deutsch
Adresse: Universität Augsburg
Lehrstuhl für Angewandte Analysis/Numerik
Institut für Mathematik
Universitätsstr. 14, D-86159 Augsburg
Tel.: 0821-598-2190, Fax: 0821-598-2339
E-mail: siebert@math.uni-augsburg.de

c) Prof.Dr. Achim Wixforth
Universitätsprofessor (C4)
geb. 26.05.1956, Nationalität: deutsch
Adresse: Universität Augsburg
Lehrstuhl für Experimentalphysik I
Institut für Physik
Universitätsstr. 1, D-86159 Augsburg
Tel.: 0821-598-3300, Fax: 0821-598-3225
E-mail: achim.wixforth@physik.uni-augsburg.de

Summary

This project within the area of PDE constrained optimization focuses on the development, analysis and implementation of optimization algorithms that combine efficient solution techniques from the numerics of PDEs, namely multilevel iterative solvers, and state-of-the-art optimization approaches, the so-called `all-at-once' optimization methods. It is well-known that multilevel techniques provide efficient PDE solvers of optimal algorithmic complexity. On the other hand, optimization methods within the all-at-once approach, such as sequential quadratic programming (SQP) methods and primal-dual Newton interior-point methods, have the appealing feature that in contrast to more traditional approaches, the numerical solution of the state equations is an integral part of the optimization routine. This is realized by incorporating the PDEs as constraints into the optimization routine. These strategies allow to save a considerable amount of computational work compared to methods that treat the PDE solution as an implicit function of the control/design variables.
Moreover, the proper combination of multilevel techniques and optimization algorithms makes it possible to extract essential structural information from the originally infinite dimensional optimization problem. This can not be done with respect to a single grid. We aim to develop and analyze multilevel preconditioners for optimization subproblems arising in SQP and primal-dual Newton interior-point methods including strategies to control the level of inexactness allowable in optimization subproblems, when using iterative subproblem solvers. Moreover, we will investigate strategies to use multilevel methods for detection of negative curvature and in path following methods.

The developed PDE constrained optimization algorithms will be applied to the optimal design of microfluidic biochips with emphasis on the optimization of the geometry of the devices (shape optimization) in order to achieve an optimal operational behavior. Here, the state equations consist of a system of partial differential equations describing the transport of microfluids driven by piezoelectrically agitated surface acoustic waves.
During the last couple of years, such biochips have attracted a considerable amount of interest, since pharmacology, molecular biology, and clinical diagnostics require the precise handling of precious, tiny samples and costly reagents in amounts of nanoliters. Biochips can transport such volumes and perform biochemical analysis of the samples. Microfluidic biochips and microarrays are used in pharmaceutical, medical and forensic applications as well as in academic research and development for high throughput screening, genotyping and sequencing by hybridization in genomics, protein profiling in proteomics, and cytometry in cell analysis. Traditional technologies rely on fluorescent dyes, radioactive markers, or nanoscale gold-beads based on positive hybridization processes. However, these methods only allow a relatively small number of DNA probes per assay, and they only yield endpoint results and do not provide information about the kinetics of the processes. With the need for better sensitivity, flexibility, cost-effectiveness and a significant speed-up of hybridization, the current technological trend is obtained by the integration of the microfluidics on the chips itself.
A new type of nanotechnological devices are surface acoustic wave driven microfluidic biochips. The experimental technique is based on piezoelectrically actuated surface acoustic waves on the surface of a chip which transport the droplet containing probe along a lithographically produced network channels to marker molecules placed at prespecified surface locations. Hence, these microfluidic biochips allow the in-situ investigation of the dynamics of hybridization processes with extremely high time resolution.

The scientific merits of the proposal are the design, analysis and implementation of efficient algorithmic tools for a class of challenging problems in nonlinear optimization and the demonstration of their performance by the application to a real-life problem that has a significant impact on material sciences and life sciences.
The broader technological impact of the project is a better design of microfluidic biochips with an improved performance for biochemical applications.

PDE Constrained Optimization

Optimization problem with constraints given by partial differential equations (PDE) can be written as follows

\begin{subequations} \begin{align} &\text{minimize} \quad J(y,u) \\ &\text... ... &\text{subject to:} \quad e(y,u) = 0 . \end{align}%EQNO: \end{subequations}
Here, $ J(\cdot,\cdot) : Y \times U \rightarrow \textbf{\rm l\hskip-0.5mm R}$ denotes the objective functional which depends on the state variables $ y $ in the state space $ Y $ and the design variables $ u $ in the design space $ U $. The equation (1b) corresponds to the PDE, whereas $ K_{1} \subset Y , K_{2} \subset U $ refer to the sets of admissible states and design variables, respectively.
The numerical solution of PDE constrained optimization problems is currently an utmost active area of research. The most challenging theoretical and practical issues can be summarized as follows:
  • The structural interaction between the optimization and the underlying PDE and the discretization processes have a tremendous impact on the solution of the problem which is reflected by the paradigms ''discretize first, then optimize'', and conversely, ''optimize first, then discretize''. For instance, following the principle ''discretize first, then optimize'', the granularity of the discretizations is a critical issue, since due to the discretization, we are faced with a certain loss of the inherent structure of the infinite dimensional problem. On the other hand, the principle ''optimize first, then discretize'' preserves the structure of the problem, but the optimality conditions can only be realized approximately. The ultimate goal is to find some sort of best compromise between these two principles.
  • In contrast to more traditional optimization methodologies that rely on a separate treatment of the optimization issue and the solution of the state equation, recent interest is on so-called ''all-at-once'' approaches. Here, the numerical solution of the state equation is an integral part of the optimization routine which may result in a considerable savings of computational time (cf., e.g., [29,30,86,145,152,176,200]).
  • The numerical approaches typically lead to large-scale nonlinear programming problems. With regard to algorithmic complexity, their numerical solution requires the use of efficient iterative schemes such as multilevel techniques.
Simplified problems in structural optimization have already been addressed by Bernoulli, Euler, Lagrange and Saint-Venant. However, it became its own discipline during the second half of the last century when the rapidly growing performance of computing platforms and the simultaneously achieved significant improvement of algorithmic tools enabled the appropriate treatment of complex problems (cf. [5,19,20,50,61,104,178,187,194,203] and the references therein). The design criteria in structural optimization are determined by a goal oriented operational behavior of the devices and systems under consideration and typically occur as a nonlinear, often non convex, objective functionals which depend on the state variables describing the operational mode and the design variables determining the shape and the topology. The state variables have to satisfy differential equations or systems thereof representing the underlying physical laws. Technological aspects are taken into account by constraints on the state and design variables which may occur both as equality and inequality constraints in the model.

The discretization of PDE constrained optimization problems typically gives rise to large-scale, equality and inequality constrained nonlinear programming problems of the form

\begin{subequations} \begin{align} &\text{minimize} \quad f(x_{1},x_{2}) \ ... ...ject to:}} \quad g(x_{1},x_{2}) \ge 0 , \end{align}%EQNO:1 \end{subequations}
where $ f : R^{n_{1}} \times \textbf{\rm l\hskip-0.5mm R}^{n_{2}} \rightarrow \textbf{\rm l\hskip-0.5mm R}$ is the discretized objective functional depending on the discrete state variables $ x_{1} \in \textbf{\rm l\hskip-0.5mm R}^{n_{1}} $ and design variables $ x_{2} \in \textbf{\rm l\hskip-0.5mm R}^{n_{2}} $, (2b) with $ h = (h_{1},h_{2})^{T} : R^{n_{1}} \times \textbf{\rm l\hskip-0.5mm R}^{n_{2}} \rightarrow \textbf{\rm l\hskip-0.5mm R}^{m} , m \ge n_{1}, $ comprises the discretized state equation $ h_{1} : R^{n_{1}} \times \textbf{\rm l\hskip-0.5mm R}^{n_{2}} \rightarrow \textbf{\rm l\hskip-0.5mm R}^{n_{1}} $ and possibly further equality constraints $ h_{2} : R^{n_{1}} \times \textbf{\rm l\hskip-0.5mm R}^{n_{2}} \rightarrow \textbf{\rm l\hskip-0.5mm R}^{m-n_{1}} $, whereas (2c) with $ g: R^{n_{1}} \times \textbf{\rm l\hskip-0.5mm R}^{n_{2}} \rightarrow \textbf{\rm l\hskip-0.5mm R}^{p} $ represents inequality constraints on the state and design variables.
If we define the feasible set $ \mathcal{F} $ according to

$\displaystyle \mathcal{F} := \{ (x_{1},x_{2}) \vert h(x_{1},x_{2}) = 0 , \ g(x_{1},x_{2}) \ge 0 \} , $

then (2a)-(2c) can be written more concisely as

$\displaystyle \min\limits_{(x_{1},x_{2}) \in \mathcal{F}} f(x_{1},x_{2}) .$ (3)

Coupling the constraints by Lagrangian multipliers $ \lambda \in \textbf{\rm l\hskip-0.5mm R}^{m} $ and $ \mu \in \textbf{\rm l\hskip-0.5mm R}^{p}_{+} $, we are led to the saddle point problem

$\displaystyle \min\limits_{(x_{1},x_{2}) \in \textbf{\rm l\hskip-0.5mm R}^{n_{1... ...^{m} \times \textbf{\rm l\hskip-0.5mm R}_{+}^{p}}L(x_{1},x_{2},\lambda,\mu) ,$ (4)

where $ L $ stands for the Lagrangian
$\displaystyle L(x_{1},x_{2},\lambda,\mu) := f(x_{1},x_{2}) + \lambda^{T} h(x_{1},x_{2}) - \ \mu^{T} g(x_{1},x_{2}) .$     (5)

The Karush-Kuhn-Tucker (KKT) conditions give rise to the nonlinear system

$\displaystyle \left[ \begin{array}{c} \nabla f(x) + \nabla h(x) \lambda +... ... \mu h(x) g(x) - z D_{z} D_{\mu} e \end{array} \right] = 0 .$ (6)

Here, $ z \in \textbf{\rm l\hskip-0.5mm R}^{p}_{+} $ is the vector of slack variables, $ D_{z} =$   diag$ (z_{1},...,z_{p}) \ , D_{\mu} =$   diag$ (\mu_{1},\ldots, \mu_{p}) $ and $ e = (1,...,1)^{T} $. If Newton's method is applied to the KKT conditions, each Newton step requires the solution of a linear algebraic system representing the optimality conditions of a related quadratic programming (QP) problem. Hence, Newton methods can be interpreted in the framework of sequential quadratic programming (SQP) which is the most successful method for solving constrained nonlinear optimization problems [35,183].
As far as the appropriate treatment of the inequality constraints is concerned, a local optimum can be approximated from within the feasible set, which is the idea behind interior-point methods. The so-called interior-point revolution in continuous optimization started in the eighties of the last century with Karmarkar's polynomial-time linear programming algorithm [163]. It was immediately found [92] that there is a close relationship to barrier functions which had been used long time before for inequality constrained nonlinear programming problems. Nowadays, interior-point methods are well established tools for constrained nonlinear optimization problems (cf., e.g., [15,44,77,78,79,85,89,92,211,212,218,228,229,230]).
An alternative is to use active set strategies where the iterates are not necessarily restricted to feasibility (cf., e.g., [22,117] for recent work in the context of constrained optimal control problems).
Reliable error estimates and mesh adaptation procedures for numerically solving PDEs have reached some state of maturity as documented by a series of monographs on this issue published during the last decade [4,10,13,74,181,213]. However, error estimates, local error indicators and adaptivity procedures for PDEs cannot be carried over directly to optimization problems with PDE constraints. In case of unconstrained optimal control problems, pioneering work has been done on the basis of the goal oriented weighted dual approach (cf., e.g., [13,18] and the references therein). However, if some of the system variables, like the control or the state, have to satisfy certain inequality constraints, which may be technologically or physically motivated, the literature on mesh adaptivity for optimal control problems for PDEs is extremely scarce [119,172,173,174]. Also, concepts known for obstacle problems cannot be adopted to the inequality constrained optimization context [3,15,133,161,165].
Mesh adaptivity for structural optimization problems has been considered in the engineering literature (cf., e.g., [120]), but is restricted to the application of standard error indicators for the error in the state variable and, in particular, does not take into account the impact of inequality constraints and their approximation on the quality of the mesh and the accuracy of the approximate solution.

Microfluidic Biochips

Biochips, of the microarray type, are fast becoming the default tool for combinatorial chemical and biological analysis in environmental and medical studies. The goal of this project is to develop, analyze and implement state-of-the-art optimization techniques for an improved design of these devices.
Programmable biochips are miniaturized biochemical labs that are physically and/or electronically controllable. The technology combines digital photolithography, microfluidics and chemistry. The precise positioning of the samples (e.g., DNA solutes or proteins) on the surface of the chip in picoliter to nano liter volumes can be done either by means of external forces (active devices) or by specific geometric patterns (passive devices).
During the last couple of years, such biochips have attracted a considerable amount of interest, since pharmacology, molecular biology, and clinical diagnostics require the precise handling of precious, tiny samples and costly reagents in amounts of nanoliters. Biochips can transport such volumes and perform biochemical analysis of the samples. Nanoliter fluidic biochips and microarrays are used in pharmaceutical, medical and forensic applications as well as in academic research and development for high throughput screening, genotyping and sequencing by hybridization in genomics, protein profiling in proteomics, and cytometry in cell analysis [111]. Traditional technologies rely on fluorescent dyes, radioactive markers, or nanoscale gold-beads based on positive hybridization processes. However, these methods only allow a relatively small number of DNA probes per assay, and they only yield endpoint results and do not provide information about the kinetics of the processes. With the need for better sensitivity, flexibility, cost-effectiveness and a significant speed-up of hybridization, the current technological trend is obtained by the integration of the microfluidics on the chips itself.
A substantial amount of research activity in this direction is in the area of microchip capillary electrophoresis (cf., e.g., [3,81,110]). Here, gel-filled capillary channels are placed on micromachined glass or polymer microreplicated chips. More recent and novel nanotechnological devices are surface acoustic wave driven nanoliter fluidic biochips and pressure-driven nano-chambered, multi-channel microfluidic biochips. The former technique is based on piezoelectrically actuated surface acoustic waves on the surface of a chip which transport the droplet containing probe along a lithographically produced network to marker molecules placed at prespecified surface locations (cf., e.g., [13,114,139,140]). The latter one adequately combines digital photolithography, microfluidics and chemistry (cf., e.g., [38,39,40,89,90,110,146]). These nanoliter fluidic biochips allow the in-situ investigation of the dynamics of hybridization processes with extremely high time resolution.
By changing the surface chemistry appropriately, a fluidic network is produced on the chip: Without mechanical tools the chip is equipped with paths on which samples (and reagents) propagate as if on tracks. This is done lithographically by a lateral modulation of the wetting properties of the free surface which leads to pronounced hydrophilic and super-hydrophobic regions with significantly different wetting angles (cf., e.g., [27,51,73,95,99]). Small amounts of reagents are confined to these tracks in contrast to mechanical barriers used in conventional microfluidics.
The active devices which will be considered in this project are nanoliter fluidic biochips where the core of the technology are nanopumps featuring surface acoustic waves generated by electric pulses of high frequency. These waves propagate like a miniaturized earthquake (nanoscale earthquake) and in this way transport liquids along the surface of the chip.
Figure 1 below gives an illustration of a nano titration chip. On the fluidic network a small portion of titrate solution (middle) is separated from a larger volume (right). Surface acoustic waves transport this quantity towards the analyte (left) at the reaction site. Once a critical concentration is attained, it can be either detected by a change of the color of the analyte or a change of the conductivity. In the latter case, this can be easily measured by a sensor that is integrated on the same chip.

Figure 1: Surface acoustic wave micropump
\includegraphics[width=12.5cm,height=7.5cm]{network2.eps}

Surface Acoustic Waves (SAW) have been used for a long time in high frequency applications (cf., e.g., [21,101,105], and [141]). Using SAW-principles it is now possible to combine microelectronics and biochemistry. Modern semiconductor technology enables the cost-effective production of devices that unify biological functionality, sensors and pumps for the transport of samples. These devices can be easily integrated in electronic systems like those that are used in point-of-care diagnostics (see [4,116,117,126]).
The micropump consists of a piezoelectric substrate equipped with so-called interdigital transducers on the surface. Radio-frequency signals are fed into those transducers and are converted to a deformation of the crystal underground. In this way, a mechanical wave is launched across the surface with wavelengths in the range of a few microns and amplitudes about only a nanometer. Liquids on the surface are subject to the vibrating force and absorb parts of the energy. The absorbtion of energy for various frequencies depends on the density and viscosity.
Surface acoustic waves of larger amplitudes move liquid droplets as a whole whereas low power surface acoustic waves induce some sort of an internal streaming. The latter case enables to construct surface acoustic wave based nanomixers. If the frequency of the surface acoustic wave is changed, different streaming patterns are induced and superimposed within the droplet that leads to a homogeneous blend of the water and the probe much faster than by more conventional diffusion type microfluidic mixing techniques.
Figure 2 illustrates the effect of nanomixing in case of the dissolution of a fluorescent dye deposited on the chip surface with agitation (acoustically induced mixing) and without agitation.

Figure 2: Surface acoustic wave nanomixer (with agitation (top), without agitation (below)
\includegraphics[width=12.0cm,height=3.0cm]{mischen_3er.ps}

By using surface acoustic wave nanopumps, different reagents can be efficiently mixed, separated or moved to different reaction sites on the chip. Compared with conventional micro titer plates, the respective volumes are reduced by several orders of magnitudes.

Previous work of the project leaders

Primal-dual Newton interior point methods

Discretization of the objective functional, the state equation and the equality constraints by finite elements with respect to a simplicial triangulation $ \mathcal{T}_{h} $ of the computational domain $ \Omega $ leads to the finite dimensional constrained minimization problem: Find $ ({\bf u}_{\bf h},$$ \mbox {\boldmath $\alpha $}) \in \textbf{\rm l\hskip-0.5mm R}^{n_{h}} \times \textbf{\rm l\hskip-0.5mm R}^{m} $ such that

$\displaystyle J_{h}({\bf u}_{\bf h},$$\displaystyle \mbox {\boldmath$\alpha $}) = \inf\limits_{{\bf v}_{\bf h},\mbox {\boldmath$\mu $}} J_{h}({\bf v}_{\bf h},\mbox {\boldmath$\mu $})$     (7)

subject to the constraints
$\displaystyle A_{h}($$\displaystyle \mbox {\boldmath$\alpha $}) {\bf u}_{\bf h}$ $\displaystyle =$ $\displaystyle {\bf b}_{h} ,$ (8)
$\displaystyle c_{h}($$\displaystyle \mbox {\boldmath$\alpha $})$ $\displaystyle :=$ $\displaystyle C ,$ (9)
$\displaystyle \alpha_{i}^{(min)} \le$ $\displaystyle \alpha_{i}$ $\displaystyle \le \alpha_{i}^{(max)} \quad , \quad 1 \le i \le m ,$ (10)

where (8) denotes the FE discretized state equation with the stiffness matrix $ A_{h}($$ \mbox {\boldmath $\alpha $}) \in \textbf{\rm l\hskip-0.5mm R}^{n_{h} \times n_{h}} $ and the load vector $ {\bf b}\in \textbf{\rm l\hskip-0.5mm R}^{n_{h}} $.
For the efficient numerical solution of (8)-(10) we have used an ''all-at-one approach'' in terms of a primal-dual Newton interior-point method (cf., e.g., [123,124,140,141,145,146,149,150,152]). The interior-point aspect is to take care of the inequality constraints (10) by parameterized logarithmic barrier functions
$\displaystyle B_{h}^{(p)}({\bf u}_{\bf h},$$\displaystyle \mbox {\boldmath$\alpha $}) := J_{h}({\bf u}_{\bf h},\mbox {\... ...lpha_{i}^{(min)} ) + \mbox{log} ( \alpha_{i}^{(max)} - \alpha_{i} ) \big] ,$      

whereas the primal-dual aspect is to couple the remaining equality constraints (8),(9) by Lagrangian multipliers $ \mbox {\boldmath $\lambda $}_{\bf h} \in \textbf{\rm l\hskip-0.5mm R}^{n_{h}} $ and $ \eta_{h} \in \textbf{\rm l\hskip-0.5mm R}$ which gives rise to the saddle point problem
$\displaystyle \inf\limits_{{\bf u}_{\bf h},\mbox {\boldmath$\alpha $}} \sup\l... ...h},\mbox {\boldmath$\alpha $},\mbox {\boldmath$\lambda $}_{\bf h},\eta_{h}) ,$     (11)

where the Lagrangian is given by

$\displaystyle L_{h}^{(p)}({\bf u}_{\bf h},$$\displaystyle \mbox {\boldmath$\alpha $},\mbox {\boldmath$\lambda $}_{\bf h},\e... ...{\bf b}_{\bf h}) + \eta_{h} ( c_{h}(\mbox {\boldmath$\alpha $}) - C) . $

Finally, the Newton aspect is to apply Newton's method to the Karush-Kuhn-Tucker conditions associated with (11)
$\displaystyle \nabla_{{\bf u}_{\bf h}} L_{h}^{(p)}$ $\displaystyle =$ $\displaystyle \nabla_{{\bf u}_{\bf h}} J_{h} + A_{h}($$\displaystyle \mbox {\boldmath$\alpha $})^{T} \mbox {\boldmath$\lambda $}_{\bf h} = 0 ,$  
$\displaystyle \nabla_{\mbox {\boldmath$\alpha $}} L_{h}^{(p)}$ $\displaystyle =$ $\displaystyle \nabla_{\mbox {\boldmath$\alpha $}} J_{h} + \ \partial_{\mbox ... ...x {\boldmath$\alpha $}) - p D_{1}^{-1} {\bf e}+ p D_{2}^{-1} {\bf e} = 0 ,$  
$\displaystyle \nabla_{\mbox {\boldmath$\lambda $}_{\bf h}} L_{h}^{(p)}$ $\displaystyle =$ $\displaystyle A_{h}($$\displaystyle \mbox {\boldmath$\alpha $}) {\bf u}_{\bf h} - {\bf b}_{\bf h} = 0 ,$  
$\displaystyle \nabla_{\eta_{h}} L_{h}^{(p)}$ $\displaystyle =$ $\displaystyle c_{h}($$\displaystyle \mbox {\boldmath$\alpha $}) - C = 0 ,$  

where $ D_{1} :=$   diag$ (\alpha_{i} - \alpha_{i}^{(min)}) , D_{2} :=$   diag$ (\alpha_{i}^{(max)} - \alpha_{i}) $ and $ {\bf e}:= (1,...,1)^{T} $.
The resulting primal-dual Hessian system is solved by a null-space approach using right-transforming iterations (originally suggested in [224]) with respect to the special block structure of the primal-dual Hessian. After computation of the Newton increments, a line search is performed and combined with a watchdog strategy for convergence monitoring to ensure convergence to a local minimum. We refer to [145] and [149] for details.

Shape optimization of electrorheological shock absorbers

Electrorheological fluids are concentrated suspensions of small electrically polarizable particles with diameters in the range of micrometers dissolved in nonconducting silicon oils. The rheological effect is based on the fact that under the influence of an outer electric field the particles form chains along the field lines and then aggregate to form larger and larger columns. The impact on the macroscopic scale consists in a rapid change of the rheological properties which happens within a few milliseconds. The viscosity increases in the direction orthogonal to the electric field such that the character of the fluid changes from liquid to almost solid. Under the action of large stresses, depending on the electric field strength (field dependent yield stress), the columnar structures break such that the viscosity decreases and the fluid behaves less anisotropic. The process is reversible, i.e., the viscosity decreases with decreasing field strength and the fluid behaves again like a Newtonian fluid for vanishing outer electric field.
Therefore, electrorheological fluids are used in all technological processes where a controlled power transmission plays a significant role. The field of applications ranges from automotive shock absorbers and actuators in hydraulic systems to tactile devices for virtual reality [76].

Figure: Electrorheological shock absorber (left) and its characteristics (right)
\begin{figure}\vbox{\hbox to\textwidth{\hspace{0.8cm} \psfig{figure=ERF_SD.ps,wi... ...5.5cm,height=5.5cm}\hfill}} \vspace{0.2cm} \hspace{2mm}\vspace{2mm} \end{figure}

In particular, we have been concerned with the optimization of the shape of the walls of an ERF shock absorber (cf. Figure 3 (left)) in a vicinity of the inlet and outlet boundary of the ERF transfer ducts (cf. [134]). A schematic diagram of such a shock absorber is shown in Figure 4 (left). In contrast to conventional shock absorbers, where the fluid chambers are filled with hydraulic oils, both in the compression and in the rebound stage, ERF shock absorbers have a much wider characteristics (damper force as a function of the velocity of the piston; cf. Figure 3 (right)). Therefore, ERF shock absorbers offer the best compromise between safety and comfort for a wide spectrum of road conditions.

Figure: Electrorheological shock absorber: schematic diagram (left) and B-Spline representation of the inlet and outlet boundaries of the right part of the fluid chamber (right)
\begin{figure}\vbox{\hbox to\textwidth{\hspace{2.5cm} \psfig{figure=Damper1.ps,w... ...m,height=7.0cm}\hfill}} \vspace{0.2cm} \hspace{2mm} \vspace{2mm} \end{figure}

The performance of the shock absorber does not only depend on the applied voltage and the velocity of the piston, but also on the geometry of the device. In particular, the geometry of the inlet and outlet boundaries of the ducts play a decisive role. In extreme cases, cavitation due to high pressure variations may occur which negatively effects the damper characteristics. Therefore, given a prescribed pressure profile $ p_{d} $, the optimization issue is to design the geometry in such a way that pressure variations are minimized. Due to axisymmetry, the computational domain $ \Omega $ reduces to the right part of the fluid chamber. The inlet and outlet boundaries are represented by B-splines using the de Boor control points $ \alpha = (\alpha_{1},...,\alpha_{m}) $ as design variables (cf. Figure 4 (right)). Consequently, the computational domain depends on the choice of the design variables, i.e., $ \Omega = \Omega(\alpha) $.
In the stationary case, the fluid flow is described by the equations

$\displaystyle - \nabla \cdot$$\displaystyle \mbox {\boldmath$\sigma $}({\bf u})$ $\displaystyle =$ $\displaystyle {\bf f}$   in$\displaystyle \Omega(\alpha) ,$ (12)
$\displaystyle \nabla \cdot {\bf u}$ $\displaystyle =$ 0   in$\displaystyle \Omega(\alpha)$ (13)

along with appropriate boundary conditions. Here, $ {\bf u}= (u_{1},u_{2}) $ is the velocity vector, $ \mbox {\boldmath $\sigma $}$ refers to the stress tensor and $ {\bf f} $ describes exterior forces acting on the fluid. The stress tensor $ \mbox {\boldmath $\sigma $}$ is related to the rate of deformation tensor $ {\bf D}({\bf u}) $ with $ ({\bf D}({\bf u}))_{ij} := (\partial u_{i}/\partial x_{j} + \partial u_{j}/\partial x_{i})/2 , 1 \le i,j \le 2, $ by a constitutive equation where the electric field $ {\bf E}$ enters as a parameter. Most constitutive equations are of extended Bingham fluid type [71,190] and use some sort of power law dependence [195]. Here, we have used a constitutive equation developed in [135]

$\displaystyle \mbox {\boldmath$\sigma $} = -p {\bf I} +2 varphi(I({\bf u}),\vert {\bf E}\vert,\mu({\bf u},{\bf E})) {\bf D}({\bf u}) ,$ (14)

where $ p $ is the pressure and $ \varphi$ stands for a viscosity function depending on the shear rate $ I({\bf u})$, the electric field strength $ \vert {\bf E}\vert$, and the angle $ \mu({\bf u},{\bf E}) $ between the velocity field $ {\bf u}$ and the electric field $ {\bf E}$.
The shape optimization problem can be stated as follows

$\displaystyle \inf\limits_{{\bf u},p,\alpha} J({\bf u},p,\alpha) , \J({\bf u}... ...pha) := \int\limits_{\Omega(\alpha)} \vert \nabla p - \nabla p_{d} \vert^{2} dx$ (15)

subject to the state equations (12),(13) and the inequality constraints

$\displaystyle \alpha_{i}^{min} \le \alpha_{i} \le \alpha_{i}^{max} \quad , \quad 1 \le i \le m$ (16)

on the design variables which are motivated by technological constraints on the shape of the inlet and outlet boundaries.
Assuming a known angle $ \mu({\bf u},{\bf E}) $ between the velocity vector and the exterior electric field, the flow model (12)-(14) has been discretized in space by $ Q_2-P_1 $ Taylor-Hood elements with respect to a simplicial triangulation of a reference computational domain $ \Omega(\hat{\alpha}), $ where $ \hat{\alpha} \in \textbf{\rm l\hskip-0.5mm R}^{m} $ is chosen within the set of feasible design parameters such that $ \Omega(\alpha) = \Phi(\Omega(\hat{\alpha})) $. The resulting nonlinear programming problem has been solved by the primal-dual Newton interior point method featuring an augmented Lagrangian approach [94] for the solution of the discretized state equations within the optimization loop (for details see [134,136]).

Figure: Optimized inlet boundary (left) and zoom (right) (the red line marks the final optimized configuration)
\begin{figure}\vbox{\hbox to\textwidth{ \put(3.0,0.0){\psfig{figure=OptOutflow.e... ...idth=4.0cm,height=3.5cm}} }} \vspace{0.2cm} \hspace{2mm}\vspace{2mm}\end{figure}

We note that the differences between the initial configuration and the optimized configurations are small, but give rise to reductions in the pressure variations between 10 and 20 %, depending on the operating conditions.

2.2.4 Modeling and simulation of active biochips

2.2.4.1 Piezoelectrically actuated acoustic surface waves

Mathematical models for SAW biochips are based on the linearized equations of piezoelectricity as given by

$\displaystyle \rho u_{i,tt} - c_{ijkl} u_{k,lj} - e_{kij} \Phi_{,kj} $ $\displaystyle =$ $\displaystyle 0 ,$  
$\displaystyle e_{ikl} u_{k,li} - \varepsilon_{ij} \Phi_{,ji} $ $\displaystyle =$ $\displaystyle 0 ,$  

where the $ u_i$ denote the components of the mechanical displacement, $ \Phi$ stands for the electric potential, and $ c_{ijkl}$, $ e_{kij}$ and $ \epsilon_{ij}$ refer to the components of the material tensors.
In a functional analytic setting, the piezoelectric equations can be stated as the following system of operator equations
$\displaystyle \rho {\bf u}_{,tt} + \mathcal{A} {\bf u} + \mathcal{B} \Phi $ $\displaystyle =$ $\displaystyle {\bf f} ,$    in $\displaystyle \mathcal{V}^{*} ,$  
$\displaystyle \mathcal{B}^{*} {\bf u} - \mathcal{C} \Phi $ $\displaystyle =$ $\displaystyle g ,$    in $\displaystyle \mathcal{W}^{*} ,$  

where $ \mathcal{V}$ and $ \mathcal{W}$ are appropriately chosen Hilbert spaces with $ \mathcal{V}^{*}$ and $ \mathcal{W}^{*}$ denoting their duals. The operators $ \mathcal{A}$ and $ \mathcal{C}$ inherit symmetry and ellipticity properties directly from the corresponding properties of $ {\bf c}$ and $ \mbox {\boldmath $\varepsilon $}$. Elimination of $ \Phi$ results in the Schur complement system

$\displaystyle \rho {\bf u}_{,tt} + \mathcal{S} {\bf u} = {\bf F} ,$    in $\displaystyle \mathcal{V}^{*} . $

The operator $ \mathcal{S}$ also plays a dominant role in the time harmonic approach,

$\displaystyle {\bf u}({\bf x}, t) = $   Re$\displaystyle \left( {\bf u}( {\bf x}) \exp(- i \omega t) \right) , $

resulting in the time harmonic operator relation

$\displaystyle \mathcal{S} {\bf u} - \rho \omega^2 {\bf I}{\bf u} = {\bf F} ,$    in $\displaystyle \mathcal{V}^{*} ,$

where $ \omega $ stands for the frequency of the excitation. The Riesz-Schauer theory provides conditions for the solvability of this equation [,]. In particular, these equations admit unique solutions except for countably many generalized real eigenvalues $ \rho \omega^2$ of $ \mathcal{S}$ that do not accumulate in $ \textbf{\rm l\hskip-0.5mm R}$. Even if $ \rho \omega^2$ is an eigenvalue of $ \mathcal{S}$, the solvability is guaranteed under certain conditions on the boundary data. However, in these cases the solutions are in general not unique.
After a finite element discretization by Lagrangian type finite elements with respect to a simplicial triangulation of the computational domain, we are led to the algebraic saddle point problem
$\displaystyle {\bf A}{\bf u} + {\bf B}\Phi $ $\displaystyle =$ $\displaystyle {\bf f} ,$  
$\displaystyle {\bf B}^T {\bf u} - {\bf C}\Phi $ $\displaystyle =$ $\displaystyle {\bf g} .$  

This saddle point problem is numerically solved by preconditioned Krylov subspace iterations featuring blockdiagonal multilevel preconditioners obtained with respect to a hierarchy of finite element meshes (cf. [81]).

For an SAW chip based on lithium niobate (LiNbO$ _3$) of length $ 1.2$mm and height $ 0.6$mm and with an operating frequency of $ f=\frac{\omega}{2\pi} = 100MHz$ Figure 6 displays the amplitudes of the electric potential which is in the range of nanometers.

Figure 6: Electric potential wave
\includegraphics[width=12.5cm,height=7.0cm]{100MHzLi4refpot2.eps}

The SAWs are strictly confined to the surface of the substrate. Their penetration depth into the piezoelectric material is in the range of one wavelength. The velocity of the SAW is independent of the applied frequency. In the case of $ \tt {YXl } 128^\circ$ LiNbO$ _3$, the SAW velocity is given by $ v=3992 \frac{m}{s}$, cf. []. Thus, for an excitation at the frequency $ f=100$MHz the theoretical wavelength of the SAW is given by $ l = \frac{v}{f}\approx 40 \mu m$. The simulations show the same wavelength for the SAW.
The excitation of an IDT on the surface of a piezoelectric material leads to the generation of bulk acoustic waves (BAWs) as well as surface acoustic waves. These bulk waves can also be observed in Figure 6. Technologically, they are desirably employed in solid-state circuits []. However, for SAW devices their presence is unwanted, since the interference of BAWs with SAWs can lead to a complete loss of functionality of the device. The used approach is sufficiently general to simulate every kind of piezoelectric resonator. In Figure 7 we have used an $ \tt {YXl} 38$° cut of LiNbO$ _3$ to generate a strong bulk acoustic wave at frequency $ f=200MHz$.

Figure 7: Bulk wave excitation
\includegraphics[width=12.5cm,height=7.5cm]{bulk.eps}

Rayleigh surface waves characteristically show an elliptical displacement, i.e., the displacements in the $ x_1$- and $ x_2$-direction are $ 90$° out of phase with each other. Additionally, the amplitude of the surface displacement in the $ x_2$-direction is larger than the one along the SAW propagation axis $ x_1$. The simulations reveal these properties (see Figures 8 and 9). In Figure 8, the displacements in the $ x_1$- and $ x_2$-direction for a certain surface area are depicted. The $ x_2$-displacements are flipped vertically for easier comparability.

Figure 8: Phase shift of flipped displacements
\includegraphics[width=12.5cm,height=2.5cm]{phaseshift.eps}

In Figure 9 a certain surface area is magnified and the vectors indicate the surface displacements.

Figure 9: Displacement vectors for the SAW
\includegraphics[width=12.5cm,height=7.0cm]{dampedvec.eps}



Modeling and simulation of microfluidic flows on biochips

The modeling of the micro-fluidic flow is based on the compressible Navier-Stokes equations (see [170,42]). In the model, compressible and non-linear effects are the driving force of the resulting flow. We consider the following situation: We start with a given equilibrium position $ \Omega $ of the flow field together with a given displacement $ {\boldsymbol{u}}\colon\partial\Omega\times\textbf{\rm l\hskip-0.5mm R}^+\to\textbf{\rm l\hskip-0.5mm R}^d$ describing the actual position of the time-dependent domain $ \Omega(t)$, occupied by the fluid. The displacement $ {\boldsymbol{u}}$ will be given by the modeling of active biochips described in the previous section. In the presented approach, we neglect the damping effect of the fluid on the solution of the piezoelectric equations.

We are now looking for the velocity $ {\boldsymbol{v}}$, pressure $ p $ and density $ \rho$ such that

\begin{equation*}\begin{aligned}\rho\left(\frac{\partial{\boldsymbol{v}}}{\parti... ...boldsymbol{v}}) &= 0\end{aligned}\qquad\text{in } \Omega(t), t>0\end{equation*}

where $ \eta$ and $ \zeta$ are the coefficients of standard and bulk viscosity, see e.g. [170], and

$\displaystyle {\boldsymbol{v}}({\boldsymbol{x}}+ {\boldsymbol{u}}({\boldsymbol{x}}, t), t) = \frac{\partial{\boldsymbol{u}}}{\partial t}({\boldsymbol{x}}, t)$   on $\displaystyle \partial\Omega.$ (17)

Note, that due to the shift $ {\boldsymbol{x}}\to {\boldsymbol{x}}+ {\boldsymbol{u}}({\boldsymbol{x}},t)$ the boundary condition (18) on the fixed boundary $ \partial\Omega$ turns into a boundary condition on $ \partial\Omega(t)$ and states that the velocity of flow field is in equilibrium with the velocity of the moving boundary. Hereafter, $ \eta$ and $ \zeta$ denote the normal respectively the bulk viscosity.

In order to close the system we prescribe initial values

$\displaystyle {\boldsymbol{v}}= 0,\quad p = p_0$   for $\displaystyle t=0 $

together with the constitutive equation

$\displaystyle p=p(\rho) = a\rho^\gamma$   with $\displaystyle a,\gamma>0$ (18)

coming from the barotropic (or isentropic) approximation, see e.g. [64]. It states that the entropy per unit mass of the fluid remains constant in space and time.

Model derivation

The Navier-Stokes system given above is not solved directly, mainly because the underlying physical processes take place on extremely different time scales. The elastic displacement of the substrate traversed by SAWs is a process with a time scale of nanoseconds, while the resulting flow due to acoustic streaming reaches an equilibrium on a time scale of milliseconds.

We use a technique from physics known as approximation theory to derive the model. First, we consider an expansion of the unknowns $ {\boldsymbol{v}}$, $ p $ and $ \rho$ in a parameter $ \epsilon$, to be chosen appropriately:

$\displaystyle p$ $\displaystyle = p_0+\epsilon p^\prime+\epsilon^2 p^{\prime\prime}+\ensuremath{\mathcal{O}\!\left(\epsilon^3\right)},$  
$\displaystyle \rho$ $\displaystyle =\rho_0+\epsilon\rho^\prime+\epsilon^2\rho^{\prime\prime}+\ensuremath{\mathcal{O}\!\left(\epsilon^3\right)},$  
$\displaystyle {\boldsymbol{v}}$ $\displaystyle = 0 + \epsilon {\boldsymbol{v}}^\prime + \epsilon^2 {\boldsymbol{v}}^{\prime\prime} + \ensuremath{\mathcal{O}\!\left(\epsilon^3\right)}.$  

Define $ p_1:=\epsilon p^\prime$, $ p_2:=\epsilon^2 p^{\prime\prime}$, and analogously $ \rho$, $ {\boldsymbol{v}}$. We assume that all terms with index 0 are known and constant in time and space. Choosing $ \epsilon$ as the maximal displacement of the elastic walls, we obtain in a first step a set of first order equations by collecting all terms of order $ \mathcal{O}(\epsilon)$:

$\displaystyle \begin{aligned} \rho_0\frac{\partial{\boldsymbol{v}}_1}{\partial ... ...\rho_0\nabla\cdot {\boldsymbol{v}}_1 &= 0,\\ p_1 &= c_0^2 \rho_1 \end{aligned}$   in $\displaystyle \Omega, t>0 $

together with the boundary condition

$\displaystyle {\boldsymbol{v}}_1({\boldsymbol{x}}, t) = \frac{\partial{\boldsymbol{u}}}{\partial t}({\boldsymbol{x}}, t) $

on the fixed boundary $ \partial\Omega$ for $ t>0$ and initial conditions $ {\boldsymbol{v}}_1 = 0$, $ p_1 = 0$ for $ t=0$. The constant $ c_0^2:=a\gamma\rho_0^{\gamma-1}$ represents the small signal sound speed in the fluid. This linear system describes the propagation of damped acoustic waves in the fluid.

In a second step, we then collect all second order terms while considering the solution of the first problem as given data:

$\displaystyle \begin{aligned} \rho_0\frac{\partial{\boldsymbol{v}}_2}{\partial ... ...cdot{\boldsymbol{v}}_2 &= - \nabla\cdot(\rho_1{\boldsymbol{v}}_1) \end{aligned}$   in $\displaystyle \Omega, t>0 $

together with the boundary condition

$\displaystyle {\boldsymbol{v}}_2 = -[\nabla{\boldsymbol{v}}_1]{\boldsymbol{u}}$   on $\displaystyle \partial\Omega, t>0 $

and initial conditions $ {\boldsymbol{v}}_2 = 0$, $ p_2 = 0$ for $ t=0$.

Contrary to the first order acoustic equations, we expect second order velocity and pressure to have nonzero mean value in time, as the data of the problem exhibit nonzero mean values. This is the mathematical background of the effect known in physics as acoustic streaming, which describes the creation of stationary flow patterns in reaction to a high frequency acoustic wave.

We now assume that the excitation $ {\boldsymbol{u}}$ is a harmonic process with period $ T:=2\pi/\omega$. This period, being on the order of nanoseconds, is small compared to the time scale of the relaxation processes in acoustic streaming. Applying the time-averaging operator

$\displaystyle \langle a\rangle := \frac{1}{T}\int_{t_0}^{t_0+T}a dt. $

to the set of second order equations, we arrive at

$\displaystyle \begin{aligned} - \eta\Delta {\boldsymbol{v}}_2 - \left(\zeta + \... ... \bigl\langle - \nabla\cdot(\rho_1{\boldsymbol{v}}_1)\bigr\rangle \end{aligned}$   in $\displaystyle \Omega $

and

$\displaystyle {\boldsymbol{v}}_2 = -\left\langle[\nabla{\boldsymbol{v}}_1]{\boldsymbol{u}}\right\rangle$   on $\displaystyle \partial\Omega.$

In the averaging process, we have also dropped the terms with temporal derivatives as we are only interested in the stationary equilibrium state observed in experiments.

Note that the boundary condition $ -\left\langle[\nabla{\boldsymbol{v}}_1]{\boldsymbol{u}}\right\rangle$ and the right hand side $ \bigl\langle - \nabla\cdot(\rho_1{\boldsymbol{v}}_1)\bigr\rangle$ are compatible, since

$\displaystyle \int_{\partial\Omega}{\boldsymbol{v}}_2\cdot\nu do$ $\displaystyle = - \int_{\partial\Omega} \left\langle[\nabla{\boldsymbol{v}}_1]{\boldsymbol{u}}\right\rangle\cdot\nu do$  
  $\displaystyle = \int_{\partial\Omega}\left\langle\frac{-1}{\rho_0}\rho_1{\bolds... ...a\times({\boldsymbol{\zeta}}_1\times{\boldsymbol{v}}_1)\right\rangle\cdot\nu do$  
  $\displaystyle = -\int_\Omega\frac{1}{\rho_0}\langle\nabla\cdot(\rho_1{\boldsymbol{v}}_1)\rangle dx = \int_\Omega\nabla\cdot{\boldsymbol{v}}_2dx.$  

Here $ {\boldsymbol{\zeta}}_1:=\int{\boldsymbol{v}}_1dt$ is an extension of $ {\boldsymbol{u}}$ to the interior of $ \Omega $.



Subproblem 1: Simulation of the damped acoustic field

The mathematical structure of the acoustic problem (after suitable scaling) is

\begin{equation*}\begin{aligned}\mu_v\frac{\partial{\boldsymbol{v}}}{\partial t}... ...{v}}(\cdot,0) = p(\cdot,0) &= 0 && \text{in }\Omega.\end{aligned}\end{equation*}

with the constants $ \mu_v,\mu_p,\nu_1>0,\nu_2>0$. The parameters $ \nu_i$ represent inverse Reynolds numbers for the problem, while $ \mu_p$ is proportional to $ 1/\rho_0c_0^2$. The boundary boundary condition is $ {\boldsymbol{g}}(t) = \frac{\partial{\boldsymbol{u}}}{\partial t}$.

This system is discretized by finite elements in space and the implicit Euler or Crank-Nicholson scheme in time. In the $ n$-th time-step this then results in a linear system of the form

$\displaystyle \begin{pmatrix}\ensuremath{\boldsymbol{A}}& \ensuremath{\boldsymb... ...rix} = \begin{pmatrix}\ensuremath{\boldsymbol{F}}^{n+1} G^{n+1} \end{pmatrix}$ (19)

for the discrete velocity $ \ensuremath{\boldsymbol{V}}^{n+1}$ and discrete pressure $ P^{n+1}$. The matrix $ \ensuremath{\boldsymbol{A}}$ corresponds to the linear elliptic operator in the velocity space, $ \ensuremath{\boldsymbol{B}}$ is the discretization of $ \nabla$, and $ \ensuremath{\boldsymbol{C}}$ is a scaled mass matrix in the pressure space. The right-hand-side of the system contains information about the solution of the old time step and Dirichlet boundary values. The corresponding Schur-complement operator $ \ensuremath{\boldsymbol{S}}:=\ensuremath{\boldsymbol{B}}^*\ensuremath{\boldsymbol{A}}^{-1}\ensuremath{\boldsymbol{B}}+\ensuremath{\boldsymbol{C}}$ is continuous and coercive on the discrete pressure space and thus (21) can efficiently be solved by a preconditioned conjugate gradient method.

As in the modeling of the piezoelectric material above, we also have the possibility of assuming a time harmonic dependency of the solution if given boundary data is time harmonic. In this situation we make the ansatz

\begin{displaymath}\begin{split}{\boldsymbol{v}}({\boldsymbol{x}},t) &= \widehat... ...ga t + \widehat{p}_s({\boldsymbol{x}})\sin\omega t, \end{split}\end{displaymath} (20)

where the frequency $ \omega>0$ corresponds to the governing frequency of the boundary datum $ {\boldsymbol{g}}({\boldsymbol{x}},t)$. With this ansatz we arrive at the following set of equations in the new unknowns $ \widehat{{\boldsymbol{v}}}_c$, $ \widehat{{\boldsymbol{v}}}_s$, $ \widehat{p}_c$, and $ \widehat{p}_s$:

$\displaystyle \begin{aligned} \omega\mu_v\widehat{{\boldsymbol{v}}}_s - \nu_1\D... ...\mu_p\widehat{p}_c + \nabla\cdot\widehat{{\boldsymbol{v}}}_s &= 0 \end{aligned}$

in $ \Omega $ together with the boundary conditions

$\displaystyle \widehat{{\boldsymbol{v}}}_c = \widehat{{\boldsymbol{g}}}_c$   and$\displaystyle \quad \widehat{{\boldsymbol{v}}}_s = \widehat{{\boldsymbol{g}}}_s$   on $\displaystyle \partial\Omega.$

Denoting by $ \ensuremath{\boldsymbol{M}}_{{\boldsymbol{v}}}$ the mass matrix in the velocity space, by $ \ensuremath{\boldsymbol{M}}_{p}$ the mass matrix in the pressure space and in capitals the discrete quantities, the finite element approximation of the time-harmonic equations leads to the linear system of the form

$\displaystyle \begin{bmatrix} \ensuremath{\boldsymbol{A}}& \omega\mu_v \ensurem... ... \widehat{{\boldsymbol{F}}}_s -\widehat{G}_c -\widehat{G}_s \end{bmatrix}$

where the right hand side comes from the non-homogeneous Dirichlet boundary conditions for the velocity. Equivalently, we can write this system in symmetric indefinite form

$\displaystyle \begin{bmatrix} -\omega\mu_v \ensuremath{\boldsymbol{M}}_{{\bolds... ... \widehat{{\boldsymbol{F}}}_c -\widehat{G}_s -\widehat{G}_c \end{bmatrix}. $

There is no significant difference concerning the numerical effort needed to solve the discretized system. For the algorithm, we chose the symmetric form and used a preconditioned CG method, which converged for the cases studied.

Concerning the difference between the time discretisation and the harmonic ansatz just presented, the numerical effort to solve these is also comparable. The harmonic ansatz has the bonus of yielding precise values for the right hand sides and boundary values of the acoustic streaming problem treated below. However, it is naturally not applicable to any problem involving nonharmonic driving forces, e.g. a SAW device operated in pulsed mode.



Subproblem 2: Simulation of Acoustic Streaming

The problem of acoustic streaming has the mathematical form

\begin{equation*} \begin{aligned} - \nu_1\Delta {\boldsymbol{v}}- \nu_2\nabla(\... ...{v}}&= {\boldsymbol{h}}&& \text{on }\partial\Omega, \end{aligned}\end{equation*}

with the above constants $ \nu_1>0$, $ \nu_2>2$, boundary conditions and right hand sides

$\displaystyle {\boldsymbol{h}}$ $\displaystyle = -\left\langle[\nabla{\boldsymbol{v}}_1]{\boldsymbol{u}}\right\rangle,$  
$\displaystyle {\boldsymbol{f}}_v$ $\displaystyle = -\beta_v\Bigl\langle(\nabla\cdot{\boldsymbol{v}}_1){\boldsymbol{v}}_1 + [\nabla{\boldsymbol{v}}_1]{\boldsymbol{v}}_1\Bigr\rangle,$  
$\displaystyle f_p$ $\displaystyle = -\beta_p\langle\nabla\cdot(p_1{\boldsymbol{v}}_1)\rangle.$  

The constants $ \beta_v,\beta_p>0$ arise from scaling.

This is a Stokes system with non-homogeneous divergence constraint and can efficiently be solved using a stable finite element discretization such as the Taylor-Hood element.



Simulation results

In this section we present numerical results for a 2d problem. The solvers for the damped acoustic field and for the acoustic streaming are also implemented within ALBERTA. We have used following data:

  • The underlying geometry is a unit square of side length $ 1$mm.
  • Excitation occurs through an SAW running from left to right along the lower edge, frequency is $ 100$ MHz, substrate material is LiNiO$ _3$.
  • Fluid is water at $ 20^\circ$C.
  • FEM-discretization using triangular continuous Lagrange-elements (quadratic for velocity, linear for pressure, i.e. the Taylor-Hood element).
  • CPU time: $ \approx 3$ hours

Figure 10: Snapshot of the oscillating pressure field $ p_1$ at time $ 5\cdot 10^{-7}s$, corresponding to $ 50$ periods.
\includegraphics[height=8cm]{koester_p1_1000}

Figure 11: Effective force for acoustic streaming $ {\boldsymbol {f}}_v$.
\includegraphics[height=10cm]{koester_fv}

Figure 12: Resulting acoustic streaming velocity field $ {\boldsymbol {v}}_2$.
\includegraphics[height=10cm]{koester_v2}



Simulation results

The problem was discretised by the Finite Element Method using the ALBERTA package developed as a joint project with Kunibert Siebert, Alfred Schmidt, and others. The data of the first 2D test problem:

  • Geometry: parallelogram with lower edge length and height $ 1$mm. The angle of the left wall to the horizontal axis is $ 35^\circ$.
  • Excitation through a SAW along the lower edge, frequency is $ 100$ MHz, substrate material is LiNiO$ _3$.
  • Fluid is water at $ 20^\circ$C.
  • FEM-discretisation using triangular Lagrange-elements (quadratic for velocity, linear for pressure).
  • CPU time: $ \approx 1$ hour

Figure 13: Mesh used for the calculation, the distance scale is in $ mm$. The bottom edge is strongly refined.
\includegraphics[height=7cm]{koester_mesh1}

Figure 14: Snapshot of the acoustic velocity field $ {\bf v}_1$ at the bottom left corner, velocity scale on the left is in $ mm/s$. Obviously visible is the strong damping caused by the liquid.
\includegraphics[height=10cm]{koester_exp1_v1}



Experimental work

The acoustically driven micro- and nanofluidic has over the last four years been developed in the research group of Prof. Wixforth. Based on this technology, many different applications have been devised which even led to the foundation of a start up company (Advalytix, Brunnthal). Since then , a close cooperation between the research group and the company has been established, including many more cooperation partners working in this field. Mostly, we have been concentrating on so-called planar microfluidic chips, an example of which is given in section 2.1.2. of this proposal. Here, surface acoustic waves (SAW) are used to interáct with small amounts of fluids at the surface of a planar chip. Chemical instrumentalization of parts of the chip's surface is employed to laterally modulate the wetting properties of the surface, hence providing "virtual tracks" for the fluid, if actuated by the SAW. A more or less detailed understanding of the physics behind the interaction between the fluid and the SAW has been reached [98,82]. Based on the planar microfluidic technology, quite complex programmable biochips have been demonstrated. The most advanced chip, so far, is probably a system, on which a sub-microliter polymerase chain reaction (PCR) has been successfully demonstrated [97].

Aims of the project

The goal of the project is

  • the development, analysis, and implementation of adaptive multilevel based algorithms within the class of all-at-once methods for the efficient and reliable numerical solution of PDE constrained optimization problems and
  • the application of these algorithms to the optimal design of microfluidic biochips used as labs-on-a-chip for biomedical purposes.
The theoretical part of the project will focus on
  • path following techniques in primal-dual Newton interior-point methods and
  • Newton multigrid and fully nonlinear multigrid methods based on primal-dual active set strategies,
  • automatic time-stepping and adaptivity in space based on the goal oriented weighted dual approach.
The emphasis in the application to the structural optimization of microfluidic biochips is on
  • optimization of the channel geometry on the chip surface,
  • optimal positioning of the interdigital transducers, optimization of its apertures, and optimization of the frequency of the piezoelectrically agitated surface acoustic waves,
both in order to achieve maximal pumping rates for the surface acoustic wave driven fluidic transport. Further issues are the
  • optimization of the geometry of the reaction chamber,
  • optimal positioning and design of the capillary barriers at the inflow/outflow boundaries of the reaction chamber
to provide the filling of the chamber in minimal time under no bubbly flow conditions.

Methods

Primal-dual Newton interior-point methods

The interior-point aspect of primal-dual Newton interior-point methods is that classical barrier methods can be used to transform constrained problems to unconstrained ones. They typically give rise to parameterized families of approximate subproblems. In particular, coupling the inequality constraints by logarithmic barrier functions gives rise to

$\displaystyle B^{(\beta)}(x_{1},x_{2}) := f(x_{1},x_{2}) - \beta \sum\limits_{i=1}^{p}$   log$\displaystyle (g_{i}(x_{1},x_{2})) .$ (21)

Here, $ \beta > 0 $ stands for the barrier parameter [75,78,228]. We are thus led to the following parameterized family of minimization subproblems
\begin{subequations} \begin{align} &\text{minimize} \quad B^{(\beta)}(x_{1},x... ...ubject to:} \quad h(x_{1},x_{2}) = 0 \ . \end{align}%EQNO:21 \end{subequations}
Under standard assumptions on the NLP, for sufficiently small $ \beta > 0 $ there exist local minima $ (x_{1}^{(\beta)},x_{2}^{(\beta)}) $ of (24a),(24b) which converge along centered pathes to an isolated local minimum of (2a)-(2c) as $ \beta \rightarrow 0 $ (cf., e.g., [75]).
The primal-dual aspect is to take care of the equality constraint (24b) by Lagrangian multipliers $ \lambda = (\lambda_{1},\lambda_{2})^{T} , \lambda_{1} \in \textbf{\rm l\hskip-0.5mm R}^{n_{1}}, \lambda_{2} \in \textbf{\rm l\hskip-0.5mm R}^{m-n_{1}}, $ (dual variables). Hence, we obtain the saddle point problem

$\displaystyle \min\limits_{(x_{1},x_{2}) \in \textbf{\rm l\hskip-0.5mm R}^{n_{1... ...da \in \textbf{\rm l\hskip-0.5mm R}^{m}} L^{(\beta)}(x_{1},x_{2},\lambda) .$ (22)

Here, $ L^{\beta} $ stands for the Lagrangian

$\displaystyle L^{(\beta)}(x_{1},x_{2},\lambda) := B^{(\beta)}(x_{1},x_{2}) + \lambda^{T} h(x_{1},x_{2}) .$ (23)

Considering the KKT-conditions for (25), introducing the slack variable $ z \in \textbf{\rm l\hskip-0.5mm R}^{q} $ (approximate complementarity) according to

$\displaystyle z_{i} := \frac{\beta}{g_{i}(x_{1},x_{2})} \quad , \quad 1 \le i \le p , $

the Newton aspect is to apply Newton's method with respect to $ \varphi = (x_{1},x_{2},\lambda,z) $. We thus arrive at the primal-dual system

$\displaystyle K \delta_{\varphi} = - F^{(\beta)}(\varphi)$ (24)

for the Newton increments $ \delta_{\varphi} := (\delta_{x_{1}},\delta_{x_{2}},\delta_{\lambda},\delta_{z})^{T} $. Here, the Hessian $ K $ is given by
$\displaystyle K = \left( \begin{array}{cccc} \nabla^{2}_{x_{1}x_{2}} L^{(\beta)... ... \nabla_{x_{1}} g & \nabla_{x_{2}} g & 0 & D_{z}^{-1} D_{g} \end{array} \right)$      

with $ D_{g} :=$   diag$ (g_{i}(x_{1},x_{2}))_{i=1}^{p} $ and $ D_{z} :=$   diag$ (z_{i})_{i=1}^{p} $, whereas the right-hand side is the residual with respect to the KKT conditions .
Since $ D_{g} , D_{z} $ are diagonal matrices, it is easy to perform block Gauss elimination of the approximate complementarity $ z $ which leads to the condensed primal-dual system in $ \delta_{\varphi^{(C)}} = (\delta_{x_{1}},\delta_{x_{2}},\delta_{\lambda})^{T} $

$\displaystyle K^{(C)} \delta_{\varphi^{(C)}} = - F^{(C)}(\varphi^{(C)})$ (25)

with the condensed Hessian

$\displaystyle K^{(C)} = \left( \begin{array}{ccc} \tilde{\nabla}^{2}_{x_{1}... ..._{2}} h)^{T} \nabla_{x_{1}} h & \nabla_{x_{2}} h & 0 \end{array} \right) .$ (26)

Here,

$\displaystyle \tilde{\nabla}^{2}_{x_{i}x_{j}} L^{(\beta)} = \ \nabla^{2}_{x_... ...bla_{x_{i}} g)^{T} D_{g}^{-1} D_{z} \nabla_{x_{j}} g , 1 \le i,j \le 2 . $

The right-hand side in (28) is given by
$\displaystyle F_{x_{i}}^{(C)}(\varphi^{(C)})$ $\displaystyle =$ $\displaystyle \nabla_{x_{i}} f + \nabla_{x_{i}} h - \beta (\nabla_{x_{i}} g)^{T} D_{g}^{-1} e ,$  
$\displaystyle F_{\lambda}^{(C)}(\varphi^{(C)})$ $\displaystyle =$ $\displaystyle h ,$   where$\displaystyle \ e = (1,...,1)^{T} \in \textbf{\rm l\hskip-0.5mm R}^{p} .$  

Multilevel iterative solvers for the condensed primal-dual system (28) have been considered in [,145,176] using a heuristically motivated decreasing sequence of barrier parameters and a step-length selection to ensure feasibility of the iterates. Following [85], global convergence issues have been addressed by means of a watchdog strategy relying on a hierarchy of two merit functions.
Our aim is to develop, analyze and implement theoretically supported multilevel path-following continuation methods with an adaptive trust-region approach featuring a multilevel predictor-corrector strategy. We note that for parameter dependent nonlinear PDEs and variational inequalities, multilevel continuation methods have been suggested in [14,99] and [139], but have not yet been studied in the context of barrier methods for PDE constrained optimization problems.
In particular, we will study the
  • sensitivity of the multigrid scheme with respect to the barrier parameters, since it is well-known that the Hessian gets ill-conditioned when the barrier parameter approaches zero (cf., e.g., [79]),
  • detection of negative curvature by using computationally inexpensive information from coarser grids [41],
  • appropriate specification of termination criteria for the multilevel iterative solver within an affine invariant framework (cf., e.g., [65,198,219]).


Primal-dual active set strategies

Active set strategies are iterative schemes for the solution of inequality constrained optimization problems, where at each iteration step the solution of an equality constrained subproblem is required by identifying a set of active constraints (cf., e.g., [183] and the references therein). In case of bilateral box constraints on the design variables

$\displaystyle \hspace*{1.6cm} c_{i} \le x_{2,i} \le d_{i} \quad , \quad i \in I := \{1,...,n_{2}\} , $

the active set $ \mathcal{A}$, i.e., the set where the inequality constraints are active, is defined according to
$\displaystyle \mathcal{A}$ $\displaystyle =$ $\displaystyle \mathcal{A}_{1} \bigcup \mathcal{A}_{2} , \hspace{1.2cm}$  
$\displaystyle \mathcal{A}_{1} := \{ i \in I \vert x_{2,i} = c_{i} \}$ $\displaystyle ,$ $\displaystyle \mathcal{A}_{2} := \{ i \in I \vert x_{2,i} = d_{i} \} .$  

The complement

$\displaystyle \mathcal{I} := I \setminus \mathcal{A} $

is said to be the set of inactive constraints. In this case, the KKT conditions associated with the saddle point problem take the form
\begin{subequations} \begin{align} \nabla_{x_{1}} L(x_{1},x_{2},\lambda,\mu) \ ... ...\ , \ \mu\vert _{\mathcal{I}} = 0 \ . \ \end{align}%EQNO:26 \end{subequations}
Here, $ \mu\vert _{S} $ denotes the restriction of $ \mu $ to $ S \subset I $.
The complementarity conditions (30d)-(30e) can be stated as

$\displaystyle \mu \in \partial I_{K}(x_{2}) ,$ (27)

where $ \partial I_{K} $ denotes the subdifferential of the indicator function $ I_{K} $ of the convex set
$\displaystyle K := \{ x_{2} \in \textbf{\rm l\hskip-0.5mm R}^{n_{2}} \vert c_{i} \le x_{2,i} \le d_{i} , 1 \le i \le n_{2} \} .$      

Following [22,117] and using the generalized Moreau-Yosida approximation

$\displaystyle (\partial I_{K})_{\sigma}(v) := \frac{1}{\sigma} ( I - ... ...}}) = ( (\partial I_{K})^{-1} + \ \sigma I)^{-1} , \sigma > 0 , $

of the indicator function $ I_{K} $, (31) can be replaced by the computationally more feasible condition

$\displaystyle \mu = \sigma \big[ x_{2} + \sigma^{-1} \mu - P_{K}(x_{2} + \sigma^{-1} \mu) \big] ,$ (28)

where $ P_{K} $ denotes the projection onto $ K $ as given by
$\displaystyle P_{K}(w) := \left\{ \begin{array}{r@{\quad,\quad}l} c_{i} & w_{i}... ...} \le d_{i} , \\ d_{i} & w_{i} > d_{i} \end{array} \right. , i \in I .$      

The primal-dual active set method proceeds as follows: Given startiterates $ x_{1}^{(0)},x_{2}^{(0)},\lambda^{(0)},\mu^{(0)} $, for $ \nu \ge 1 $ we determine the set $ \mathcal{A}^{(\nu)} $ of active constraints according to

$\displaystyle \mathcal{A}^{(\nu)} = \mathcal{A}_{1}^{(\nu)} \bigcup \mathcal{A}_{2}^{(\nu)} , $

$\displaystyle \mathcal{A}_{1}^{(\nu)} := \{ i \in I \vert x_{2,i}^{(\nu-1)}... ...n I \vert x_{2,i}^{(\nu-1)} + \sigma^{-1} \mu_{i}^{(\nu-1)} > d_{i} \} , $

and define $ \mathcal{I}^{(\nu)} := I \setminus \mathcal{A}^{(\nu)} $. We set

$\displaystyle x_{2,i}^{(\nu)} := c_{i} , i \in \mathcal{A}_{1}^{(\nu)} ... ...nu)} \quad , \quad \mu_{i}^{(\nu)} := 0 , i \in \mathcal{I}^{(\nu)} .$ (29)

Taking this into account, for $ \hat{x}_{2} := (x_{2,i})_{i \in \mathcal{I}^{(\nu)}} $ we define $ \hat{f} : \textbf{\rm l\hskip-0.5mm R}^{n_{1}} \times \textbf{\rm l\hskip-0.5m... ...at{x}_{2}) \longmapsto \hat{f}(x_{1},\hat{x}_{2}) := f(x_{1},x_{2}^{(\nu-1/2)})$, and $ \hat{h} : \textbf{\rm l\hskip-0.5mm R}^{n_{1}} \times \textbf{\rm l\hskip-0.5m... ...t{x}_{2}) \longmapsto \hat{h}(x_{1},\hat{x}_{2}) := h(x_{1},x_{2}^{(\nu-1/2)}) $, where $ n_{3} :=$   card$ \mathcal{I}^{(\nu)} $ and
\begin{displaymath}\hspace*{1.6cm} x_{2,i}^{(\nu-1/2)} := \left\{ \begin{arr... ...,i}^{(\nu)} & i \in \mathcal{A}^{(\nu)} \end{array} \right. .\end{displaymath}      

Then, we compute $ (x_{1}^{(\nu)},\hat{x}_{2}^{(\nu)},\lambda^{(\nu)}) \in \textbf{\rm l\hskip-0.... ...s \textbf{\rm l\hskip-0.5mm R}^{n_{3}} \times \textbf{\rm l\hskip-0.5mm R}^{m} $ as the solution of

$\displaystyle \min\limits_{(x_{1},\hat{x}_{2}) \in \textbf{\rm l\hskip-0.5mm R}... ...da \in \textbf{\rm l\hskip-0.5mm R}^{m}} \hat{L}(x_{1},\hat{x}_{2},\lambda) ,$ (30)

where the Lagrangian is given by
$\displaystyle \hspace*{1.0cm} \hat{L}(x_{1},\hat{x}_{2},\lambda) := \ \hat{f}(x_{1},\hat{x}_{2}) + \lambda^{T} \hat{h}(x_{1},\hat{x}_{2}) .$      

Finally, setting $ x_{2,i}^{(\nu)} := \hat{x}_{2,i}^{(\nu)} , i \in \mathcal{I}^{(\nu)}, $ the new iterate $ x_{2}^{(\nu)} $ for the design variable is completely specified.
The KKT conditions for (34) are given by
$\displaystyle \hspace*{0.2cm} \nabla_{x_{1}} \hat{L}(x_{1},\hat{x}_{2},\lambda)$ $\displaystyle =$ $\displaystyle \nabla_{x_{1}} \hat{f}(x_{1},\hat{x}_{2}) + \nabla_{x_{1}} \hat{h}(x_{1},\hat{x}_{2}) = 0 ,$  
$\displaystyle \hspace*{0.2cm} \nabla_{\hat{x}_{2}} \hat{L}(x_{1},\hat{x}_{2},\lambda)$ $\displaystyle =$ $\displaystyle \nabla_{\hat{x}_{2}} \hat{f}(x_{1},\hat{x}_{2}) + \nabla_{\hat{x}_{2}} \hat{h}(x_{1},\hat{x}_{2}) = 0 ,$  
$\displaystyle \hspace*{0.2cm} \nabla_{\lambda} \hat{L}(x_{1},\hat{x}_{2},\lambda)$ $\displaystyle =$ $\displaystyle \hat{h}(x_{1},\hat{x}_{2}) = 0 .$  

The iteration is stopped when the active set becomes stationary, i.e., $ \mathcal{A}^{(\nu)} = \mathcal{A}^{(\nu-1)} $.
Applying Newton's method to the KKT conditions, we arrive at a block-structured linear system of much the same form as (29).


In view of (32), the active set strategy can be interpreted as the discrete version of a semismooth Newton method in function space [118,210]. Based on this interpretation, it is the aim of the proposal to develop, analyze and implement Newton multigrid methods as well as full nonlinear multigrid methods.
Specific issues that will be addressed include the

  • synopsis of the continuous and discrete Newton methods with regard to a mesh independence principle (cf., e.g., [221] for elliptic PDEs),
  • specification of appropriate termination criteria for the inner multilevel iterative solvers within a Newton multilevel approach (cf., e.g., [66,67] for nonlinear elliptic PDEs),
  • design of suitable transfer operators ('fine-to-coarse' and 'coarse-to-fine') in a full nonlinear multilevel approach, when the transfer effects active and inactive nodal points at subsequent levels of the hierarchy of grids (cf., e.g., [133,165] in case of obstacle problems for elliptic PDEs).


Adaptive mesh refinement

In this project, we plan to investigate different target quantities and study their impact on the mesh adaptation and solution accuracy. We emphasize that the inequality constraints introduce a nonlinear and even nonsmooth behavior to the problem. Due to this fact there are several impacts on the (a posteriori) error estimator based mesh adaptivity:

  • For the inequality constrained variable, typically standard error estimates (based on variational equalities) do not apply. The presence of the constraints may result in a reduced approximation order.
  • The design of the error estimator for the mesh adaptation should yield a tool which reduces to a (sharp) estimator for the unconstrained problem, i.e., in the case where the inequality constraints are always inactive; otherwise, a strong overestimation may result in inefficient and even improper mesh adaptation.
  • The target quantities should be connected to the optimization problem. This may result in nonlinear and even nonsmooth functionals.
In this project, we will develop
  • residual type a posteriori error estimators providing local error terms for the global discretization errors in the state and design variables,
  • error estimators relying on the dual weighted approach involving problem oriented target quantities,
and use them as a basis for mesh refinement of the finite element meshes for the state and the design variables.
Since the optimality conditions for the parameterized barrier functions are related to the complementarity condition that holds true at optimality as the barrier parameter $ p $ converges to zero, the impact of the inequality constraints on the error estimation will be measured by means of the primal-dual central path associated with the primal-dual Newton interior point method which is used in the numerical solution. A similar remark applies to path-following primal-dual active set strategies.

Shape optimization of microfluidic biochips

A current trend in the design of active SAW-driven microfluidic biochips is the production of plastic chips with a network of channels and reservoirs that are placed on top of the SAW device. The channels and reservoirs are milled out of the plastic plate and a plastic foil is finally vapor deposited on top of it. Figure 15 shows such a plastic chip with two interconnected reservoirs featuring capillary barriers at its inflow/outflow boundaries for the purpose of a controlled filling process (see Figure 17). The chips can be used in various biomedical applications such as the early diagnosis of liver damage (e.g., hepatitis) by measuring enzyme activities in the blood serum of patients. The concentration of the enzyme is determined indirectly by measuring the optical density. For doing the biochemical analysis, certain reagents must be transported to the probe in a controlled way along the channels of the network to the reaction chamber. The transport of the fluid and the precise dosage are realized by the piezoelectrically actuated SAWs and by the design of the capillary barriers.
For the mathematical description of the transport process and the filling of the reservoirs, we will provide a two-phase flow model (air/liquid 1) and a multi-phase flow model (air/liquid 1/liquid 2/liquid 3). Here, liquid 1 is supposed to contain the reagents, whereas liquid 2 is carrying the probe and liquid 3 results from the mixing of liquid 1 and liquid 2. We will consider the optimization of the channel geometry and the optimal positioning of the interdigital transducers (IDT), the optimization of its apertures and of the frequency of the SAWs. At a later stage, we will be concerned with the shape optimization of the reservoir as well as with the optimal positioning of the capillary barriers and the optimization of its apertures which requires the use of Signorini type boundary conditions for modeling the stopping properties of the capillary barriers.
For the efficient numerical solution of these optimization problems, we will use the multilevel based all-at-once methods to be developed in this project.

Figure 15: Lab-on-a-chip (channel network with two reservoirs)
\includegraphics[width=6.5cm,height=6.5cm]{Chip_komplett.ps}

Optimization of the geometry of the channel

We assume the position of the IDTs and the frequency of the SAWs to be fixed. Hence, in a pre-processing step we are able to compute the displacement $ {\bf u}$ on the walls $ \Gamma_{SAW} \subset \partial \Omega $ traversed by the SAWs using the methods as described in 2.2.2.1. The walls of the channel are represented by a spline surface involving the product of cubic spline blending functions with respect to an array $ \mbox {\boldmath $\alpha $}:= (\alpha_{i_{1},i_{2}})_{i_{1}=1,i_{2}=1}^{n_{1},n_{2}} $ of control points serving as the design variables. Due to technical reasons, the design variables are subject to bilateral constraints $ \alpha_{i_{1},i_{2}}^{min} , \alpha_{i_{1},i_{2}}^{max} , 1 \le i_{\nu} \le n_{\nu} , 1 \le \nu \le 2 $. The computational domain $ \Omega $, which is the flow region inside the channel, thus depends on the choice of the design variables, i.e., $ \Omega = \Omega($$ \mbox {\boldmath $\alpha $}) $. The design objective is to achieve an optimal pumping rate. The state equation is given by the two-phase air/liquid 1 flow model.

Optimization of the position and aperture of the IDT and of the frequency of the SAW

In this case, we keep the geometry of the channel $ \Omega_{ch} $ fixed, but choose the position $ {\bf x}_{\bf IDT} $ of the IDT (e.g., monodirectional IDT in the channel or IDT at the channel corner), its aperture $ a $ and the frequency $ \omega $ of the SAW as design variables, i.e., $ \mbox {\boldmath $\alpha $}:= ({\bf x}_{\bf IDT},a,\omega) $. The design objective and the state equation are the same as before. As far as the inequality constraints are concerned, suitable bounds for the design variables $ \alpha_{4} = a $ and $ \alpha_{5} = \omega $ are $ a^{min} = 16 dBm , a^{max} = 28 dBm $ and $ \omega^{min} = 75 MHz , \omega^{max} = 190 MHz $.

Shape optimization of the reaction chamber

Figure 16 displays the fluid flow into the reservoir $ \Omega_{res} $ with two capillary barriers $ \Gamma_{cb}^{(\nu)} \subset \partial \Omega_{res} , 1 \le \nu \le 2, $ (see also Figure 17) at various stages of the filling process. In particular, Figure 16 (left) shows a snapshot of the first exit time of the (green) dye, whereas Figure 16 (middle) represents a snapshot of the first time instant where the (green) dye covers the outermost edges of the wall of the reaction chamber.

Figure 16: Filling of the reaction chamber
\includegraphics[width=13.5cm,height=4.2cm]{filling.ps}

Experimental evidence suggests that the diamond-like geometry of the reservoir is not optimal for rapid filling.
At a later stage of the project, we will consider the optimization of the geometry of the reservoir in order to achieve the filling of the chamber in minimal time $ t_{fil}^{*} $. The mathematical description of the filling process requires the use of a multi-phase flow model in not necessarily circular channels and reservoirs, where the capillary barriers can be modeled by Signorini type boundary conditions

$\displaystyle p_{1} - p_{2} - p_{crit}^{(\nu)} \ge 0 , {\bf n}\cdot {\bf v}\ge 0 , (p_{1} - p_{2} - p_{crit}^{(\nu)}) ({\bf n}\cdot {\bf v}) = 0$   a.e. on$\displaystyle \ \Gamma_{cb}^{(\nu)} . $

Here, $ p_{1} $ and $ p_{2} $ denote the pressures of two separated phases, and $ p_{crit}^{(\nu)} $ stands for a critical pressure, depending on the aperture $ a_{cb}^{(\nu)} $ of the capillary barriers, the surface tension $ \sigma $ and the contact angle $ \Theta $ between the liquid in phase 2 and the wall.
For fixed positions and apertures of the capillary barriers, the design objective is to optimize the geometry of the reservoir $ \Omega_{res} $ in such a way that the filling time $ t_{fil} $ is minimized under the additional constraint that no bubbly flow conditions occur. As before, the design variables are given by an array $ a_{i_{1},i_{2}} , 1 \le i_{\nu} \le n_{\nu} , 1 \le \nu \le 2, $ of control points for a spline surface representation of the wall of the reservoir.

Structural optimization of the capillary barriers

For a fixed geometry of the reservoir, another design objective is to determine the positions $ {\bf x}_{cb}^{(\nu)} $ and the apertures $ a_{cb}^{(\nu)} $ of the capillary barriers in order to minimize the filling time. The state equations are the same as before including Signorini type boundary conditions for the capillary barriers, whereas the design variables are given by $ \mbox {\boldmath $\alpha $}= ({\bf x}_{cb}^{(\nu)},a_{cb}^{(\nu)})_{\nu=1}^{2} $.

Figure 17: Capillary barriers
\includegraphics[height=4.5cm]{Kapillarstop2.ps} \includegraphics[height=4.5cm]{Kapillarstop3.ps}



Experimental work

In this project, we plan to investigate a combination of acoustic nanopumps and conventional microfluidic systems, i.e., three dimensional complex systems of channels and reactors, in which biological assays can be performed. Such 3d biochips usually consist of cheap substrate materials like, e.g., plastic, silicon, or glass. Deep etching, molding or hot embossing techniques are employed to fabricate the channels and reactor compartments in such 3d microfluidic chips. However, most of these chips rely on external pumping or integrated electro-osmotic actuation schemes, that cause additional costs and problems. Moreover, mixing in such devices is usually a difficult task because of the laminar flow at low Reynolds numbers. Here, we have recently successfully demonstrated [204] that the coupling of SAW into a 3d microfluidic chip leads to pronounced and efficient mixing, like in capillary gaps [98] or in open droplets [97]. Here, the SAW on a piezoelectric substrate is coupled into the chip substrate (plastic in the latter case), where it is converted into a wave in the substrate material, and then coupled into the fluid in the channel, where it induces acoustic streaming and mixing.
The coupling of SAW into different substrate materials and finally into the fluid in a 3d channel requires a detailed knowledge and understanding of the sound propagation in such "layered" systems. Not only diffraction and mode conversion has to be taken into account, also the complex shape of the channels and reactor compartments need to be considered.
A typical example of such a 3d biochip is given in Fig. 1, where we show a microfluidic reactor for a biological assay. It holds an inlet and an outlet (filled with blue colored water) and a reaction chamber (filled with yellow fluid). The inlet is coupled to an SAW pump on a piezoelelectric chip underneath the plastic chip (not visible in the figure), which - if activated - pumps the reagent (blue) into the reaction chamber. There, it is stirred by another SAW coming from below the yellowish region (Fig. 18).

Figure 18: 3d microfluidic chip coupled to a piezoelectric SAW chip residing below the plastic reactor
\includegraphics[height=4.5cm]{Dosierung1.ps} \includegraphics[height=4.5cm]{Dosierung2.ps}

In this figure, another important functional block of the investigated microfluidic system is already visible: To prevent premature mixing of the two reagents and to be able to meter a certain amount of "blue" reagent into the "yellow" reactor, a capillary barrier has been employed at the connection between the inlet and the reactor, and the outlet and the reactor. This capillary barrier should withhold the reagent 1 without the activation of the SAW pump, and should open under the influence of the SAW induced acoustic streaming. In that sense, it represents a pressure driven valve, which needs to be optimized for optimum functionality. This optimization is an important task to be solved in this proposal.
To be able to precisely meter a certain amount of fluid, employing the SAW pump, a detailed understanding of the flow profile, the velocity fields and hence the pressure in the channel is necessary.
As has already been shown, the acoustic coupling between SAW and fluid leads to a streaming velocity field, that results in an angle between the surface of the piezoelectric chip and the acoustic "jet" into the fluid. This "jet can be nicely visualized by using a bit of ink, as shown in Figure 19.

Figure 19: Layered microfluidic system as described in the text
\includegraphics[width=10.5cm,height=5.5cm]{Jet.eps}

An SAW is coupled from a piezoelectric substrate (LiNbO3) through a plastic substrate into the fluid. The resulting acoustic streaming under an angle is visualized by a little bit of ink in the water.
Recent experimental investigations have resulted in a modified chip layout, taking into account this "refraction" angle. It could be shown that an oblique inlet channel (tilted under exactly this angle) leads to a much higher streaming velocity than just coupling the ultrasound into a parallel channel. First experiments revealed an increase of the flow velocity by about a factor of five, even in a non-optimized channel geometry. Again, the optimization of the channel geometry with respect to the maximum flow velocity is an important scope of the present proposal.

Bibliography