Spatially Adaptive Stochastic Methods for Fluid-Structure Interactions Subject to Thermal Fluctuations in Domains with Complex Geometries
SSPATIALLY ADAPTIVE STOCHASTIC METHODS FORFLUID-STRUCTURE INTERACTIONS SUBJECT TO THERMALFLUCTUATIONS IN DOMAINS WITH COMPLEX GEOMETRIES
PAT PLUNKETT ∗ , JON HU † , CHRIS SIEFERT ‡ , AND
PAUL J. ATZBERGER § Abstract.
We develop stochastic mixed finite element methods for spatially adaptive simula-tions of fluid-structure interactions when subject to thermal fluctuations. To account for thermalfluctuations, we introduce a discrete fluctuation-dissipation balance condition to develop compati-ble stochastic driving fields for our discretization. We perform analysis that shows our condition issufficient to ensure results consistent with statistical mechanics. We show the Gibbs-Boltzmann dis-tribution is invariant under the stochastic dynamics of the semi-discretization. To generate efficientlythe required stochastic driving fields, we develop a Gibbs sampler based on iterative methods andmultigrid to generate fields with O ( N ) computational complexity. Our stochastic methods providean alternative to uniform discretizations on periodic domains that rely on Fast Fourier Transforms.To demonstrate in practice our stochastic computational methods, we investigate within channelgeometries having internal obstacles and no-slip walls how the mobility/diffusivity of particles de-pends on location. Our methods extend the applicability of fluctuating hydrodynamic approachesby allowing for spatially adaptive resolution of the mechanics and for domains that have complexgeometries relevant in many applications. Key words.
Stochastic Eulerian Langrangian Method, Immersed Boundary Method, Adap-tive Numerical Methods, Multigrid, Stochastic Numerical Methods, Stochastic Partial DifferentialEquations.
Introduction.
We develop general computational methods for applications in-volving the microscopic mechanics of spatially extended elastic bodies within a fluidthat are subjected to thermal fluctuations. Motivating applications include the studyof the microstructures of complex fluids [17], lipid bilayer membranes [28, 32, 48], andmicro-mechanical devices [37, 29]. Even in the deterministic setting, the mechanicsof fluid-structure interactions pose a number of difficult and long-standing challengesowing to the rich behaviors that can arise from the interplay of the fluid flow andelastic stresses of the microstructures [19, 42]. To obtain descriptions tractable foranalysis and simulations, approximations are often introduced into the fluid-structurecoupling. For deterministic systems, many spatially adaptive numerical methods havebeen developed for approximate fluid-structure interactions [25, 35, 26, 39, 2, 30]. Inthe presence of thermal fluctuations, additional challenges arise from the need to cap-ture in computational methods the appropriate propagation of fluctuations through-out the discretized system to obtain results consistent with statistical mechanics. Inpractice, challenges arise from the very different dissipative properties of the discreteoperators relative to their continuum differential counterparts. These issues haveimportant implications for how stochastic fluctuations should be handled in the dis-crete setting. Even when it is possible to formulate stochastic driving fields in awell-founded manner consistent with statistical mechanics, these Gaussian random ∗ University of California, Department of Mathematics † DOE Sandia ‡ DOE Sandia, Sandia National Laboratories is a multiprogram laboratory managed and oper-ated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for theU.S. Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000. § University of California, Department of Mathematics , Santa Barbara, CA 93106; e-mail:[email protected]; phone: 805-893-3239; Work supported by DOE CM4 and NSF CAREERGrant DMS-0956210. 1 a r X i v : . [ c ond - m a t . m e s - h a ll ] N ov P. PLUNKETT, P.J. ATZBERGER fields have often many degrees of freedom and non-trivial spatial correlations thatcan be difficult to sample without significant computational expense. Many finitedifference methods on uniform periodic meshes have been developed for fluctuatinghydrodynamics [9, 16, 41, 18, 7, 47, 6]. One of the main reasons that fluctuating hy-drodynamics is treated on uniform periodic domains is so that stochastic driving fieldscan be generated using Fast Fourier Transforms (FFTs) [7, 6]. Here, we take a dif-ferent approach by developing stochastic methods based on Finite Element Methodsfor fluctuating hydrodynamics and provide an alternative to Fast Fourier Transformsfor the generation of stochastic driving fields. Our approach allows for non-uniformspatially adaptive discretizations on non-periodic domains with geometries more nat-urally encountered in many applications.We develop Finite Element Methods with properties that facilitate the introduc-tion of stochastic driving fields and their efficient generation. We show our discretiza-tion approach provides operators that satisfy certain symmetry and commutationconditions that are important when subject to the incompressibility constraint forhow thermal fluctuations propagate throughout the discrete system. We formulatethe stochastic equations for our fluid-structure system subject to thermal fluctuationsin Section 1. We introduce for a given spatial discretization our general procedure forderiving compatible stochastic driving fields that model the thermal fluctuations in amanner consistent with statistical mechanics in Section 2. To obtain the stochasticdriving fields with the required spatial correlation structure, we develop stochasticiterative methods based on multigrid to generate the Gaussian random fields withcomputational complexity O ( N ) in Section 3. We present validation of our stochasticnumerical methods with respect to the hydrodynamic coupling and thermal fluctua-tions in Section 5. To demonstrate our approach in practice, we present simulationsof a few example systems in Section 6.Overall, our approach extends the range of problems that can be treated numer-ically with fluctuating hydrodynamic methods by allowing for arbitrary geometrieswith walls having no-slip boundary conditions and by allowing for spatially adaptiveresolution. Many of the central ideas used for our numerical approximation of thefluctuating hydrodynamic equations should also be applicable in the approximationof other parabolic Stochastic Partial Differential Equations (SPDEs). We expect ourstochastic numerical methods for fluctuating hydrodynamics to be useful in applica-tions where the domain geometry plays an important role.
1. Fluid-Structure Hydrodynamics and Fluid-Structure Interactions.
We describe the mechanics of fluid-structure interactions subject to thermal fluctua-tions using the Stochastic Eulerian Lagrangian Method (SELM) [6]. In the inertialregime this is given by momentum equations for the fluid coupled to momentumequations for the microstructures [6]. We consider here the regime in which the fluid-structure coupling is strong and the microstructures are mass density matched withthe fluid [6, 4]. This regime is closely related to the Stochastic Immersed BoundaryMethod [7, 4, 36, 14]. In this regime, we use the time-dependent Stokes equations forthe fluid coupled to an equation of motion for the microstructures ρ ∂ u ∂t = µ ∆ u − ∇ p + f s + f thm in Ω ∇ · u = 0 in Ω u | ∂ Ω = 0 . (1.1) TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS X are given by the following equationof motion and coupling condition that models the bidirectional coupling between thefluid and microstructures d X dt = Γ u (1.2) f s = Λ [ −∇ Φ( X )] . (1.3)The thermal fluctuations are taken into account by the Gaussian random field f thm which when decomposed into a mean and fluctuating part f thm = ¯ f thm + ˜ f thm hasthe form ¯ f thm = (cid:104) f thm (cid:105) = k B T ∇ X · Λ (1.4) (cid:104) ˜ f thm ( x )˜ f Tthm ( y ) (cid:105) = 2 µ ∆ C ( x − y ) (1.5) C ( x − y ) = k B T ρ − δ ( x − y ) . (1.6)These stochastic driving fields were derived for the mechanical system using the SELMframework in [6]. A notable difference with the original formulation of the Stochas-tic Immersed Boundary Method (SIBM) is the presence of the thermal drift term inequation 1.4 which arises from the more systematic treatment through stochastc av-eraging to obtain in this regime the stochastic fluid-structure equations which handlesthe generalised configuration-momentum coordinates used in such descriptions [6, 4].In the notation, the u is the fluid velocity field, X is the collective microstructureconfiguration, p the pressure, and µ the dynamic shear viscosity. The hydrodynamicequations 1.1 account for the microstructure interaction through f s . For short, welet f = f s + f thm . In the motion of the microstructures in response to the fluidflow is given by equation 1.2. Similar to the Immersed Boundary Method [36, 14],the operator Γ provides a model for how the microstructure locally responds to thefluid flow. The influence that the microstructure has on the nearby fluid is given byequation 1.3. The Λ operator models the neighborhood of surrounding fluid that isaffected by forces acting on the microstructure. These operators can be chosen quitegenerally provided they satisfy the adjoint condition [36, 14, 6] (cid:104) Λ V , u (cid:105) Ω = (cid:104) V , Γ u (cid:105) M (1.7)where u has the same form as the fluid velocity field and V the same form as themicrostructure velocity. The (cid:104) f , g (cid:105) Ω = (cid:82) Ω f ( x ) · g ( x ) dx denotes integration over thespatial domain of the fluid and (cid:104) F , G (cid:105) M = (cid:80) F j · G j denotes the dot-product over themicrostructure degrees of freedom. With these assumptions, the equations 1.1– 1.3describe the mechanics of the fluid-structure interactions subject to thermal fluctua-tions in the physical regime where the microstructure is strongly coupled to the fluidand mass density-matched with the fluid [6, 4].In the regime where the hydrodynamics relaxes rapidly relative to the time-scaleof microstructure motions, the fluid-structure dynamics can be further reduced [6, 4]to obtain d X dt = H [ −∇ Φ(( X ))] + k B T ∇ · H + h (1.8) (cid:104) hh T (cid:105) = 2 k B T H. (1.9)We refer to this as the overdamped Stokes regime which is similar to Brownian-Stokesian dynamics [10, 21]. The effective hydrodynamic coupling tensor of the mi-crostructure is given by H = Γ ℘ L ℘ T Λ with L = µ − ∆ − and ℘ = I − ∇ ∆ − ∇· is P. PLUNKETT, P.J. ATZBERGER the solenoidal projection operator imposing the hydrodynamic incompressibility [15].The thermal fluctuations are taken into account through the stochastic driving term h and the thermal drift-divergence term k B T ∇ · H . The above equations can be shownto have the Gibbs-Boltzmann distribution invariant with detailed-balance under thestochastic dynamics [6, 4].
2. Finite Element Semi-Discretization : Mixed-Method.
We developmixed finite element methods for semi-discretization of the fluid-structure systemwhere the fluid is governed by the time-dependent Stokes equations. We perform astochastic reduction of the semi-discretized equations in the limit where the hydrody-namics is assumed to relax rapidly on the time-scale of the motions of the microstruc-tures. This reduction provides a numerical discretization of the fluid-structure systemin the overdamped regime.In the weak formulation of the fluid equations, we consider the fluid velocity u in the space H (Ω) and the pressure p in the space L (Ω). The H (Ω) denotes thetriple product of the Sobolev space under the L -norm with weak derivatives up toorder one and zero trace on the domain boundary [11]. The weak formulation of thefluid equations is ρ (cid:104) ∂ t u , ϕ (cid:105) = − µa ( u , ϕ ) + b ( ϕ , p ) + (cid:104) f , ϕ (cid:105) , for all ϕ ∈ H (Ω) ,b ( u , ψ ) = 0 , for all ψ ∈ L (Ω) . (2.1)The (cid:104)· , ·(cid:105) denotes the L -inner product. The a : H (Ω) × H (Ω) → R and b : H (Ω) × L (Ω) → R are the continuous bilinear forms a ( u , v ) = (cid:104)∇ u , ∇ v (cid:105) = (cid:90) Ω ∇ u : ∇ v d x ,b ( u , q ) = (cid:104) q, ∇ · u (cid:105) = (cid:90) Ω q ∇ · u d x . We use the Ritz-Galerkin approximation of the variational problem corresponding torestriction to the finite dimensional subspaces V h ⊂ H (Ω) and P h ⊂ L (Ω). Ourspecific choice of spaces V h , P h will be discussed in Section 2.1. The exact solutionto the variational problem in equation 2.1 is approximated by u h ∈ V h and p h ∈ P h satisfying the finite-dimensional problem ρ (cid:104) ∂ t u h , ϕ i (cid:105) = − µa ( u h , ∇ ϕ i ) + b ( ϕ i , p h ) + (cid:104) f , ϕ i (cid:105) , for all ϕ i ∈ V h b ( u h , ψ i ) = 0 , for all ψ i ∈ P h . (2.2)This can be expressed in matrix form as ρM ˙ U = − µLU + GP + FDU = 0 . (2.3)We represent functions using a basis V h = span { ϕ i } and a basis P h = span { ψ i } with u h = (cid:88) i U i ϕ i and p h = (cid:88) i P i ψ i . The matrices have entries M ij = (cid:104) ϕ i , ϕ j (cid:105) , L ij = a ( ϕ i , ϕ j ) , (2.4) G ij = − b ( ϕ j , ψ i ) , D ij = b ( ϕ i , ψ j ) . (2.5) TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS F can be decomposed into a contribution from the microstructurecoupling and thermal fluctuations F = F str + F thm . The microstructure coupling term F str is given by F str i = (cid:104) f s , ϕ i (cid:105) = (cid:104) Λ [ −∇ Φ] , ϕ i (cid:105) which by linearity of Λ we can expressas F str i = B T X ( −∇ Φ) . In the weak formulation, the operator B T X now plays the analogous role as Λ inequation 1.3.In summary, for the full fluid-structure system, the finite element semi-discretizationcan be expressed as ρM ˙ U = − µL U + GP + B T X ( −∇ Φ) + F thm D U = 0˙ X = B X U ¯ F thm = (cid:104) F thm (cid:105) = k B T ∇ X · B T X (cid:28) ˜ F thm (cid:16) ˜ F thm (cid:17) T (cid:29) = G = µ (cid:0) LC + ( LC ) T (cid:1) C = k B T ρ − M − . (2.6)As in the continuum formulation, the thermal fluctuations are decomposed into themean and fluctuating parts F thm = ¯ F thm + ˜ F thm . The operator B X now plays theanalogous role as Γ and B T X the analogous role as Λ in equations 1.1– 1.3. Thisensures that the important adjoint condition between the force-spreading and velocity-averaging operators hold [36, 6, 7]. We caution that the stochastic term F thm was notderived by simply projecting the full stochastic field f thm . Following the approachin [6, 7], we took into account the particular properties of the discrete operators ofthe system to ensure fluctuations propagate throughout the discretized fluid-structuresystem in a manner consistent with statistical mechanics. We establish in more detailthe statistical mechanics of our finite element discretization in Section 2.2.For the overdamped regime, we follow an approach similar to [4, 6] to reduce thesemi-discretized equations 2.6 in the limit of rapid hydrodynamic relaxation to obtaina stochastic semi-discretization of the overdamped equations 1.8. This yields˙ X = H ( −∇ Φ) + k B T ∇ X · H + R (2.7) (cid:104) RR T (cid:105) = 2 k B T H. (2.8)The effective hydrodynamic coupling tensor is H = B X SB T X . The operator S repre-sents solving for the fluid velocity U in the following discretized Stokes equations A (cid:20) U P (cid:21) = (cid:20) B T X ( −∇ Φ)0 (cid:21) , where A = (cid:20) µL − GD (cid:21) . (2.9)This can be expressed more explicitly as S = [ I A − [ I T which gives U = SB T X ( −∇ Φ).For general finite elements a semi-discretization of the stochastic fluid-structuresystem in the regimes of inertial dynamics and overdamped dynamics is given by equa-tions 2.6 and 2.7. For use in practice, this semi-discretization requires a specific choiceof appropriate approximating function spaces and corresponding finite elements.
P. PLUNKETT, P.J. ATZBERGER V h and P h . The Stokes problem given by equation 2.9 for our mixed finite elementdiscretization is well-known to be a saddle problem. To ensure numerical stability anda well-posed variational problem requires a careful choice for the two function spaces V h and P h [11]. The Stokes problem can be split into two sub-problems. In the first,one attempts to find the fluid velocity u h within the subspace Z ⊂ V h of vector fieldssatisfying exactly the incompressibility constraint, Z = { v ∈ V h | b ( v , q ) = 0 , ∀ q ∈ P h } .This is done by considering the variational problem [11] a ( u h , ϕ ) = −(cid:104) f , ϕ (cid:105) , ∀ ϕ ∈ Z ⊂ V h . (2.10)In the second, one attempts to solve for the pressure p in P h by considering thevariation problem [11] b ( ϕ , p ) = − a ( u h , ϕ ) − (cid:104) f , ϕ (cid:105) , ∀ ϕ ∈ V h . (2.11)Provided both of these two sub-problems can be solved, we have a consistent solutionto the Stokes problem. A central issue is that for mixed problems the bilinear form b is often not coercive. This prevents a direct application of the Lax-Milgram Theoremto ensure well-posedness [11]. As an alternative, Babuˇska and Brezzi [8, 27, 12] foundfor the two bilinear forms a , b on the function spaces V h and P h a set of conditionsthat are weaker than coercitivity but still sufficient to ensure well-posedness. The firstcondition concerns the bilinear form a and amounts to a form of coercitivity but nowweaker only requiring this property when restricting to the subspace Z as suggestedby the sub-problem in equation 2.10 α (cid:107) v (cid:107) ≤ a ( v , v ) , ∀ v ∈ Z ⊂ V h . (2.12)The second condition requires the bilinear form b satisfy for the two function spacesinf q h ∈P h sup v h ∈V h b ( v h , q h ) (cid:107) v h (cid:107) V h (cid:107) q h (cid:107) P h ≥ β > . (2.13)This is referred to as the Babuˇska-Brezzi condition or the inf-sup condition . Babuˇskaand Brezzi showed for mixed methods these conditions are sufficient to ensure well-posedness of the variational problems and desirable numerical properties [8, 27, 12].The Babuˇska-Brezzi condition provides important guidance on which finite el-ements should be used for the fluid and pressure in the Stokes problem. For theStokes problem, many of the most obvious choices of function spaces do not workvery well in practice and in fact violate the Babuˇska-Brezzi condition given by equa-tion 2.13 [3, 11]. For instance using V h = ( P k ) , P h = P k where P k is the spaceof piecewise-polynomials of degree k does not satisfy the Babuˇska-Brezzi conditionfor any choice of k >
1. Even the choice V h = ( P ) , P h = P does not satisfy theconditions to provide a stable method.To satisfy the Babuˇska-Brezzi conditions, we use finite elements for the velocityfield that are enriched by an additional bubble mode V h = ( P -bubble) and forthe pressure P h = P , as introduced in [3]. The P -bubble consists of the usualpiecewise linear functions but enriched with a quartic “bubble element” located atthe barycenter of each tetrahedral element but vanishing at the faces, see Figure 2.1.The enrichment by bubble-elements provides enough stabilisation that the Babuˇska-Brezzi conditions are satisfied and yield for the Stokes equations a convergent method TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS Fig. 2.1 . P -MINI Elements [3]. To obtain a discretization of the Stokes equations satisfyingthe Babuˇska-Brezzi condition we use a combination of piecewise linear elements (P1-elements) forthe pressure and quartic Bubble-elements for the velocity. The mesh degrees of freedom (DOF)consist of the usual nodal variables for P1-elements (labeled in red) and nodal variables at the centerof each element to determine the bubble mode (labeled in blue). While the bubble-elements contributea large number of degrees of freedom, they have the convenient property that their support is onlywithin the interior of an element and are decoupled from one another [3]. which has order h [11]. Given that the bubble-elements do not overlap since they arecontained within the interoir of each element, the overhead associated with the bubbleenrichment is marginal since they contribute in a manner to the overall linear systemthat is decoupled. This also has desirable properties in the stochastic setting whengenerating fluctuations as we shall discuss. The combination ( P -bubble) / P has aminimal footprint. We shall refer throughout to this pair as the P -MINI Elements,as in [3]. We show that ourstochastic semi-discretization provides an approximation yielding fluctuations consis-tent with statistical mechanics. We show that the Gibbs-Boltzmann distribution isinvariant under the stochastic dynamics. The Gibbs-Boltzmann distribution is givenby Ψ GB ( z ) = 1 Z exp [ − E ( z ) /k B T ] . (2.14)The z = ( u , X ) is the state of the system, E is the energy, k B is Boltzmann’s con-stant, T is the system temperature, and Z is a normalization constant for the distri-bution [38]. The energy of the fluid-structure system is given by the kinetic energy ofthe fluid and the potential energy of the microstructures E [ z ] = ρ (cid:90) | u | dx + Φ( X ) . (2.15)For the discrete system we consider the energy E [ z ] = ρ U T M U + Φ( X ) . (2.16)The M is the mass matrix defined in equation 2.4. For this energy, we have atequilibrium that the fluid velocity U has fluctuations with mean zero and covariance C = k B T ρ − M − . (2.17) P. PLUNKETT, P.J. ATZBERGER
For the stochastic fluid-structure equations 2.6, the Fokker-Planck equation for Ψ GB is given by ∂ Ψ GB ∂t = −∇ · J (2.18) J = (cid:20) ρ − M − (cid:0) µL + B T X + k B T ∇ X · B T X (cid:1) B X (cid:21) Ψ GB −
12 ( ∇ · G )Ψ GB − G∇ Ψ GB . The G denotes the covariance operator for the stochastic driving fields given by equa-tion 2.6. The Gibbs-Boltzmann distribution is invariant provided that ∇ · J = A + A + ∇ · A + ∇ · A = 0 (2.19)where A = (cid:2) ρ − M − (cid:0) B T X + k B T ∇ X · B T X (cid:1) · ∇ U E + B X · ∇ X E (cid:3) ( − k B T ) − Ψ GB A = (cid:2) ∇ U · (cid:0) ρ − M − (cid:0) B T X + k B T ∇ X · B T X (cid:1)(cid:1) + ∇ X · B X (cid:3) Ψ GB A = −
12 ( ∇ · G )Ψ GB A = (cid:20) µL + [ G UU ∇ U E + G UX ∇ X E ] (2 k B T ) − [ G XU ∇ U E + G XX ∇ X E ] (2 k B T ) − (cid:21) Ψ GB . To simplify the notation we have suppressed explicitly denoting the functions on whichthe operators act, which can be inferred from equation 2.6. We have also suppressedthe incompressibility constraint since this can be handled readily by considering thedynamics decomposed into a component on the space of solenoidal fields and itsorthogonal complement as was done in [4]. To compare terms, we use the gradientsobtained from the energy in equation 2.16 given by ∇ U E = ρM U (2.20) ∇ X E = ∇ X Φ . (2.21)By direct substitution of the gradients into A given by equation 2.20– 2.21, we findafter cancellation A = − ( ρ − M − ( ∇ X · B T X ) · ∇ U E )Ψ GB = − ( U T ( ∇ X · B T X ))Ψ GB .Since B T X only depends on X , we have A = ( ∇ X · B X U )Ψ GB . This gives that A + A = 0. We remark that this is a direct consequence of requiring the couplingoperators to be linear and adjoints Γ = Λ † in the sense of equation 1.7 where Γ = B X [6, 4]. The term A accounts for probability fluxes driven by state dependentchanges in the covariance of the stochastic driving fields. In this case the covariance G does not depend on the system state, so the divergence gives A = 0. The term A arises from the interplay between dissipation and fluctuations in the dynamics.By looking at the dependence of the differentiated expressions most of the termsare seen to be zero. By our choice of the covariance G u , u given by equation 2.6which was motivated by the form of the discrete energy in equation 2.16 and thetarget equilibrium covariance for fluctuations given in equation 2.17, we have thatthe non-zero terms balance to yield A = 0. The choice of G and its relation to C can be thought of as imposing a discrete fluctuation-dissipation balance for oursemi-discretization [7, 6]. These results establish that ∇ · J = 0 and that the Gibbs-Boltzmann distribution is invariant under our semi-discretized stochastic dynamics.In the overdamped limit, the fluid-structure system is governed only by the mi-crostructure potential energy E [ X ] = Φ( X ) . (2.22) TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS GB are given by ∂ Ψ GB ∂t = −∇ · J (2.23) J = ( H ( −∇ X Φ( X )) + k B T ∇ · H ) Ψ GB −
12 ( ∇ · k B T H )Ψ GB −
12 (2 k B T H ) ∇ X Ψ GB . From the form of Ψ GB , we have − (2 k B T H ) ∇ X Ψ GB = H ( ∇ X Φ( X ))Ψ GB which givesthat J = 0. This shows that the Gibbs-Boltzmann distribution is invariant withdetailed-balance under the stochastic dynamics of our semi-discretization given inequation 2.7.These results establish that our semi-discretizations both in the inertial and over-damped regimes yield stochastic dynamics consistent with statistical mechanics. Animportant issue in practice is to generate efficiently the stochastic driving fields F thm and h with the required covariances.
3. Generation of the Stochastic Driving Fields.
The thermal fluctuationsin the overdamped regime require generating a Gaussian random field R with meanzero and covariance (cid:104) RR T (cid:105) = 2 k B T H = 2 k B T B X SB T X . (3.1)In general, this is difficult since S is a dense matrix and calculating the square rootusing methods such as Cholesky yields a dense factor making field generation com-putationally expensive. The issue is further complicated by the fact that the actionof S involves taking the inverse of A which is not a sign definite matrix. We take analternative approach for the efficient generation of the stochastic field R by factoringthe Stokes operator S to reduce the problem to generating variates with a covariancerelated only to the discrete Laplacian stiffness matrix L . A useful identity for the discrete Stokes operator is S ( µL ) S = S. (3.2)This follows since for an arbitrary choice of F in equation 2.7 we have the solution U = SF so that S ( µL ) SF = SF + SGP from substitution of U back into theStokes equation 2.7. The SG = 0 since for an arbitrary P if we set F = GP thenthe solution is clearly U = SF = SGP = 0 by uniqueness the Stokes problem for U . This identity is useful since it reduces generating the stochastic driving field R with covariance 2 k B T H to that of generating variates ξ with covariance C = − ( µL ). More specifically, we can generate the stochastic driving field by using R = √ k B T B X S ξ which yields the required statistics through (cid:104) RR T (cid:105) = 2 k B T B X S (cid:104) ξξ T (cid:105) S T B T X = 2 k B T B X S ( µL ) SB T X = 2 k B T B X SB T X = 2 k B T H.
To generate efficiently the variates ξ with covariance − L , we develop stochastic iter-ative methods.We remark that for this factorization there are some further advantages whenusing the P -MINI elements. For the Laplacian stiffness matrix L , the collectionof entries corresponding to the bubble elements is diagonal. This allows for trivial0 P. PLUNKETT, P.J. ATZBERGER generation of variates at the bubble elements by using independent Gaussians andweighting by the square root of the diagonal entry of the stiffness matrix. This meansthat stochastic field generation can be reduced to the problem of generating variateswith covariance − L restricted to the non-bubble nodal degrees of freedom which arefar fewer in number than the bubble nodes, see Figure 2.1. This provides in thestochastic setting a particular advantage of the P -MINI elements over higher ordermethods such as the Taylor-Hood elements [3, 11]. We develop stochastic iterative methodsto obtain Gaussian random variates with a target mean and covariance structure. Ageneral one-step linear iterative method can be developed into a Gibbs sampler byintroducing a stochastic term η n into the iteration step as ξ n +1 = M ξ n + N b + η n . (3.3)We assume throughout that η n is a Gaussian variate with mean zero and covariance G = (cid:104) η n η n,T (cid:105) . The stochastic iterations can be expressed in terms of the probabilitydensity ρ n ( ξ ) and transitions at iteration n as ρ n +1 ( ξ ) = (cid:90) π ( ξ , w ) ρ n ( w ) d w (3.4) π ( ξ , w ) = 1 √ π det G exp (cid:104) ( ξ − M w − N b ) T G − ( ξ − M w − N b ) (cid:105) . (3.5)This yields a set of Gaussians ξ n having mean value µ n and covariance C n sat-isfying µ n +1 = M µ n + N b (3.6) C n +1 = M C n M T + N b µ n,T M T + M µ n b T N T + N bb T N T + G. (3.7)In the case that the target mean µ = 0, we can choose b = 0 and this simplifiesto C n +1 = M C n M T + G. (3.8)The autocorrelation Φ m = (cid:104) ξ n ( ξ n + m ) T (cid:105) for how the random variates ξ n decorrelateover m iterations satisfies Φ m +1 = M Φ m . (3.9)This relation has the consequence that in the stochastic setting the decay in corre-lation has the same behavior as the decay in error in the deterministic setting wheniteratively solving the problem L x = b . To obtain random variates with a targetcovariance C , we have from equation 3.8 that the covariance of η must satisfy G = C − M CM T = − AC − CA T + ACA T (3.10)where A is defined by M = I − A . For generating general Gaussian random variates ξ ,the stochastic iterative method 3.3 provides a useful approach provided the iterativemethod converges efficiently and the random variates η have a covariance G that iseasier to generate than C . For traditional iterative methods, such as SOR, Gauss-Seidel, and Jacobi iterations, stochastic iterative counterparts have been introducedin [24, 23, 1, 49]. The convergence of such stochastic iterative methods can be furtherimproved by using preconditioning strategies such as multigrid [24, 23]. For stochasticfield generation in our semi-discretized fluctuating hydrodynamic equations, we showhow the multigrid preconditioner can be adopted to sample the stochastic drivingfields. We base our sampler on Gauss-Seidel iterations. TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS The Gauss-Seidel iterative method for solving the system A v = b in the deterministic setting isgiven by v n +1 = ( D − L ) − U v n + ( D − L ) − b . (3.11)We use the splitting A = D − L − U . The D denotes the diagonal elements of thematrix A , L the lower triangular elements, and U the upper triangular elements.This corresponds to the general iterative method of equation 3.3 with the choice M = ( D − L ) − U and N = ( D − L ) − . To obtain Gaussian random variates withmean µ = 0 and covariance C = A − , we must generate for each iteration a stochasticvariate η n having covariance G satisfying equation 3.10. For Gauss-Seidel iterations,this requires covariance G = ( D − L − U ) − − ( D − L ) − U ( D − L − U ) − L ( D − L ) − T . (3.12)It has been shown in [24] that this can be expressed as G = ( D − L ) − (cid:2) ( D − L − U )( D − L − U ) − ( D − U )+ U ( D − L − U ) − ( D − L − U ) (cid:3) ( D − U ) − (3.13)= ( D − L ) − D ( D − U ) − . This gives the factorization G = QQ T with Q = ( D − L ) − D / . (3.14)We have used that A is symmetric positive definite with U = L T and that the trans-pose of the inverse is the inverse of the transpose. This provides a way to generate therandom variates through η = Q ζ where ζ is a standard multi-variate Gaussian withmean zero and unit covariance. This is an efficient procedure, especially for sparsematrices A , since the action of ( D − L ) − can be obtained through back-substitutionin the same manner as in deterministic Gauss-Seidel iterations [13]. This provides ourstochastic Gauss-Seidel iterative method for the generation of random variates. In the deterministic setting multigrid is often usedto improve the convergence of iterative methods [13, 34]. Multigrid makes use ofdifferent levels of refinement (cid:96) to perform iterations. In the multigrid iterations threefundamental operators are utilized. The first is the smoother operator S (cid:96) which isused to iteratively approximate the linear system of equations on a given level (cid:96) .For transmission of information between different levels of resolution, a restrictionoperator I (cid:96) − (cid:96) and prolongation operator I (cid:96)(cid:96) − are defined. The restriction operator I (cid:96) − (cid:96) maps data from a more refined level (cid:96) to data on a less refined level (cid:96) −
1. Theprolongation operator I (cid:96)(cid:96) − maps data from a less refined level (cid:96) − (cid:96) , see Figure 3.1. Multigrid provides improvements in the convergencein the deterministic setting by the different rates that the smoother relaxes modeson different levels. This same feature that relaxes errors translates to the stochasticsetting by improving the rate at which the generated random variates decorrelate inthe number of iterations.To ensure a viable stochastic multigrid method it is useful to have a few additionalproperties that are not strictly required in the deterministic setting. We take the linear2 P. PLUNKETT, P.J. ATZBERGER
Fig. 3.1 . Multigrid uses subproblems on different levels of refinement to iteratively solve thesystem of equations. Three operators are utilized: (i) smoother operator, (ii) restriction operator,and (iii) prolongation operator. The smoother operator serves on each level to iteratively relaxvalues toward the solution. The restriction operator and prolongation operator serve to transfer databetween levels (left). Multigrid iterations are performed by a protocol combining these operations(right) [13]. equations A ( (cid:96) ) v = b ( (cid:96) ) at refinement level (cid:96) to be related to those at the most refinedlevel (cid:96) ∗ by A ( (cid:96) ) = I (cid:96)(cid:96) ∗ A I (cid:96) ∗ (cid:96) , b ( (cid:96) ) = I (cid:96)(cid:96) ∗ b . (3.15)An important property when adopting the multigrid method to the stochastic settingis to ensure a consistent variational principle for the linear equations at different re-finement levels. The variational principle at the most refined level is that the solutionis a minimizer of the energy E ( v ) = 12 v T A v − v T b . (3.16)Consistency requires that the energy defined by E (cid:96) ( v ( (cid:96) ) ) := E ( I (cid:96) ∗ (cid:96) v ( (cid:96) ) ) provide thevariational principle at level (cid:96) . This is satisfied when the prolongation and restrictionoperators are adjoints I (cid:96)(cid:96) ∗ = (cid:0) I (cid:96) ∗ (cid:96) (cid:1) T , which yields E (cid:96) ( w ) = 12 w T (cid:16) I (cid:96) ∗ (cid:96) (cid:17) T A (cid:16) I (cid:96) ∗ (cid:96) (cid:17) w − w T (cid:16) I (cid:96) ∗ (cid:96) (cid:17) T b (3.17)= 12 w T A ( (cid:96) ) w − w T b ( (cid:96) ) . (3.18)The variational property has some important consequences. For the target Gaussiandistribution ρ ( v ), the stochastic smoother samples at level (cid:96) the marginal probabilitydistribution ρ ( (cid:96) ) ( w ) = (cid:82) { v = I (cid:96) ∗ (cid:96) w } ρ ( v ) d v , which is given by ρ ( (cid:96) ) ( w ) = 1 √ π det A ( (cid:96) ) exp (cid:20) − w T A ( (cid:96) ) w + w T b ( (cid:96) ) (cid:21) . (3.19)This variational property can be shown to be sufficient to ensure the probabilitydistribution of the target multi-variate Gaussian is the invariant distribution of thestochastic multigrid iterations [24]. TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS geometric multigrid in which the restriction and prolongation op-erators are constructed by considering geometric correspondences between differentlevels of the discretization such as spatial averages or interpolations [13]. The secondis algebraic multigrid in which the restriction and prolongation operators are con-structed by considering the algebriac structure of the linear system such as groupingclusters of entries of the matrix A [20, 43, 45, 44, 46]. We use an algebraic multigrid method with smoothed aggregation following an approach similar to [20].4 P. PLUNKETT, P.J. ATZBERGER
4. Algorithm : Summary of SELM based on FEM Stokes.Input : Polyhedral domain Ω, initial body configuration X , body potential Φ, ter-minal time t end , timestep size ∆ t , Peskin δ footprint h , kinematic viscosity µ , tem-perature T . Output : Body configuration X t at each t = ∆ t, t, t, . . .
1. Construct the adapted tetrahedralization: {T , V } = AdaptedMesh (Ω , X , h ).2. Construct the linear operators: { L, D } = FluidOperators ( {T , V } ).3. Set A = (cid:20) µL D T D (cid:21) .4. Set t, t remesh = 0.5. Begin computing the body trajectory. While t < t end , do steps (a)-(g):(a) Construct the velocity interpolator/force spreading operator: B = VelocityInterpolator ( {T , V } , X t ).(b) Generate the N ( , ( µL ) − ) sample: ξ = GaussianSample ( µL ).(c) Generate the body forces: f = B T ( −∇ Φ( X t )).(d) Perform the fluid solve: U = FluidSolve ( A , f ∆ t + √ k B T ∆ tµL ξ )(e) Interpolate fluid velocity to particles: X t +∆ t = X t + B U (f) Test for remeshing. If | X it +∆ t − X it remesh | > h for any i , do steps i-iv:i. Reconstruct the (adapted) tetrahedralization: {T , V } = AdaptedMesh (Ω , X t +∆ t , h ).ii. Reconstruct the linear operators: { L, D } = FluidOperators ( {T , V } ).iii. Set A = (cid:20) µL D T D (cid:21) .iv. Set t remesh = t + ∆ t .(g) Update the current time: set t = t + ∆ t . TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS
5. Validation Studies.
We perform several studies to establish the validity ofthe computational methods. We first investigate the empirical covariance structureobtained from our iterative Gibbs samplers to make comparisons with the target co-variance structure. We then consider the decorrelation rates exhibited by the directstochastic Gauss-Seidel iterations in comparison to the stochastic multigrid itera-tions. We next explore the effective hydrodynamic coupling tensors corresponding tothe case when the fluid and structures are coupled based on the
Immersed BoundaryMethod [36, 14]. We then investigate the equilibrium fluctuations of the fluid-structuresystem by computing the statistics of diffusing particles subject to a harmonic po-tential. To validate our overall method we compare these results with the predictedGibbs-Boltzmann distribution.
A key component of our computational methodsfor fluctuating hydrodynamics is to generate efficiently the stochastic fields drivingfields. As discuss in Section 3, this requires methods to generate Gaussian variateswith target covariance − L − . We consider the fluctuating hydrodynamics when con-fined within a spherical domain with Dirichlet no-slip boundary conditions. For bothdirect stochastic Gauss-Seidel iterations and stochastic multigrid iterations we gen-erated 10,000 samples, shown in Figure 5.1. The covariance matrix is shown for allof the finite element degrees of freedom of the discrete system including the bubblemodes. The magnitude of the covariance is shown on a logarithmic scale. We findthat the empirical covariance structure from both of our stochastic iterative methodsis in very good agreement with the target covariance structure with an average errorof less than 4%. An important consideration is how efficiently the stochastic iterativemethods can sample nearly independent Gaussian random variates.6 P. PLUNKETT, P.J. ATZBERGER
Exact Covariance
Fig. 5.1 . Covariance Structure. The matrix of covariance for all finite element degrees offreedom are shown on a logarithmic scale for 10,000 samples from our stochastic Gauss-Seideliterations (GS) and stochastic multigrid iterations (AMG) (top row). These are compared to thetarget covariance − L − which is the negative inverse Laplacian. The diagonal entries show theself-correlation while the off-diagonal entries show correlations between distinct degrees of freedom.The lower right entries of the correlation matrix is shown in more detail on the bottom. Bothstochastic iterative methods yield results with good agreement with the target covariance structurewith an average error of 4%. TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS The stochastic driving fields appearing in the fluc-tuating hydrodynamic equations 1.1 and 1.8 are treated as uncorrelated in time. Instochastic numerical methods this requires on successive time steps the generationof independent Gaussian random variates to account for the thermal fluctuations.When using stochastic iterative methods to sample random variates there are alwayscorrelations between generated variates but these diminish over successive iterations.This means the quality of the random variates for fluctuating hydrodynamic simu-lations requires taking a sufficient number of iterations. The efficiency depends onhow rapidly these correlations decay on successive iterations. We investigate thisby computing empirically a random variate its autocorrelation statistics as given inequation 3.9. We consider the fluid within a spherical domain with no-slip Dirich-let boundary conditions with a finite element discretization having 460,904 degreesof freedom. The autocorrelation is considered for direct stochastic Gauss-Seidel it-erations and stochastic multigrid iterations counting both the number of iterationsand the number of Gauss-Seidel visitations to individual finite element nodes. Thevisitations count more closely reflects the computational expense in acheiving a givenlevel of decorrelation with a given method. We find for the stochastic Gauss-Seideliterations that many more iterations are required to acheive decorrelated variates, seeFigure 5.2. In fact, to achieve a correlation between the random variates of less than1% requires on the order of 100 Gauss-Seidel iterations. In contrast, the stochasticmultigrid method achieves a correlation less than 1% in less than 10 iterations. To fur-ther characterise the performance of the stochastic iterative methods, we consideredhow the random variates decorrelate as the Gauss-Seidel smoother visits individualfinite element nodes. We find that standard Gauss-Seidel iterations achieve correla-tions between the random variates less than 1% after 10 vistiations. In contrast,the stochastic multigrid method makes better use of Gauss-Seidel visitations throughcoarsening of the system and acheives correlations between the random variates lessthan 1% after 10 visitations (two orders of magnitude less). These results showthe significant improvement provided by using stochastic multigrid to sample randomvariates. The differences between the stochastic sampling methods are expected tobecome even more significant when increasing the system size.8 P. PLUNKETT, P.J. ATZBERGER C o rr( X , X n ) Iterations nCorrelation over IterationsSlope: -3.1e-03 Slope: -1.4e-01Gauss-SeidelAMG0.0010.010.1 0 200 400 0 5 10 15 C o rr( X , X n ) Visitations nCorrelation over VisitationsSlope: -2.8e-09 Slope: -5.7e-07Gauss-SeidelAMG0.0010.010.1 0 2e8 4e8 0 2e6 4e6
Fig. 5.2 . Decorrelation of random variates. Sampling was perform for the stochastic drivingfields on a spherical domain with Dirichlet boundary conditions for the target covariance of thediscrete inverse Laplacian − L − . We used both the stochastic Gauss-Seidel and stochastic multigridmethods. The decorrelation in the number of iterations is shown on the left. The decorrelation inthe number of Gauss-Seidel visitations to the individual nodes of the finite element are shown onthe right. We find that the stochastic multigrid method greatly improves performance. To achieve acorrelation of less than between random variates, direct stochastic Gauss-Seidel requires on theorder of iterations. In contrast, stochastic multigrid requires less than iterations. For thenumber of visitations, the stochastic Gauss-Seidel requires on the order of steps. In contrast, thestochastic multigrid iterations only require steps. The rate of exponential decay for stochasticGauss-Seidel is estimated to be − . × − and the rate for the stochastic multigrid iterations − . × − . TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS Many approaches can be used to couple themicrostructure to the fluid including the Force Coupling Method [33], DiscontinousGalerkin Finite Element Methods [22, 31], or the Immersed Boundary Method [36, 14].The simplest of these is to couple the microstructure and fluid using the ImmersedBoundary Method [36, 14]. This corresponds to the specific choice of coupling oper-ators Λ X ( F ) = F ( x ) δ a ( x − X ) (5.1)Γ u = (cid:90) Ω u ( x ) δ a ( x − X ) d x . (5.2)The δ a is the Peskin- δ function which is non-zero over the distance a [36]. To adoptthis approach in the current finite element setting the operators are discretized byapproximating the integral for Λ using the finite element basis to obtain the numericalcoupling operator B X ( · ). The other coupling operator Γ is approximated by usingthe adjoint condition in equation 2.6 which yields the numerical operator B T X , seeAppendix A. H (r) / H ( ) RPY00.20.40.60.81 0 1 2 3 4 5 6 7 8 r/aHydrodynamic Coupling Tensor Hydrodynamics
Close agreement with RPY
Oseen R PY Fig. 5.3 . Hydrodynamic Coupling. We find the Immersed Boundary Coupling resultsin hydrodynamic coupling that closely resembles the Rotne-Prager-Yamakowa Tensor, see rightpanel [5, 40, 50]. Both tensors in the near-field provide a regularization of the hydrodynamic cou-pling. In the far-field both tensors give the same results as the Oseen tensor. The hydrodynamic flowgenerated under the Immersed Boundary Coupling in response to forces acting on the two particlesis shown on the right.
For this approximate approach to the physics of fluid-structure interaction, weinvestigate the hydrodynamic response when forces are applied to a pair of particlesand their velocities. This corresponds to averaging the fluid using Γ and solving thesteady-state Stokes equations when forces are spread to the fluid using Λ. This re-sponse is characterized by the effective hydrodynamic coupling tensor H SELM . Wefind that the effective hydrodynamic coupling tensor H SELM very closely agrees withthe Rotne-Prager-Yamakawa tensor H RP Y [40, 50], see Figure 5.3. This is in agree-ment with our prior work in the context of finite volume methods reported in [5].These results establish the approximate way in which fluid-structure interactions are0
P. PLUNKETT, P.J. ATZBERGER handled give the expected results in the far-field and in the near-field exhibit a well-characterized regularization.
TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS Fig. 5.4 . Spherical domain discretized by a tetrahedralization with non-uniform spatial resolu-tion. The full spherical domain is shown on the left, a cross-section revealing the interior elementsis shown in the middle, a wireframe representation of the elements is shown on the right. This spe-cific discretization was used in the validation studies of a Brownian particle diffusing in a harmonicpotential. Given how the harmonic potential is expected to confine particles toward the middle, agreater level of spatial refinement is used near the domain center.
To investigate thefluctuations of the discretized fluid-structure system, we consider the simulation ofdiffusing particles that are subjected to forces from a harmonic potential that tethereach particle to the origin. From equilibrium statistical mechanics the probabilitydistribution of the particle positions is predicted to be the Gibbs-Boltzmann distri-bution Ψ( X ) = 1 Z exp (cid:18) − E [ X ] k B T (cid:19) , E [ X ] = K X . (5.3)The E denotes the energy of a configuration, K = 7 . × − ag/ns the Hookeanspring stiffness, T = 300 K the temperature, and Z denotes the partition function andnormalizes the probability density [38]. The physical units we use are attograms (ag),nanoseconds (ns), nanometers (nm), and Kelvin (K). The particles are taken to diffusewith an effective radius of 10nm on a spherical domain of radius 1000nm, see figure5.4. Samples were collected from 18 separate trajectories each of which has at least1000 timesteps each. The first 10 timesteps were discarded in each case. We find thatthe computational simulations show very good agreement with the predictions of sta-tistical mechanics, see Figure 5.5. These results show that the stochastic driving fieldswe have introduced and the stochastic iterative sampling methods provide appropriatefluctuations in the discrete system to account for the thermal fluctuations.
6. Application.6.1. Diffusion of Particles within a Microfluidic Device Geometry.
Weshow how our methods can be used to capture effects on particle diffusion/mobilitywhen confined within a microfluidic channel having a complex geometry. We showhow our methods can be used for hydrodynamic flows through post-like obstacles tocapture the particle-wall hydrodynamic interactions and related correlation effects inthe thermal fluctuations of diffusing particles. The specific device geometry and ourfinite element discretization using the P -MINI elements of Section 2.1 is shown inFigure 6.1.We consider the diffusion of particles in a regime of moderate flow through thedevice. In this case, the main role of diffusion is in the directions lateral to the flowserving to change over time the effective streamline followed by a particle. Given thegeometry and small dimensions of the device relative to the particle size, significant2 P. PLUNKETT, P.J. ATZBERGER
Fig. 5.5 . Probability density of a particle’s location when diffusing in a harmonic potential.The probability distribution is shown for the particle location when expressed in spherical coordinates ( r, ϕ, θ ) (top) and Cartesian coordinates ( x, y, z ) (bottom). The simulations results (red) show goodagreement with the Gibbs-Boltzmann distribution (green). Fig. 6.1 . Diffusion of particles in a domain with geometry motivated by microfluidic device(left). The domain consists of a long channel with cylindrical posts obstructing flow near the center(right). The channel side walls and cylindrical posts are treated as having no-slip boundary condi-tions. The long-axis of the channel has periodic boundary conditions at the ends. The domain isdiscretized using a tetrahedralization and the P -MINI finite elements. hydrodynamic coupling can occur between the particle and walls. The diffusivity maychange significantly depending on the particular particle location within the channel.To investigate this effect, we performed simulations by starting particles at severallocations within the microfluidic channel. These included (i) placing particles near thechannel wall, (ii) placing particles near the cylindrical posts, and (iii) placing particlesnear the channel center. Simulation trajectories of particles in these locations is shownin Figure 6.2. The trajectories starting near the wall and cylindrical posts seem toexhibit smaller fluctuations than when started near the channel center away from suchboundaries. To investigate this further, we consider in the absence of flow the particle TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS D = 2 k B T M. (6.1)The particle mobility M gives the velocity response V of a particle to an applied force F V = M F . (6.2) Fig. 6.2 . Diffusivity of particles started at select locations along a cross-section of the chan-nel. Particles are subject to hydrodynamic flow through the device channel from left to right. Thestreamlines show the influence of the channel walls / cylindrical posts and traceout how passive par-ticles would be transported by the flow. For particles subject to thermal fluctuations the diffusivitycan causes particles to change streamlines. For instance, the bottom magenta trajectory starts ona streamline going below the post but the diffusivity moves the particle to a streamline carrying iton the other side of the post. In another example, the trajectory in black (bottom panel) shows aparticle which migrates very little when it comes into close proximity to the leftmost channel postwhich reduces both the fluid flow and particle diffusivity.
We consider how the mobility changes based on the particle location along a lineof positions passing across the channel near to the cylindrical posts, shown in Figure6.3. As the particles become close to one of the cylindrical posts we expect a dropin the mobility. This was found to be significant over the positions considered witha drop in mobility of around 30% relative to locations away from the post. Whilethe most significant drop in mobility occurs in the direction for the particle to movetoward the cylindrical post, we find the other directions also exhibit a non-negligibledrop. We also found cross-terms in the mobility which indicate that the cylindrical4
P. PLUNKETT, P.J. ATZBERGER post geometry induces correlations between the particle motions in different direc-tions, shown in Figure 6.3 on the right. In particular, we find that the componentsin the direction from the particle to the post and in the lateral direction give non-negligible correlations. This is in agreement with the streamlines that are found whenthere is fluid flow through the device, see inset in Figure 6.3. These results indicatethat the diffusivity of particles can exhibit significantly different qualitative behaviorsdepending on the particle location within the device. The results show the overallpotential of our computational methods to capture both the thermal fluctuations andhydrodynamic coupling in a consistent manner within non-trivial device geometries.While we demonstrated the methods for particles, our approach can also be used tosimulate the elastic responses to flow of more complex spatially extended microstruc-tures such as polymeric filaments or flexible membranes. Overall, we expect ourcomputational methods to be applicable quite broadly to the simulation of systemsinvolving hydrodynamically coupled microstructures confined within domains withcomplex geometries and subject to thermal fluctuations. -0.200.20.40.60.81-1 -0.8 -0.4 0 0.4 0.8 1 M / M x/LMobility (cross-coupling) -0.1-0.0500.050.1 -0.8 -0.4 0 0.4 0.8 H / H x/L M / M x/LMobilityM M M M M M cross-section Fig. 6.3 . Mobility dependence on the particle location within the channel. The mobility is shownfor a line intersecting the channel near the cylindrical posts (inset). The mobility significantly dropsto around near the channel walls then recovers but drops again by about when encounteringthe cylindrical posts. The region near the cylindrical posts is also found to induce some cross-couplingin the mobility components, particularly M . These results show the potential for our computationalmethods to capture in the particle motions the effects of walls and related confinement effects. Thefull three dimensional geometry is shown in Figure 6.1 TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS
7. Conclusions.
We have developed finite element discretizations for fluctuat-ing hydrodynamic simulations of fluid-structure interactions subject to thermal fluc-tuations. We introduced a discrete fluctuation-dissipation principle to obtain semi-discretizations of the stochastic equations. Our approach accounts explicitly for therole of discretization errors on the thermal fluctuations so that the Gibbs-Boltzmanndistribution is invariant exactly under the stochastic dynamics of the semi-discretizedsystem. To compute efficiently these required stochastic driving fields, we developed aGibbs sampler based on stochastic iterations of multigrid. The computational meth-ods are expected to extend significantly the applications that can be treated withfluctuating hydrodynamic approaches. The computational methods allow for spa-tially adaptive resolution of the mechanics and the treatment of complex geometriesoften relevant in application.
8. Acknowledgments.
The authors would like to acknowledge support fromDOE ASCR CM4. The author P.J.A. acknowledges support from research grant NSFCAREER DMS-0956210 and W.M. Keck Foundation. The Trilinos packagess MLand Epetra were used for the application simulations. The authors thank AlexanderRoma, Boyce Griffith, Mike Parks, and Micheal Minion for helpful suggestions.
REFERENCES[1]
Stephen L. Adler , Over-relaxation method for the monte carlo evaluation of the partitionfunction for multiquadratic actions , Phys. Rev. D, 23 (1981), pp. 2901–.[2]
Ann S. Almgren, John B. Bell, Phillip Colella, Louis H. Howell, and Michael L. Wel-come , A conservative adaptive projection method for the variable density incompressiblenavierstokes equations , Journal of Computational Physics, 142 (1998), pp. 1–46.[3]
Douglas N Arnold, Franco Brezzi, and Michel Fortin , A stable finite element for thestokes equations , Calcolo, 21 (1984), pp. 337–344.[4]
P. Atzberger and G. Tabak , Systematic stochastic reduction of inertial fluid-structure in-teractions subject to thermal fluctuations , (submitted), (2013).[5]
P. J. Atzberger , Stochastic euler-lagrange methods for fluid-structure interactions with ther-mal fluctuations and shear boundary conditions , arXiv:0910.5739, (2009).[6]
Paul J. Atzberger , Stochastic eulerian lagrangian methods for fluidstructure interactionswith thermal fluctuations , Journal of Computational Physics, 230 (2011), pp. 2821–2837.[7]
P. J. Atzberger, P. R. Kramer, and C. S. Peskin , A stochastic immersed boundarymethod for fluid-structure dynamics at microscopic length scales , Journal of Computa-tional Physics, 224 (2007), pp. 1255–1292.[8]
Ivo Babuˇska , The finite element method with lagrangian multipliers , Numerische Mathematik,20 (1973), pp. 179–192.[9]
John B. Bell, Alejandro L. Garcia, and Sarah A. Williams , Numerical methods for thestochastic landau-lifshitz navier-stokes equations , Phys. Rev. E, 76 (2007), pp. 016708–.[10]
J. F. Brady and G. Bossis , Stokesian dynamics , Annual review of fluid mechanics.Vol.20—Annual review of fluid mechanics. Vol.20, (1988), pp. 111–57.[11]
S. Brenner and L. R. Scott , The Mathematical Theory of Finite Element Methods , Springer,2008.[12]
F. Brezzi , On the existence, uniqueness and approximation of saddle-point problems arisingfrom lagrangian multipliers , ESAIM: Mathematical Modelling and Numerical Analysis, 8(1974).[13]
William L. Briggs, Van Emden Henson, and Steve F. McCormick , A Multigrid Tutorial ,Society for Industrial and Applied Mathematics, 2000.[14]
T. T. Bringley and C. S. Peskin , Validation of a simple method for representing spheres andslender bodies in an immersed boundary method for stokes flow on an unbounded domain ,Journal of Computational Physics, 227 (2008), pp. 5397–5425–.[15]
A. J. Chorin , Numerical solution of navier-stokes equations , Mathematics of Computation,22 (1968), pp. 745–&–.[16]
G. De Fabritiis, M. Serrano, R. Delgado-Buscalioni, and P. V. Coveney , Fluctuatinghydrodynamic modeling of fluids at the nanoscale , Phys. Rev. E, 75 (2007), pp. 026307–. P. PLUNKETT, P.J. ATZBERGER[17]
M. Doi and S. F. Edwards , The Theory of Polymer Dynamics , Oxford University Press, 1986.[18]
A. Donev, A. Vanden-Eijnden, E. Garcia, and J. Bell , On the accuracy of finite-volumeschemes for fluctuating hydrodynamics , Communications in Applied Mathematics andComputer Science, Vol. 5 (2010), p. 149197.[19]
E.H. Dowell and K.C. Hall , Modeling of fluid-structure interaction , Ann. Rev. Fluid Mech.,33 (2001), pp. 445–490.[20]
Howard C Elman, Victoria E Howle, John N Shadid, and Ray S Tuminaro , A parallelblock multi-level preconditioner for the 3d incompressible navierstokes equations , Journalof Computational Physics, 187 (2003), pp. 504 – 523.[21]
D. L. Ermak and J. A. McCammon , Brownian dynamics with hydrodynamic interactions , J.Chem. Phys., 69 (1978), pp. 1352–1360.[22]
B. Froehlea and P. Persson , A high-order discontinuous galerkin method for fluid-structureinteraction with efficient implicit-explicit time stepping , Computer Methods in AppliedMechanics and Engineering (submitted), (2013).[23]
Jonathan Goodman and Alan D. Sokal , Multigrid monte carlo method for lattice fieldtheories , Phys. Rev. Lett., 56 (1986), pp. 1015–1018.[24] ,
Multigrid monte carlo method. conceptual foundations , Phys. Rev. D, 40 (1989),pp. 2035–2071.[25]
B.E. Griffith, X. Luo, D.M. McQueen, and C.S. Peskin. , Simulating the fluid dynamicsof natural and prosthetic heart valves using the immersed boundary method. , Int J ApplMech., 1 (2009), pp. 137–177.[26]
Boyce E. Griffith, Richard D. Hornung, David M. McQueen, and Charles S. Peskin , An adaptive, formally second order accurate version of the immersed boundary method ,Journal of Computational Physics, 223 (2007), pp. 10–49.[27]
Priv-Doz Dr-Ing F Hartmann , The discrete babuˇska-brezzi condition , Ingenieur-Archiv, 56(1986), pp. 221–228.[28]
W. Helfrich , Elastic properties of lipid bilayers: Theory and possible experiments , Z. Natur-forsch., 28c (1973), pp. 693–703.[29]
C. M. Ho and Y. C. Tai , Micro-electro-mechanical-systems (mems) and fluid flows , AnnualReview of Fluid Mechanics, 30 (1998), pp. 579–612.[30]
Louis H. Howell, John, and John B. Bell , An adaptive-mesh projection method for viscousincompressible flow , SIAM J. Sci. Comput, 18 (1996), pp. 18–996.[31]
Adrin J. Lew and Gustavo C. Buscaglia , A discontinuous-galerkin-based immersed bound-ary method , Int. J. Numer. Meth. Engng., 76 (2008), pp. 427–454.[32]
R. Lipowsky , The conformation of membranes , Nature, 349 (1991), pp. 475–481.[33]
M.R. Maxey and B.K. Patel , Localized force representations for particles sedimenting instokes flow , International Journal of Multiphase Flow, 27 (2001), pp. 1603–1626.[34]
Steve F. McCormick , Multilevel Adaptive Methods for Partial Differential Equations , Societyfor Industrial and Applied Mathematics, 1989.[35]
Michael L. Minion , A projection method for locally refined grids , Journal of ComputationalPhysics, 127 (1996), pp. 158–178.[36]
C. S. Peskin , The immersed boundary method , Acta Numerica, 11 (2002), pp. 479–517.[37]
Arvind Raman, John Melcher, and Ryan Tung , Cantilever dynamics in atomic force mi-croscopy , Nano Today, 3 (2008), pp. 20–27.[38]
L. E. Reichl , A Modern Course in Statistical Physics , John Wiley and Sons, 1998.[39]
Alexandre M Roma, Charles S Peskin, and Marsha J Berger , An adaptive version of theimmersed boundary method , Journal of Computational Physics, 153 (1999), pp. 509–534.[40]
Jens Rotne and Stephen Prager , Variational treatment of hydrodynamic interaction inpolymers , J. Chem. Phys., 50 (1969), pp. 4831–4837.[41]
Nitin Sharma and Neelesh A. Patankar , Direct numerical simulation of the brownian mo-tion of particles by using fluctuating hydrodynamic equations , Journal of ComputationalPhysics, 201 (2004), pp. 466–486.[42]
Michael J. Shelley and Jun Zhang , Flapping and bending bodies interacting with fluid flows ,Annu. Rev. Fluid Mech., 43 (2011), pp. 449–465.[43]
Ray S. Tuminaro and Charles Tong , Parallel smoothed aggregation multigrid: Aggregationstrategies on massively parallel machines , in in SuperComputing 2000 Proceedings, 2000.[44]
Petr Vanek, Marian Brezina, and Radek Tezaur , Two-grid method for linear elasticityon unstructured meshes , SIAM Journal on Scientific Computing, 21 (1999), pp. 900–923.[45]
Petr Vanek, Petr Van Ek, Jan Mandel, and Marian Brezina , Convergence of algebraicmultigrid based on smoothed aggregation , Computing, 56 (1998), pp. 179–196.[46]
Petr Vanˇek, Jan Mandel, and Marian Brezina , Algebraic multigrid by smoothed aggrega-tion for second and fourth order elliptic problems , Computing, 56 (1996), pp. 179–196.TOCHASTIC MULTIGRID FOR BROWNIAN-STOKESIAN DYNAMICS [47] Y Wang and P. Atzberger , Fluctuating hydrodynamic methods for fluid-structure interac-tions in confined channel geometries
Yaohong Wang, Jon Karl Sigurdsson, Erik Brandt, and Paul J. Atzberger , Dynamicimplicit-solvent coarse-grained models of lipid bilayer membranes: Fluctuating hydrody-namics thermostat , Phys. Rev. E, 88 (2013), pp. 023301–.[49]
Charles Whitmer , Over-relaxation methods for monte carlo simulations of quadratic andmultiquadratic actions , Phys. Rev. D, 29 (1984), pp. 306–.[50]
Hiromi Yamakawa , Transport properties of polymer chains in dilute solution: Hydrodynamicinteraction , J. Chem. Phys., 53 (1970), pp. 436–443.[51]
Jinyun Yu , Symmetric gaussian quadrature formulae for tetrahedronal regions , ComputerMethods in Applied Mechanics and Engineering, 43 (1984), pp. 349–353.