Stochastic Eulerian-Lagrangian Methods for Fluid-Structure Interactions with Thermal Fluctuations and Shear Boundary Conditions

# Stochastic Eulerian-Lagrangian Methods for Fluid-Structure Interactions with Thermal Fluctuations and Shear Boundary Conditions

Paul J. Atzberger University of California, Department of Mathematics , Santa Barbara, CA 93106; e-mail: atzberg@math.ucsb.edu; phone: 805-893-3239; Work supported by NSF DMS-0635535.
###### Abstract

A computational approach is introduced for the study of the rheological properties of complex fluids and soft materials. The approach allows for a consistent treatment of microstructure elastic mechanics, hydrodynamic coupling, thermal fluctuations, and externally driven shear flows. A mixed description in terms of Eulerian and Lagrangian reference frames is used for the physical system. Microstructure configurations are represented in a Lagrangian reference frame. Conserved quantities, such as momentum of the fluid and microstructures, are represented in an Eulerian reference frame. The mathematical formalism couples these different descriptions using general operators subject to consistency conditions. Thermal fluctuations are taken into account in the formalism by stochastic driving fields introduced in accordance with the principles of statistical mechanics. To study the rheological responses of materials subject to shear, generalized periodic boundary conditions are developed where periodic images are shifted relative to the unit cell to induce shear. Stochastic numerical methods are developed for the formalism. As a demonstration of the methods, results are presented for the shear responses of a polymeric fluid, lipid vesicle fluid, and a gel-like material.

Key words. Statistical Mechanics, Complex Fluids, Soft Materials, Stochastic Eulerian Lagrangian Methods, SELM, Stochastic Immersed Boundary Methods, SIB, Fluctuating Hydrodynamics, Fluid-Structure Coupling, Polymeric Fluid, FENE, Vesicles, Gels.

## 1 Introduction

Soft materials and complex fluids are comprised of microstructures which have mechanics and interactions characterized by energy scales comparable to thermal energy. This feature results in interesting bulk material properties and phenomena which often depend sensitively on temperature and applied stresses [Bird1987Vol1, ; Barnes1997, ; Doi1986, ]. Example materials include colloidal suspensions, foams, polymeric fluids, surfactant solutions, lipid vesicles, and gels  [Lubensky1997, ; Hamley2003, ; bird1987, ; Gompper2006, ; Doi1986, ; Coussot2007, ; Raub2007, ]. Microstructures of such materials include flexible filaments, bubbles, colloidal particles, lipid chains, and polymers. These microstructures are typically surrounded by a solvating fluid which further mediates interactions through solvation shells [Meyer2006, ; Leckband2001, ] and hydrodynamic coupling [Bird1987Vol1, ; Brady1988, ; Doi1986, ]. In addition, given the energy scales of the microstructure mechanics and interactions, thermal fluctuations often play an important role both in microstructure organization and kinetics [bird1987, ; Brady1988, ; Doi1986, ]. A fundamental challenge in the study of soft materials is to relate bulk material properties to microstructure mechanics, interactions, and kinetics.

For the study of soft materials we introduce a modeling and simulation approach which consistently accounts for microstructure elastic mechanics, hydrodynamic coupling, and thermal fluctuations. The modeling approach is based on a mixed Eulerian and Lagrangian description. The microstructure configurations are modeled in a Lagrangian reference frame, while an Eulerian reference frame is used to account for conserved quantities, such as momentum, of the system. When coupling these disparate descriptions an important issue is to formulate methods which do not introduce artifacts into the conservation laws, such as artificial dissipation of energy or loss of momentum. These properties are especially important when introducing stochastic driving fields to account for thermal fluctuations. We discuss a general approach for developing such coupling schemes, focusing primarily on one such realization referred to as the Stochastic Immersed Boundary Method [Atzberger2007a, ; Peskin2002, ].

To facilitate studies of the rheological properties of soft materials we introduce methods to account for externally driven shear flows. To account for shearing of the material, we generalize the usual periodic boundary conditions so that periodic images are shifted relative to the unit cell to induce shear. Our approach is based on boundary conditions introduced for Molecular Dynamics simulations which are referred to as Lees-Edwards boundary conditions [LeesEdwards1972, ]. These conditions present a number of challenges in the context of numerically solving the hydrodynamic equations. We develop numerical methods which utilize jump conditions in the velocity field at domain boundaries and utilize a change of variable to facilitate handling of the shifted boundaries. Further issues arise when accounting for the thermal fluctuations. For the introduced discretizations we develop stochastic driving fields which yield stochastic numerical methods which are consistent with the principles of statistical mechanics.

We consider primarily two physical regimes. In the first, the relaxation dynamics of the fluid modes is explicitly resolved. In the second, the fluid modes are treated as as having relaxed to a quasi-steady-state distribution. In the first regime we develop efficient stochastic numerical methods for the generation of the corresponding fluctuating fields of the fluid. In the second, we develop efficient stochastic numerical methods which account for the correlated stochastic driving fields which account for the effective thermal fluctuations which drive the microstructure dynamics.

As a demonstration of the proposed stochastic numerical methods, simulations are performed for specific systems. These include studying the shear responses of (i) a polymeric fluid, (ii) a vesicle fluid, and (iii) a gel-like material. To relate microstructure interactions and kinetics to bulk material properties we develop estimators for an associated macroscopic stress tensor. The estimators take into account the n-body interactions in the microstructure mechanics and the generalized boundary conditions. For the polymeric fluid, this notion of stress is used to investigate the dependence of the shear viscosity and normal stresses on the rate of shear. The vesicle fluid is subject to oscillating shear and simulations are preformed to characterize the frequency response in terms of the elastic storage modulus and viscous loss modulus over a wide range of frequencies. As a further demonstration of the methods, the time dependent shear viscosity of a gel-like material is studied through simulations.

The ability to simulate explicitly the microstructure mechanics, hydrodynamic coupling, and thermal fluctuations provides an important link between bulk material properties and phenomena at the level of the microstructures. The presented framework and related stochastic numerical methods are expected to be applicable in the modeling and simulation of a wide variety of soft materials and complex fluids. The general approaches introduced for coupling the Eulerian and Lagrangian descriptions and for the incorporation of thermal fluctuations are expected to allow for the development of many different types of Stochastic Eulerian Lagrangian Methods.

## 2 Stochastic Eulerian-Lagrangian Modeling Approach

We use a mixed Eulerian-Lagrangian description. The conserved quantities of the entire system including both the fluid and microstructures will be accounted for in an Eulerian reference frame. The microstructure configurations will be accounted for in a Lagrangian reference frame, see Figure LABEL:figure_EL_schematic. To introduce the basic approach and simplify the presentation we consider here only a rather special case. We shall assume the solvent fluid is an incompressible Newtonian fluid of constant density and the microstructures are density matched with the fluid. In this case, the primary conserved quantity of interest is the local momentum of the material. The basic framework and principles that will be presented are more generally applicable allowing for additional conserved quantities to be taken into account, such as the local mass density and energy [Donev2009, ]. A more abstract and general presentation of the formalism will be the focus of another paper.

The basic Eulerian-Lagrangian formalism describing the state of the fluid and microstructures is given by the following equations

 Dp(x,t)Dt = ∇⋅σ(x,t)+Λ(x,t)+λ(x,t)+g(x,t) \hb@xt@.01(2.1) ∂X(q,t)∂t = Γ(q,t)+γ(q,t)+Z(q,t). \hb@xt@.01(2.2)

The accounts for the momentum of the material occupying location and denotes the material derivative. The denotes the configuration of the material at time and parameterized by . The local material stress is denoted by . We use the convention that accounts only for the dissipative stress contributions in the system. The operators , couple the Eulerian and Lagrangian descriptions of the state of the material. The operator accounts for momentum gained or lost locally in the system as the material deforms from non-dissipative stresses and body forces. The operator determines the rate of deformation of the material from the momentum of the system. The and are Lagrange multipliers associated with time-independent kinematic constraints imposed on the system, such as rigidity of a body or incompressibility. Thermal fluctuations are taken into account through the stochastic fields and .

We consider systems where the total energy is given by

 E[p,X]=∫ρ02|u(x)|2dx+Φ[X], \hb@xt@.01(2.3)

where is the velocity of the material at location , is the constant mass density of the material, and is the potential energy for a given configuration. The force associated with this energy is denoted by .

For the operators which couple the Eulerian and Lagrangian descriptions to be physically consistent, the following should hold: (i) the coupling operators should not introduce any loss or gain of energy, (ii) momentum should only change through forces acting within the system. More precisely, these conditions require

 ∫F(q)⋅Γ(q)dq = ∫[ρ−10p(x)]⋅Λ[F](x)dx \hb@xt@.01(2.4) ∫ΩΛ[F](x)dx = ∫F(q)dq. \hb@xt@.01(2.5)

The conditions are required to hold for any realization of , , and . The condition LABEL:equ_coupling_cond_energy ensures the coupling operators conserve energy. The condition LABEL:equ_coupling_cond_momentum ensures that in the absence of constraints the total momentum change of the system is equal to the total force acting on the system.

In the notation, we find it convenient to write the operator as explicitly depending on both and , which for conservative forces is technically redundant. To simplify the discussion, it has been assumed that the stress contributions denoted by are entirely dissipative and that there is no net in-flux of momentum from boundary stresses .

With these conditions satisfied by the coupling operators, we discuss how to account for thermal fluctuations using the stochastic fields and . It is convenient when accounting for thermal fluctuations to introduce coupling operators so that all configurations are equally probable at statistical steady-state, when the . It can be shown that this corresponds to dynamics determined by the constraints and coupling operators which introduces an incompressible flow on phase space. The requirement of an incompressibile flow on phase space can be expressed as

 ∫δΛδp(x,x)dx+∫δλδp(x,x)dx+∫δΓδX(q,q)dq+∫δγδX(q,q)dq = 0 \hb@xt@.01(2.6)

This condition can be relaxed to allow for more general choices of coordinates, coupling operators, and constraints. If this condition is not satisfied a more general treatment of the thermal fluctuations is required to take into account in the invariant distribution the local compression or dilation of volume under the phase space flow [Tuckerman1999, ].

To simplify the discussion, we assume that the dissipative processes can be accounted for by a negative definite self-adjoint linear operator in , so that , and that conditions LABEL:equ_coupling_cond_energy - LABEL:equ_coupling_cond_uniform are satisfied. With these assumptions, the thermal fluctuations can be accounted for using for and Gaussian stochastic fields which are mean zero and -correlated in time [Oksendal2000, ; Gardiner1985, ]. The main issue then becomes to determine an appropriate spatial covariance structure for these stochastic fields. By requiring that the Boltzmann distribution be invariant under the stochastic dynamics of equations LABEL:equ_EL_E - LABEL:equ_EL_L, it is required that , and that

 G(x,t,y,s)=⟨(g(x,t))(g(y,s))T⟩=−2kBTρ0δ(t−s)Lδ(x−y), \hb@xt@.01(2.7)

see Appendix LABEL:appendix_invariance_Boltzmann.

Similar formulations as equations LABEL:equ_EL_E - LABEL:equ_EL_L, with , , are the starting point for the derivation of a wide variety of computational approaches used for systems in which fluids interact with rigid or elastic bodies. These include Arbitrary Lagrangian-Eulerian Methods (ALE) [Donea1983, ; Donea2004, ], Fluctuating Immersed Material (FIMAT) Dynamics [Patankar2004, ], Immersed Finite Element Methods (IFEM) [Wang2009, ; Zhang2004, ], and Immersed Boundary Methods (IBM) [Atzberger2007a, ; Peskin2002, ]. The approaches we introduce allow for the incorporation and simulation of thermal fluctuations in such methods, which collectively we refer to as Stochastic Eulerian-Lagrangian Methods (SELMs).

## 3 Semi-Discretization of the Momentum Equations, Microstructures, and the Eulerian-Lagrangian Coupling

We now consider semi-discretizations of the SELM equations. The momentum equations will be spatially discretized on a uniform mesh. The will denote the momentum at the mesh site indexed by and the composite vector of such values will be denoted by . The deformation state which describes the microstructure configurations will be discretized using a finite number of degrees of freedom denoted by indexed by and the composite vector denoted by . As an energy for this discretized system we use

 E[p,X]=∑m12ρ−10|pm|2Δxd+Φ(X) \hb@xt@.01(3.1)

where is the mesh spacing and is the number of dimensions. The semi-discretization in space of the momentum and configuration equations can be expressed as

 ~Dp~Dt = Lp+Λ+λ+g \hb@xt@.01(3.2) ∂X[j](t)∂t = Γ[j]+γ[j] \hb@xt@.01(3.3)

where and denote respectively the spatially discretized approximation of the material derivative and . The , , , denote the composite vector of values on the mesh and , , denote values associated with the configurational degree of freedom. We assume the discrete dissipative operator is symmetric and negative semi-definite. For the coupling operators of the discretized equations the corresponding consistency conditions LABEL:equ_coupling_cond_energy - LABEL:equ_coupling_cond_momentum are

 Γ[p]TF = ρ−10pTΛ[F]Δxd \hb@xt@.01(3.4) ∑mΛ[F]mΔxd = ∑jF[j]. \hb@xt@.01(3.5)

The superscript T denotes the matrix transpose. The first condition ensures for the discretized system that the coupling preserves the energy and the second that changes in momentum only occur from forces acting within the system. The phase space incompressibility condition corresponding to equation LABEL:equ_coupling_cond_uniform in the discretized setting becomes

 ∇p⋅(Λ+λ)+∇X⋅(Γ+γ) = 0 \hb@xt@.01(3.6)

This condition ensures the uniform distribution for the configurations is invariant under the stochastic dynamics of equation LABEL:equ_EL_E_discr-LABEL:equ_EL_L_discr when the potential energy is constant, i.e. . When and are linear operators the energy condition LABEL:equ_coupling_cond_energy amounts to the coupling operators being adjoints (up to a scalar),

 Γ=ΛTρ0Δxd. \hb@xt@.01(3.7)

Provided conditions LABEL:equ_coupling_cond_energy_discr - LABEL:equ_coupling_cond_uniform_discr are satisfied, the thermal fluctuations can be taken into account using a Gaussian stochastic field on the lattice, without any direct stochastic forcing required in the microstructure equations LABEL:equ_EL_L. If these conditions are violated by the discretization the numerical approximation may introduce artificial loss or gain of energy or momentum in the system. In order to be consistent with the fluctuation-dissipation principle of statistical mechanics such discretizations would require additional sources of stochastic forcing to obtain the appropriate Boltzmann ensemble.

An appropriate covariance structure for the stochastic driving field for discretizations satisfying conditions LABEL:equ_coupling_cond_energy_discr - LABEL:equ_coupling_cond_uniform_discr can be determined by requiring invariance of the Boltzmann distribution under the stochastic dynamics of equation LABEL:equ_EL_E_discr - LABEL:equ_EL_L_discr. This yields for the semi-discrete system, see Appendix LABEL:appendix_invariance_Boltzmann,

 G=⟨ggT⟩=−2LC. \hb@xt@.01(3.8)

The covariance of the equilibrium fluctuations is given by the entries

 C=ρ0kBTΔxdI, \hb@xt@.01(3.9)

where is the identity matrix, see Appendix LABEL:appendix_invariance_Boltzmann.

One such realization of this SELM approach is the Stochastic Immersed Boundary Method developed in [Atzberger2007a, ]. In this case the coupling operators are given by

 ΛIBF = M∑j=1F[j](X(t))δa(xm−X[j](t)) \hb@xt@.01(3.10) [ΓIBu][j] = ∑mδa(xm−X[j](t))um(t)Δxd. \hb@xt@.01(3.11)

where , and is a special kernel function approximating the Dirac -function, see Appendix LABEL:appendix_delta_func. The dynamics are subject to the constraint that the fluid is incompressible . For the semi-discretized SIB method of [Atzberger2007a, ] the conditions LABEL:equ_coupling_cond_energy_discr and LABEL:equ_coupling_cond_momentum_discr can be readily verified to hold exactly. However, the condition LABEL:equ_coupling_cond_uniform_discr only approximately holds and is exact only in the continuum limit. As we have

 ∇X[j]⋅Γ = ∑m−∇δa(xm−X[j])umΔxd→∫−∇δa(y−X[j])u(y)dy = ∫δa(y−X[j])∇⋅u(y)dy=0.

In the last line we used that the fluid is incompressible . The divergence of the terms , , , are zero in this case. In practice, the method still yields reasonable results since the exhibited fluctuations deviate from Boltzmann statistics only up to the discretization error [Atzberger2007a, ; Atzberger2008, ].

## 4 Soft Materials and Complex Fluids Subject to Shear

We now discuss how the SELM approach can be used for the study of rheological properties of soft materials and complex fluids. We then discuss specific stochastic numerical methods for performing simulations in practice. For the type of materials we consider, it will be assumed that the solvent hydrodynamics is described well by the constitutive laws of Newtonian fluids in the physical regime where the Reynolds number is small. We also assume that the explicitly represented microstructures occupy only a relatively small volume fraction and are effectively density matched with the solvent fluid. In this regime, the momentum of the system will be accounted for using the time-dependent Stokes equations

 ρ0∂u∂t = μΔu−∇p+Λ+g \hb@xt@.01(4.1) ∇⋅u = 0 \hb@xt@.01(4.2)

where is the local velocity of the fluid body at in the Eulerian reference frame, is the fluid density, is the dynamic viscosity, and is the pressure. This corresponds to the dissipative stress and in equation LABEL:equ_EL_E. While the Reynolds number is small, the partial time derivative is retained in the Stokes flow in equation LABEL:equ_time_dep_stokes since the thermal fluctuations introduce small characteristic time scales into the dynamics.

To introduce shear we generalize the usual periodic boundary conditions. Our basic approach is motivated by the molecular dynamics methods introduced by Lees-Edwards [LeesEdwards1972, ; Evans1979, ; EvansMorriss1984, ]. In this work, molecules in the base unit cell have modified interactions with molecules in periodic images. To simulate a bulk material undergoing a shear deformation at a given rate, the periodic images are treated as shifting in time relative to the unit cell, see Figure LABEL:figure_meshLeesEdwardsI. This has the effect of modifying both the location of periodic images of molecules and their assigned velocities. This has some advantages over other approaches, where an affine-like deformation is imposed on the entire material body [EvansMorriss1990, ; Hoover1980, ; Hoover2008, ]. In contrast, for the Lees-Edwards approach the shear deformation is only imposed at the boundaries allowing within the unit cell for the molecular interactions to determine the form of the shear response.

Motivated by this molecular dynamics condition we develop a corresponding methodology for the SELM approach. For momentum accounted for by the time-dependent Stokes equations we introduce the following generalized periodic boundary conditions

 u(x,y,L,t)=u(x−vt,y,0,t)+vex. \hb@xt@.01(4.3)

For concreteness we consider the case where a shear is imposed in the z-direction giving rise to velocities in the x-direction. The is the side length of the periodic cell in the z-direction, is the velocity of the top face of the unit cell relative to the bottom face, denotes the rate of shear deformation, and is the standard unit vector in the direction. The interactions between microstructures of the system can be readily handled in the same manner as in the molecular dynamics simulation. This is done by shifting the location of any microstructure of a periodic image involved in an interaction.

While conceptually straight-forward, these boundary conditions present significant challenges in practice for the numerical discretization of the momentum equations. The conditions introduce both a jump discontinuity at periodic boundaries and a shift which potentially leads to misalignment of discretization nodes at the domain boundaries, see Figure LABEL:figure_meshLeesEdwardsI. For commonly employed approaches such as spectral Fourier methods the jump discontinuity results in a degradation of accuracy through the resulting Gibbs’ phenomena [Goldberg1976, ]. For uniform finite difference methods on the unit cell the mesh misalignment requires modified stencils or interpolations at the domain boundary. When incorporating stochastic driving fields to account for thermal fluctuations these issues are further compounded.

To address these issues, we develop discretization methods which utilize a moving coordinate frame which deforms with the unit cell, see Figure LABEL:figure_meshLeesEdwardsII. Let the velocity field in this frame be denoted by , where parameterizes the deformed unit cell. Let denote the map from the moving coordinate frame to the fixed Eulerian coordinate frame . The time-dependent Stokes equations in the deforming coordinate frame become

 dw(d)dt = ρ−10μ[ed−δd,3˙γtex]T∇2w(d)[ed−δd,3˙γtex]−∇p+F+J \hb@xt@.01(4.4) ∇⋅w − \hb@xt@.01(4.5)

where parameterizes the deformed unit cell, denotes the rate of the shear deformation, the standard basis vector in the direction with . In the notation the parenthesized superscript denotes a vector component and denotes the Kronecker -function. We also use the notational convention

 [∇2w(d)]i,j = ∂2w(d)∂qi∂qj \hb@xt@.01(4.6) [∇w]d,j = ∂w(d)∂qj. \hb@xt@.01(4.7)

In the equations, the terms are introduced to account for the jump introduced by the boundary conditions LABEL:equ_StokesJumpBndCond. This allows in the new coordinate frame for use of the usual periodic boundary conditions

 w(q1,q2,L,t)=w(q1,q2,0,t). \hb@xt@.01(4.8)

We now discuss a discretization for equations LABEL:equ_stokesShearFrame_momentum and LABEL:equ_stokesShearFrame_continuity and the corresponding source terms . The following central finite difference approximations will be used

 ∂w(d)∂qi → w(d)(q+ei)−w(d)(q−ei)2Δx \hb@xt@.01(4.9) ∂2w(d)∂qi∂qj → w(d)(q+ei+ej)−w(d)(q−ei+ej)4Δx2 − w(d)(q+ei−ej)−w(d)(q−ei−ej)4Δx2, i≠j ∂2w(d)∂q2i → w(d)(q+ei)−2w(d)(q)+w(d)(q−ei)Δx2. \hb@xt@.01(4.11)

These approximations are substituted into equations LABEL:equ_components_Laplace_wLABEL:equ_components_grad_w to approximate the operators in equation LABEL:equ_stokesShearFrame_momentumLABEL:equ_stokesShearFrame_continuity.

We remark that the moving coordinate frame makes the description of the momentum field have some features of a Lagrangian frame of reference. We none-the-less retain the Eulerian terminology treating this distinction loosely since the deformation corresponds to a somewhat arbitrary coordinate frame introduced for numerical convenience and does not directly follow from the details of the fluid flow. Our discretization approach shares features with Arbitrary Eulerian Lagrangian (ALE) Methods [Donea1983, ; Donea2004, ].

An important issue when using such deforming reference frames is that the discretization stencils may become excessively distorted [Donea1983, ; Donea2004, ]. We avoid this issue by exploiting the periodic symmetry of the system in the and directions. Let the displacement in the -direction of the top of the unit cell relative to the bottom of the unit cell be denoted by . For shear rate and cell size the displacement at time is given by . The periodicity in the and directions has the consequence that for any coordinate frame with there is another coordinate frame with which has aligned mesh sites, see Figure LABEL:figure_meshLeesEdwardsII. By adopting the convention that the coordinate frame with is always used when evaluating stencils the distortion is controlled.

To obtain approximations for the source terms the discretization stencils are applied at the shear boundaries of the unit cell. For any stencil weights involving values at mesh sites which cross the boundary the modified image value is used . The contributions of the stencil weights multiplied by are collected over all boundary mesh sites to obtain the source terms , . This allows for the usual finite difference stencils to be used on the unit cell with regular periodic boundary conditions. When including the source terms this gives the equivalent result of imposing the jump boundary condition LABEL:equ_StokesJumpBndCond. This formulation has a number of advantages when numerically solving the discretized equations and when introducing thermal fluctuations.

The Stokes equations LABEL:equ_stokesShearFrame_momentumLABEL:equ_stokesShearFrame_continuity discretized in this manner on a uniform periodic mesh can be expressed as

 ∂w∂t = L(t)w+F+J+g \hb@xt@.01(4.12) D(t)w = K \hb@xt@.01(4.13)

where denotes the finite difference operator approximating the Laplacian in the moving coordinate frame, and the approximation of the Divergence operator in the moving frame, see equation LABEL:equ_stokesShearFrame_continuity. The discretization can be shown to be symmetric and negative semi-definite for each .

An important property of is that for each time the corresponding stencils are translation invariant with respect to lattice shifts of the mesh. This has the important consequence that the matrix representations are circulant and therefore diagonalizable by Fast Fourier Transforms [Strang1988, ]. As a result, the incompressibility constraint can be handled using FFTs to obtain an exact projection method [Chorin1968, ]. This allows for the discretized approximation of the Stokes equations to be expressed as

 ∂w∂t = ℘(t)[L(t)w+F+J]+g \hb@xt@.01(4.14)

where is the operator which projects to the null space of . The incompressibility condition is then satisfied for all time provided .

We discuss stochastic numerical methods for two particular physical regimes: (i) the relaxation of the hydrodynamic modes of the system is resolved explicitly, (ii) for the current configuration the hydrodynamic modes are treated as having relaxed to statistical steady-state. We remark that the case of resolving the hydrodynamic relaxation of the system is amenable to stochastic numerical methods similar to those introduced in [Atzberger2007a, ]. We discuss this case only briefly and focus primarily on the newly introduced stochastic numerical methods for handling the second case.

## 5 Regime I : Resolution of Hydrodynamic Relaxation

We now discuss in practice how the stochastic fields may be generated in the regime where the relaxation of the hydrodynamic modes is resolved explicitly. For this purpose we express equation LABEL:equ_dw_dt in differential form

 dw = ℘(t)[L(t)w+F+J]dt+QdBt. \hb@xt@.01(5.1)

The denotes the stochastic driving field accounting for thermal fluctuations corresponding to and denotes the composite vector of a standard Brownian motion process at each of the mesh sites. Throughout our discussion the stochastic differential equations will be given the Ito interpretation [Gardiner1985, ].

Using , we see that denotes a matrix square-root of the covariance of the stochastic driving field . Given the discretizations introduced in Section LABEL:sec_soft_materials_subject_shear, the dissipative operator depends on time, which requires, see Appendix LABEL:appendix_fluct_dissip,

 G(t)=−2℘(t)L(t)C. \hb@xt@.01(5.2)

This has the consequence that the covariance of the stochastic driving field is time dependent.

In the discretized system the numerical stencils dependent on time, However, since the shear deformation is volume preserving the discretized summation introduced to model the kinetic energy of the discrete system in equation LABEL:equ_SELM_energy_discr evaluated in the deformed coordinates is in fact not dependent on time. Since we made this choice, we have the important consequence that the Boltzmann equilibrium fluctuations of the velocity field associated with this energy are stationary (independent of time). In other words, the covariance of the equilibrium fluctuations on the discretized lattice for the energy given in equation LABEL:equ_SELM_energy_discr is independent of time, . This holds even though the underlying discretization and corresponding operators , , depend on time.

To obtain an explicit form for we need to compute taking into account the incompressibility constraint LABEL:equ_incomp_time_op. The equilibrium covariance under these constraints is given by

 C=23kBTρ0ΔxdI. \hb@xt@.01(5.3)

The factor arises from application in Fourier space of the projection operator which equivalently enforces the incompressibility. The factor appears in the denominator since the velocity is considered, instead of the momentum . The notation for the stochastic driving field is used loosely when switching between the momentum and velocity equations.

The time dependent covariance structure of the stochastic driving field in equation LABEL:equ_timeDependStokesStoch is of the form . An important issue is whether this will indeed yield a consistent treatment of the thermal fluctuations so that the resulting stochastic dynamical system has the required equilibrium fluctuations. We establish a Fluctuation-Dissipation principle for such time dependent systems in Appendix LABEL:appendix_fluct_dissip.

## 6 Generating the Stochastic Driving Field I

In order for equation LABEL:equ_time_G and LABEL:equ_barC to be useful in practice, we must have efficient methods by which to generate the stochastic driving fields with the required covariance structure. A significant challenge in practice is to generate efficiently the Gaussian stochastic driving field with the required covariance structure . A commonly used approach is to generate a variate with uncorrelated standard Gaussian components and set for an appropriately chosen matrix . The resulting variate then has covariance , with a proper choice of .

However, to carry this out in practice encounters two challenges: (i) given the factor must be determined, (ii) the matrix-vector multiplication must be carried out. For (i) the Cholesky algorithm is typically used with a computational cost of , where is the number of components of . For (ii) the resulting factors are generally not sparse, which when generating each variate incurs a computational cost of . To get a sense of the costs, for a three dimensional mesh, the number of components of scales cubically as , where is the domain size and is the mesh resolution. The associated costs for generating the variates using this approach even for moderate spatial resolutions is prohibitively expensive.

To obtain a more efficient computational method we use specific features of the discretization introduced in Section LABEL:sec_soft_materials_subject_shear. One useful feature of the discretization we use is that the equilibrium covariance matrix is proportional to the identity matrix with . This allows equation LABEL:equ_time_G to be expressed as

 G(t)=−2α℘(t)L(t). \hb@xt@.01(6.1)

We also use the following specific properties of the operators , , and obtained from the discretization. The first is that each of the operators corresponds to use of numerical stencils which are translation invariant on the mesh. This has the important consequence that all of these operators are diagonalizable in the Fourier basis. This has the further important consequence that all of these operators commute. The second is that is an exact projection operator, so that and . Finally, we use that the discrete approximation of the Laplacian is symmetric negative semi-definite so that it can be factored as for some matrix .

By using these properties of the operators we can express the covariance of the stochastic driving field as

 G(t)=(√2α℘U(t))(√2α℘U(t))T. \hb@xt@.01(6.2)

In this form the required matrix square-root is readily obtained as . We remark this is different than the Cholesky factor obtained from which is required to be lower triangular [Trefethen1997, ]. Since the operators and are diagonalizable in Fourier space, the matrix action of the operators and on any vector can be computed using the Fast Fourier Transform with a cost of . In summary, our method allows in practice for the random variates of the stochastic driving field to be computed from very efficiently, with a computational cost of only . This is in contrast to the traditional Cholesky approach with a computational cost of .

## 7 Regime II : Under-resolution of Hydrodynamic Relaxation (Quasi-Steady-State Limit)

For many problems the equations of motion can be simplified by exploiting a separation of time-scales between the time-scale on which the hydrodynamic modes relax to a statistical steady-state and the time-scale associated with the motion of the microstructures. In this case the fluid equations can be approximated by

 w=−~L(t)−1[Λ+J]+a. \hb@xt@.01(7.1)

The and the inverse is defined for the operator restricted to the linear space . The term is introduced to account for the thermal fluctuations in this regime. We refer to this as the Quasi-Steady-State Stokes approximation [Brady1988, ]. Using this in equation LABEL:equ_EL_L_discr, we obtain the following closed system of equations for the motion of the microstructures

 dX(t)dt=H\tiny SELM(t)[F]+¯J+A \hb@xt@.01(7.2)

where

 H\tiny SELM(t) = −Γ~L(t)−1Λ \hb@xt@.01(7.3) ¯J = −Γ~L(t)−1J. \hb@xt@.01(7.4)

The will be used to account for the thermal fluctuations. We consider the specific case when the operator is linear in and is linear in . In this case is a tensor which we refer to as the ”effective hydrodynamic coupling tensor.”

In this regime the thermal fluctuations arise from the hydrodynamic modes which are relaxed to statistical steady-state. A key challenge is to determine the appropriate statistics of which accounts for the time integrated thermal fluctuations of the hydrodynamics which impact the microstructure dynamics. For this purpose we rewrite equation LABEL:equ_X_closed_form in differential form as

 \hb@xt@.01(7.5)

neglecting for the moment , and representing the contributions of by . We derive the covariance structure by requiring consistency with the principle of Detailed-Balance of statistical mechanics [Reichl1998, ]. The Fokker-Planck equation associated with equation LABEL:equ_X_diff_form is

 ∂Ψ(X,t)∂t = −∇⋅J \hb@xt@.01(7.6) J = H\tiny SELM(t)FΨ−12S(t)∇XΨ. \hb@xt@.01(7.7)

The is the probability density for the microstructures to have configuration at time . The equilibrium fluctuations of the system are required to have the Boltzmann distribution

 ΨBD(X)=1Zexp(−Φ(X)/kBT) \hb@xt@.01(7.8)

where is a normalization constant which ensures the distribution integrates to one [Reichl1998, ]. Substituting this above and using that is linear in gives

 J = (H\tiny SELM(t)−12kBTS(t))FΨBD \hb@xt@.01(7.9)

where . The principle of Detailed-Balance requires at thermodynamic equilibrium that . Requiring this to hold for all possible gives

 S(t)=2kBTH\tiny SELM(t). \hb@xt@.01(7.10)

For to provide a covariance for a real-valued stochastic driving term, the hydrodynamic coupling tensor must be symmetric and positive semi-definite. In the case that and are linear operators this is ensured by condition LABEL:equ_coupling_cond_energy_discr, which from expression LABEL:equ_H_el_def gives

 qTH\tiny SELM(t)q=−vT(~L(t)−1)vΔxd≥0. \hb@xt@.01(7.11)

To obtain this result we let and use that is symmetric negative definite. To obtain an approach useful in practice requires efficient methods for the generation of the stochastic driving term with covariance .

### 7.1 Generating the Stochastic Driving Field II

As discussed in Section LABEL:sec_gen_stoch_I, a significant challenge in practice is to generate efficiently the Gaussian stochastic driving terms with the required covariance structure. We discuss an approach for SELM methods when the coupling operators and are linear. In this case

 H\tiny SELM(t) = −Γ~L(t)−1ΓTΔxd \hb@xt@.01(7.12)

by condition LABEL:equ_coupling_cond_energy_discr. Using properties of the operators discussed in Section LABEL:sec_soft_materials_subject_shear, we can express the hydrodynamic coupling tensor as

 H\tiny SELM(t) = (Γ(t)V(t)Δxd/2)(Γ(t)V(t)Δxd/2)T. \hb@xt@.01(7.13)

We have used that the operators and commute and since is an exact projection that , . Since is symmetric negative definite on the linear space , we can factor . The factor is readily obtained since and are diagonalizable in the Fourier basis. From equations LABEL:equ_H_factor and LABEL:equ_def_S_DB we can factor the covariance as with

 R(t)=(2kBTΔxd)1/2Γ(t)V(t). \hb@xt@.01(7.14)

This expression for the factor can be used to compute the required Gaussian stochastic driving term with a computational cost of , where is the total number of mesh sites in the momentum field discretization and assuming the action of can be computed with a cost of with .

The random variates are generated by utilizing the underlying discretization mesh of the momentum equations. This is accomplished by generating on the mesh uncorrelated standard Gaussian random variates . Since is diagonal in the Fourier basis, the action is computed in Fourier space with a cost of only . The operator is then applied. If the operator makes use of only localized values of the mesh it can be computed with computational cost of . The last step in generating the random variate requires a scalar multiplication which incurs a computational cost of . This procedure generates the stochastic driving term with a computational cost of . For a sufficiently large number of microstructure degrees of freedom , this method is significantly more efficient than the traditional approach based on Cholesky factorization of which costs .

#### 7.1.1 Effective Hydrodynamic Coupling Tensor : H\tiny SELM

We now discuss an approach for analyzing the effective hydrodynamic coupling tensors which appear in the quasi-steady-state formulation of the SELM approach. From equation LABEL:equ_H_el_def many types of hydrodynamic coupling tensors are possible depending on the kinetic constraints and choice of coupling operators and . For concreteness we discuss the specific case corresponding to the Stochastic Immersed Boundary Method (SIB) [Peskin2002, ; Atzberger2007a, ]. In the case of the SIB method, the specific coupling operators and are given by LABEL:equ_op_SIB_Lambda and LABEL:equ_op_SIB_Gamma. From equation LABEL:equ_H_el_def the effective hydrodynamic coupling tensor is given by

 \hb@xt@.01(7.15) [H\tiny IB(t)[F]][j] = −∑mδa(xm−X[j](t))[~L(t)−1(M∑j=1F[j]δa(xm−X[j](t)))]mΔxd.

In the notation, the superscript denotes for the composite vector the components associated with the microstructure degree of freedom. The denotes the vector components associated with the mesh site with index . An analysis of variants of this tensor for point particles and slender bodies was carried-out in [Bringley2008, ; Atzberger2007c, ].

Since is linear in the microstructure forces, without loss of generality we can consider the case of only two microstructure degrees of freedom. We denote these as , and the displacement vector by . In making comparisons with other hydrodynamic coupling tensors we find it helpful to make use of approximate symmetries satisfied by . From equation LABEL:equ_def_H_IB, depends on up to a shift of relative to the nearest mesh site, and is similarly rotationally symmetry about the axis of . This allows the tensor components for all configurations to be related to a canonical configuration with . For any configuration this is accomplished by introducing the rotation matrix so that and considering . In our comparisons we consider , where the average is taken over all rotations and shifts with respect to the nearest mesh site.

In practice, to numerically compute we sample random configurations of and . A useful expression for the tensor components is . In this notation, are the standard basis vectors in direction and is the microstructure velocity. For a computational implementation of the SELM method, this can be used by applying the force to the microstructure degrees of freedom and measuring the components of the realized microstructure velocities .

When using a SELM approach the hydrodynamic coupling tensor has features which depend on the discretization of the momentum equations, discretization of the microstructures, and the specific choice of coupling operators. For the specific choice of the IB coupling operators and discretization on a uniform mesh we discuss how the effective hydrodynamic coupling tensor compares with other hydrodynamic coupling tensors. We consider two specific tensors, the Oseen Tensor [Brady1988, ] and the Rotne-Prager-Yamakawa Tensor [Rotne1969, ; Yamakawa1970, ]. The Oseen Tensor for a pair of particles experiencing equal and opposite forces can be expressed in terms of the displacement vector as

 H\tiny{OS}(z) =

Similarly, the Rotne-Prager-Yamakawa Tensor can be expressed in terms of the displacement vector as

 H\tiny{RPY}(z) =

In the notation, denotes the dynamic fluid viscosity, and denotes the effective particle size in terms of the radius of a sphere.

In Figure LABEL:figure_compare_IB_RPY_OS the is compared with the Oseen Tensor and Rotne-Prager-Yamakawa Tensor . It is found that the effective hydrodynamic coupling tensor of the Immersed Boundary Method agrees well with both of the tensors in the far-field . An interesting finding is that in the near-field shows very close agreement to , see inset in Figure LABEL:figure_compare_IB_RPY_OS.

## 8 Applications

The SELM approach is expected to be applicable in the study of many different types of complex fluids and soft materials. As a demonstration of the proposed stochastic numerical methods, simulation studies are carried out for a few specific systems. These include studying: (i) the dependence of the shear viscosity on the shear rate in a FENE polymeric fluid, (ii) the frequency response of the elastic storage modulus and viscous loss modulus of a lipid vesicle fluid subject to oscillatory shear, (iii) the rheological responses over time of a gel-like material subject to a constant rate of shear. We now discuss each of these simulation studies in detail.

### 8.1 Estimating Effective Macroscopic Stress

An important challenge in the study of complex fluids and soft materials is to relate bulk material properties to phenomena on the level of the microstructures of the material. To characterize properties of a material, experimental measurements are often made as a sample of material is subject to shear [Bird1987Vol1, ; Bird1987Vol2, ]. To link microstructure mechanics, interactions, and kinetics to macroscopic material properties we develop estimators for an effective macroscopic stress tensor. Our estimators are based on similar approaches used to obtain the Irving-Kirkwood-Kramer formulas [Irving1950, ; Doi1986, ; Bird1987Vol1, ; Bird1987Vol2, ].

When using the SELM approach, the microstructures are modeled using n-body interactions and the domain is subject to generalized boundary conditions. For example, two body interactions can arise from bonds between monomer particles and three body interactions can arise from bond angle terms included in the potential energy. Estimators for the stress must take these features into account.

To obtain a notion of macroscopic stress we define a normal direction and a plane which cuts the unit cell. We then determine on average the forces exerted by the particles which lie above this plane on the particles which lie below this plane. We define the effective stress associated with this plane as the total of this exerted force divided by the area of the plane. To define an effective macroscopic stress we average over all possible planes within the unit cell having the specified normal direction, see Figure LABEL:figure_stressEstimator.

More precisely, the effective macroscopic stress arising from n-body interactions is estimated using

 σ(n)ℓ,z=1L⟨∫baΛ(n)ℓ,z(ζ)dζ⟩. \hb@xt@.01(8.1)

The is the length of the domain in the -direction and denotes averaging over the ensemble. The denotes the microscopic stress arising from the -body interactions associated with a given stress plane and is defined by

 Λ(n)ℓ,z(ζ) = 1A∑q∈Qnn−1∑k=1k∑j=1f(ℓ)q,jk∏j=1H(ζ−x(z)qj)n∏j=k+1H(x(z)qj−ζ). \hb@xt@.01(8.2)

The is the set of -tuple indices describing the -body interactions of the system, denotes the force acting on the particle of the interaction, and denotes the particle involved in the interaction. As a matter of convention in the indexing we require that implies . This expression corresponds to a sum over all the forces exerted by particles of the material above the cross-section at on the particles of the material below. Each term of the summation over corresponds to a specific number of particles of the -body interaction lying below the cross-section at , see Figure LABEL:figure_stressEstimator.

When integrating the microscopic stress, a useful identity is that

 ∫baΠkj=1H(ζ−x(z)qj)⋅Πnj=k+1H(x(z)qj−ζ)dζ=x∗,(z)qk+1−x∗,(z)qk \hb@xt@.01(8.3)

where

 x∗,(z)qj=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩b,if % x(z)qj≥bx(z)qj,if a≤x(z)qj≤ba,if x(z)qj≤a. \hb@xt@.01(8.4)

By integrating equation LABEL:equ_micro_stress_estimator we obtain

 ∫baΛ(n)(ℓ),z(ζ)dζ=1A∑q∈Qnn−1∑k=1k∑j=1f(ℓ)q,j⋅(x∗,(z)qk+1−x∗,(z)qk). \hb@xt@.01(8.5)

This can be further simplified by switching the order of summation of and and using the telescoping property of the summation over . From equation LABEL:equ_micro_stress_estimator this yields the following estimator for the stress contributions of the -body interactions

 σ(n)ℓ,z=1AL∑q∈Qnn−1∑j=1⟨f(ℓ)q,j⋅(x∗,(z)qn−x∗,(z)qj)⟩. \hb@xt@.01(8.6)

This defines an effective macroscopic stress tensor contribution in terms of the n-body interactions of the microstructures of the material. To obtain the total contribution of the microstructure interactions to the stress, all of the contributions of the n-body interactions are summed to obtain the effective macroscopic stress tensor

 σℓ,z=∑nσ(n)ℓ,z. \hb@xt@.01(8.7)

This notion of the macroscopic stress will be used to link bulk rheological properties to the microscopic simulations.

### 8.2 Application I: Complex Fluid of Finite Extensible Non-linear Elastic (FENE) Dimers

As a demonstration of the proposed computational methodology we consider a fluid with microstructures consisting of elastic polymers. The polymers are modeled as idealized elastic dimers which have the potential energy

 ϕ(r)=12Kr20log(1−(rr0)2). \hb@xt@.01(8.8)

The denotes the polymer stiffness, denotes the length of extension of the dimer, and denotes the maximum permitted extension length [Bird1987Vol2, ]. The configuration of each dimer will be represented using two degrees of freedom , . The potential energy for the dimer is given by , where is the composite vector for the particle configuration.

When the polymeric fluid is subject to shear the thermally fluctuating polymeric microstructures are expected to significantly re-orient and deform as a consequence of the shear stresses. This along with thermal fluctuations of the microstructures is expected to play an important role in the bulk response of the polymeric fluid. To link the bulk material properties of the fluid to the microstructures, we use the effective macroscopic stress obtained from equation LABEL:equ_stress_estimator. To characterize the bulk rheological response we consider the shear viscosity and first normal stress coefficient of the polymeric fluid. We define these as [Bird1987Vol1, ; Bird1987Vol2, ]

 ηp = σ(s,v)p/˙γ \hb@xt@.01(8.9) Ψ1 = (σ(s,s)p−σ(v,v))/˙γ2. \hb@xt@.01(8.10)

The is the rate of shear of the polymeric fluid. In the notation, the superscript indicates the tensor component with the index corresponding to the direction of shear and the index corresponding to the direction of the fluid velocity. The contributions of the solvent fluid to the shear viscosity and normal stresses can be considered separately [Bird1987Vol2, ].

The SELM approach is used to study how the shear viscosity and first normal stress difference depend on the rate of shear of the polymeric fluid. Simulations are performed using the SELM method in the regime where the hydrodynamic modes are relaxed to statistical steady-state with parameters given in Table LABEL:table_FENE_param_descr. For and the coupling tensors of equation LABEL:equ_op_SIB_Lambda and equation LABEL:equ_op_SIB_Gamma are used. From an ensemble average over many computational experiments the moments of the extension vector are estimated as the shear rate is increased. The polymeric microstructure moments are seen to respond strongly as the shear stresses of the fluid increase, see Figure LABEL:figure_compare_FENE_moments. This indicates that the rheological properties of the polymeric fluid will depend significantly on the rate of shear. The SELM simulations show that the shear viscosity and the first normal stress difference do in fact vary significantly with the shear rate, see Figure LABEL:figure_FENE_rheology.

The shear viscosity is found to decrease as the shear rate increases. This appears to occur as a consequence of the dimers increasingly aligning with the direction of the fluid flow and as a consequence of the dimers approaching the maximal extension permitted by equation LABEL:equ_V_r. The increased extension results in a non-linear increase in the effective stiffness of the dimer (defined for a given extension by Taylor expanding to second order equation LABEL:equ_V_r). While the dimers become increasingly extended with stronger restoring forces this is counter-balanced by the dimers being increasingly stiff and the thermal fluctuations less frequently driving the dimer into configurations crossing the stress plane. The net effect is that the mechanical stress transmitted on average by the dimers in the direction of shear does not increase as the shear rate increases. This results in a lower effective shear viscosity (note the division by in equation LABEL:equ_eta_p_def_FENE). This is a well-known phenomena in polymeric fluids and is referred to as shear thinning. The simulations demonstrate that the SELM approach is capable of capturing at the level of the microstructures such phenomena, see Figure LABEL:figure_FENE_rheology.

### 8.3 Application II: Polymerized Lipid Vesicle Fluid

As a further demonstration of the applicability of the SELM approach we show how the stochastic numerical methods can be used to investigate the bulk material properties of a complex fluid with polymerized vesicle microstructures. We discuss how the methods can be used to compute the response of the complex fluid subject to an oscillating shear flow varied over a wide range of frequencies.

To obtain a triangulated mesh which captures the shape of a vesicle having a spherical geometry we start with an icosahedral which is circumscribed by a sphere of a given radius. We use the faces of the icosahedron as an initial triangulated mesh. To obtain a mesh which better approximates the sphere we bisect the three edges of each triangular face to obtain four sub-triangles. The newly introduced vertices are projected radially outward to the surface of the sphere. The process is then repeated recursively to obtain further refinements of the mesh. This yields a high quality mesh for spherical geometries. A vesicle represented by a mesh obtained using two levels of recursive refinement is shown in Figure LABEL:figure_vesicles_meshAlg.

To account for the mechanics of a polymerized vesicle the following interactions are used for the control points of the mesh

 ϕ1(r,ℓ) = 12K1(r−ℓ)2 \hb@xt@.01(8.11) ϕ2(τ1,τ2) = 12K1|τ1−τ2|2. \hb@xt@.01(8.12)

The denotes the displacement between two control points, denotes a preferred distance between control points, and denotes a normalized displacement vector (tangent vector) between two control points. The energy accounts for the stretching of a bond between two control points beyond its preferred extension. The energy accounts for bending of the surface locally by penalizing the misalignment of tangent vectors.

For a given triangulated mesh of control points the total energy is given by

 Φ[X] = E1[X]+E2[X] \hb@xt@.01(8.13) E1[X] = ∑(i,j)∈Q1ϕ1(rij,ℓij) \hb@xt@.01(8.14) E2[X] = ∑(i,j,k)∈Q2ϕ2(τij,τjk). \hb@xt@.01(8.15)

The denotes the composite vector of control points. The control point is denoted by . The and are index sets defined by the topology of the triangulated mesh.

The first energy term accounts for stretching of the vesicle surface and is computed by summing over all local two body interactions defined by the topology of the triangulated mesh. For the distance between the two points having index and , the energy penalizes deviations from the preferred distance . The preferred distances are defined by the geometry of a spherical reference configuration for the vesicle. To ensure the two body interactions are represented by a unique index in we adopt the convention that .

The second energy term accounts for curvature of the vesicle surface and is computed by summing over all local three body interactions defined by the topology of the triangulated mesh. The energy penalizes the the misalignment of the tangent vectors and . In the set of indices in it is assumed that the point with index is always adjacent to both and . To ensure the three body interactions are represented by a unique index in we adopt the convention that .

To investigate the bulk rheological properties, the complex vesicle fluid is subjected to an oscillatory shear with rate . We consider the dilute regime in which it is sufficient to study a single polymerized vesicle subject to oscillatory shear. To estimate the effective macroscopic stress tensor the tensor is decomposed into contributions from two body and three body interactions

 σℓ,z=σ(2)