Stochastic Eulerian Lagrangian Methods for Fluid-Structure Interactions with Thermal Fluctuations

# Stochastic Eulerian Lagrangian Methods for Fluid-Structure Interactions with Thermal Fluctuations

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 CAREER Grant DMS - 0956210. http://www.math.ucsb.edu/atzberg/
###### Abstract

We present approaches for the study of fluid-structure interactions subject to thermal fluctuations. A mixed mechanical description is utilized combining Eulerian and Lagrangian reference frames. We establish general conditions for operators coupling these descriptions. Stochastic driving fields for the formalism are derived using principles from statistical mechanics. The stochastic differential equations of the formalism are found to exhibit significant stiffness in some physical regimes. To cope with this issue, we derive reduced stochastic differential equations for several physical regimes. We also present stochastic numerical methods for each regime to approximate the fluid-structure dynamics and to generate efficiently the required stochastic driving fields. To validate the methodology in each regime, we perform analysis of the invariant probability distribution of the stochastic dynamics of the fluid-structure formalism. We compare this analysis with results from statistical mechanics. To further demonstrate the applicability of the methodology, we perform computational studies for spherical particles having translational and rotational degrees of freedom. We compare these studies with results from fluid mechanics. The presented approach provides for fluid-structure systems a set of rather general computational methods for treating consistently structure mechanics, hydrodynamic coupling, and thermal fluctuations.

F

luid-Structure Interaction, Statistical Mechanics, Fluid Dynamics, Thermal Fluctuations, Fluctuating Hydrodynamics, Stochastic Eulerian Lagrangian Method, SELM.

## 1 Introduction

The development of analytic and computational approaches for the study of fluid-structure interactions has a rich history. Motivations for past work in this area include the study of aerodynamic oscillations induced in airplane wings and propellers [22, 24], the study of animal locomotion including swimming and insect flight [44, 38, 51], and the study of physiological problems such as blood flow through heart valves [30, 49, 27]. A central challenge in work on these applications has been to develop descriptions which capture essential features of the fluid structure interactions while introducing approximations which facilitate analysis and the development of tractable numerical methods [22, 45]. Many such challenges remain and this area of research is still very active [30, 13, 45, 12]. Recent scientific and technological advances motivate the study of fluid-structure interactions in new physical regimes often involving very small length scales [57, 64, 15, 46]. At sufficiently small length scales thermal fluctuations play an important role and pose additional challenges in the study of fluid-structure systems.

Significant past work has been done on the formulation of descriptions of fluid-structure interactions subject to thermal fluctuations. Many of these analytic and numerical approaches originate from the polymer physics community [19, 23, 52, 11]. To obtain descriptions tractable for analysis and numerical simulation, these approaches typically place an emphasis on approximations which retain only the structure degrees of freedom. This often results in significant simplifications in the descriptions and in significant computational savings. This eliminates the many degrees of freedom associated with the fluid and avoids having to resolve the potentially intricate and stiff stochastic dynamics of the fluid. These approaches have worked especially well for the study of bulk phenomena in free solution and the study of complex fluids and soft materials [52, 19, 37].

Recent applications arising in the sciences and in technological fields present situations in which resolving the dynamics of the fluid may be important and even advantageous both for modeling and computation. This includes modeling the spectroscopic responses of biological materials [65, 28, 43], studying transport in microfluidic and nanofluidic devices [57, 47], and investigating dynamics in biological systems [2, 17]. There are also other motivations for representing the fluid explicitly and resolving its stochastic dynamics. This includes the development of hybrid fluid-particle models in which thermal fluctuations mediate important effects when coupling continuum and particle descriptions [18, 20], the study of hydrodynamic coupling and diffusion in the vicinity of surfaces having complicated geometries [57], and the study of systems in which there are many interacting mechanical structures [7, 50, 49]. To facilitate the development of methods for studying such phenomena in fluid-structure systems, we present a rather general formalism which captures essential features of the coupled stochastic dynamics of the fluid and structures.

To model the fluid-structure system, a mechanical description is utilized involving both Eulerian and Lagrangian reference frames. Such mixed descriptions arise rather naturally, since it is often convenient to describe the structure configurations in a Lagrangian reference frame while it is convenient to describe the fluid in an Eulerian reference frame. In practice, this presents a number of challenges for analysis and numerical studies. A central issue concerns how to couple the descriptions to represent accurately the fluid-structure interactions, while obtaining a coupled description which can be treated efficiently by numerical methods. Another important issue concerns how to account properly for thermal fluctuations in such approximate descriptions. This must be done carefully to be consistent with statistical mechanics. A third issue concerns the development of efficient computational methods. This requires discretizations of the stochastic differential equations and the development of efficient methods for numerical integration and stochastic field generation.

We present a set of approaches to address these issues. The formalism and general conditions for the operators which couple the Eulerian and Lagrangian descriptions are presented in Section 2. We discuss simplified descriptions of the fluid-structure system in different physical regimes in Section 3. A derivation of the stochastic driving fields used to represent the thermal fluctuations is also presented in Section 3. Stochastic numerical methods are discussed for the approximation of the stochastic dynamics and generation of stochastic fields in Sections 4. To validate the methodology, we perform in each regime analysis of the invariant probability distribution of the stochastic dynamics of the fluid-structure formalism. We compare this analysis with results from statistical mechanics in Section 5. To demonstrate the applicability of the methodology, we perform computational studies for spherical particles having translational and rotational degrees of freedom. We compare these computational studies with results from fluid mechanics in Section 6.

It should be mentioned that related computational methods have been introduced for the study of fluid-structure interactions  [50, 3, 10, 12, 63, 45, 67, 33, 41]. In recent papers, significant work also has been done toward incorporating the role of thermal fluctuations [6, 14, 9, 21, 7]. This includes the Stochastic Immersed Boundary Method [6], Fluctuating Immersed Material Dynamics [14], Computational Fluctuating Fluid Dynamics [9, 21], and Accelerated Stokesian Dynamics [7]. The formalism presented here can be regarded in part as a generalization of these approaches. It is expected that many of the presented results can be applied to further justify and validate these methods and to provide further extensions. The formalism presented here provides a rather general framework for the development of computational methods for applications requiring a consistent treatment of structure mechanics, hydrodynamic coupling, and thermal fluctuations.

## 2 Summary of the Stochastic Eulerian Lagrangian Method

We summarize here the Stochastic Eulerian Lagrangian Method, abbreviated as SELM. We present the general formalism and a number of alternative descriptions of the fluid-structure system. In many situations the stochastic differential equations for the full fluid-structure dynamics exhibits stiffness. To cope with this issue and to develop efficient numerical methods, simplified descriptions are discussed for various physical regimes. A more detailed discussion and derivation of SELM and the reduced equations in each of the physical regimes is given in Section 3.

To study the dynamics of fluid-structure interactions in the presence of thermal fluctuations, we utilize a mechanical description involving Eulerian and Lagrangian reference frames. Such mixed descriptions arise rather naturally, since it is often convenient to describe the structure configurations in a Lagrangian reference frame while it is convenient to describe the fluid in an Eulerian reference frame. In principle more general descriptions using other reference frames could also be considered. Descriptions for fluid-structure systems having these features can be described rather generally by the following dynamic equations

 ρdudt = Lu+Λ[Υ(v−Γu)]+λ+f\tiny thm (1) mdvdt = (2) dXdt = v. (3)

The denotes the velocity of the fluid, the uniform fluid density. The denotes the configuration of the structure and the velocity of the structure. The mass of the structure is denoted by . To simplify the presentation we treat here only the case when and are constant, but with some modifications these could also be treated as variable. The are Lagrange multipliers for imposed constraints, such as incompressibility of the fluid or a rigid body constraint of a structure. The operator is used to account for dissipation in the fluid, such as associated with Newtonian fluid stresses [1]. To account for how the fluid and structures are coupled, a few general operators are introduced, .

The linear operators are used to model the fluid-structure coupling. The operator describes how a structure depends on the fluid flow while is a negative definite dissipative operator describing the viscous interactions coupling the structure to the fluid. We assume throughout that this dissipative operator is symmetric, . The linear operator is used to attribute a spatial location for the viscous interactions between the structure and fluid. The linear operators are assumed to have dependence only on the configuration degrees of freedom , . We assume further that does not have any dependence on .

To account for the mechanics of structures, denotes the potential energy of the configuration . The total energy associated with this fluid-structure system is given by

 E[u,v,X] = ∫Ω12ρ|u(y)|2dy+12mv2+Φ[X]. (4)

The first two terms give the kinetic energy of the fluid and structures. The last term gives the potential energy of the structures.

As we shall discuss, it is natural to consider coupling operators and which are adjoint in the sense

 (5)

for any and . The and denote the spaces used to parameterize respectively the structures and the fluid. We denote such an adjoint by or . This adjoint condition can be shown to have the important consequence that the fluid-structure coupling conserves energy when in the inviscid and zero temperature limit.

In practice, the conditions discussed above can be relaxed somewhat. For our present purposes these conditions help simplify the presentation. Each of these operators will be discussed in more detail below.

To account for thermal fluctuations, a random force density is introduced in the fluid equations and in the structure equations. These account for spontaneous changes in the system momentum which occurs as a result of the influence of unresolved microscopic degrees of freedom and unresolved events occurring in the fluid and in the fluid-structure interactions.

The thermal fluctuations consistent with the form of the total energy and relaxation dynamics of the system are taken into account by the introduction of stochastic driving fields in the momentum equations of the fluid and structures. The stochastic driving fields are taken to be Gaussian processes with mean zero and with -correlation in time [54]. By the fluctuation-dissipation principle [54] these have covariances given by

 ⟨f\tiny thm(s)fT\tiny thm% (t)⟩ = −(2kBT)(L−ΛΥΓ)δ(t−s) (6) ⟨F\tiny thm(s)FT\tiny thm% (t)⟩ = (2kBT)Υδ(t−s) (7) ⟨f\tiny thm(s)FT\tiny thm% (t)⟩ = −(2kBT)ΛΥδ(t−s). (8)

We have used that and . We remark that the notation which is used for the covariance operators should be interpreted as the tensor product. This notation is meant to suggest the analogue to the outer-product operation which holds in the discrete setting [5]. A more detailed discussion and derivation of the thermal fluctuations is given in Section 3.

It is important to mention that some care must be taken when using the above formalism in practice and when choosing operators. An important issue concerns the treatment of the material derivative of the fluid, . For stochastic systems the field is often highly irregular and not defined in a point-wise sense, but rather only in the sense of a generalized function (distribution) [16, 40]. This presents issues in how to define the non-linear term arising in the material derivative, which appears to require point-wise values of . For such irregular velocity fields, this also calls into question the applicability of the theorems typically used to derive the differential equations from the conservation laws. For instance, for such velocity fields the fluid material body may no longer exhibit smooth deformations over time.

There are a number of ways to deal with this issue. The first is to consider a regularization of the fluid stresses, which are typically the source of irregularity, see equation 6. This can be motivated by the fact that the fluid stress tensors typically considered in continuum mechanics are expected to become inaccurate at molecular length-scales. Ideally, from molecular models of the fluid the small-length scale (large wave-number) responses of the fluid could be determined and provide a justified regularization. For instance, this could provide an alternative to using responses based on Newtonian stresses for all length-scales. For the SELM formalism, this would simply correspond to using for an alternative to the dissipative operator based on Newtonian stresses. The second more easily implemented approach is simply to work with the linearized material derivative, which still retains many of the essential features of the fluid dynamics and is useful for many applications [6].

In this initial presentation of SELM, we shall take the latter approach and treat . This provides a rather general description of fluid-structure systems which incorporate the role of thermal fluctuations. From this initial formalism of SELM, we shall derive a number of simplified descriptions for various physical regimes. These simplified descriptions for each regime tend to yield less stiff differential equations and have other features making them useful in the development of efficient stochastic numerical methods for the formalism.

### 2.1 Regime I

We now consider the regime in which the full dynamics of the fluid-structure system are retained but reformulated in terms of a field describing the total momentum of the fluid-structure system at a given spatial location. This description is more convenient to work with in practice since it results in simplifications in the stochastic driving fields. For this purpose we define

 p(x,t)=ρu(x,t)+Λ[mv(t)](x). (9)

The operator is used to give the distribution in space of the momentum associated with the structures for given configuration . Using this approach, the fluid-structure dynamics are described by

 dpdt = Lu+Λ[−∇XΦ(X)]+(∇XΛ[mv])⋅v+λ+g\tiny thm (10) mdvdt = (11) dXdt = v (12)

where and . The third term in the first equation arises from the dependence of on the configuration of the structures, . The Lagrange multipliers for imposed constraints are denoted by . For the constraints, we use rather liberally the notation with the Lagrange multipliers denoted here not necessarily assumed to be equal to the previous definition. The stochastic driving fields are again Gaussian with mean zero and -correlation in time [54]. The stochastic driving fields have the covariance structure given by

 ⟨g\tiny thm(s)gT\tiny thm% (t)⟩ = −(2kBT)Lδ(t−s) (13) ⟨F\tiny thm(s)FT\tiny thm% (t)⟩ = (2kBT)Υδ(t−s) (14) ⟨g\tiny thm(s)FT\tiny thm% (t)⟩ = 0. (15)

This formulation has the convenient feature that the stochastic driving fields become independent. This is a consequence of using the field for the total momentum for which the dissipative exchange of momentum between the fluid and structure no longer arises. In the equations for the total momentum, the only source of dissipation remaining occurs from the stresses of the fluid. This approach simplifies the effort required to generate numerically the stochastic driving fields and will be used throughout.

### 2.2 Regime II

We now consider a regime in which the formalism can be simplified significantly. In many situations, inertial effects often play a relatively minor role in the structure dynamics as a consequence of the small mass of the structure relative to the displaced fluid or as a consequence of viscosity of the solvent fluid [31, 42]. In such a regime, the relatively rapid dynamics associated with the momentum of the structures often presents a source of stiffness in numerical calculations. To cope with these issues we consider a reduction of the stochastic dynamics of the system in which the structure momentum is eliminated from the description. In particular, we consider the regime in which . The denotes a length-scale characteristic of the size of the immersed structure and associated flow field of the fluid. In the limit , the fluid-structure dynamics are governed by

 dpdt = ρ−1Lp+Λ[−∇XΦ(X)]+(∇X⋅Λ)kBT+λ+g\tiny thm (16) dXdt = ρ−1Γp+Υ−1[−∇XΦ(X)]+ζ+G\tiny thm. (17)

In the notation . This term arises from the thermal fluctuations associated with the momentum of the structures, which have been eliminated from the description. This term plays an important role in the system when the phase-space dynamics of has an associated vector field which is compressible. When considering the Liouville equation on phase-space, this term accounts for local changes of the phase-space volume which occurs as the configuration of the structure changes under the dynamics of the reduced description [62]. For a more detailed discussion see Section 5.2. The stochastic driving fields are again Gaussian with mean zero and with -correlation in time [54]. The stochastic driving fields have the covariance structure given by

 ⟨g\tiny thm(s)gT\tiny thm% (t)⟩ = −(2kBT)Lδ(t−s) (18) ⟨G\tiny thm(s)GT\tiny thm% (t)⟩ = (2kBT)Υ−1δ(t−s) (19) ⟨g\tiny thm(s)GT\tiny thm% (t)⟩ = 0. (20)

A more detailed discussion and derivation of the equations in this regime is given in Section 3.2.

### 2.3 Regime III

The description of the fluid-structure system can be further simplified by considering the viscous coupling between the fluid and structures in the limit . In this case, the fluid-structure dynamics are given by

 dpdt = ρ−1Lp+Λ[−∇XΦ(X)]+(∇X⋅Λ)kBT+λ+g\tiny thm (21) dXdt = ρ−1Γp (22) ⟨g\tiny thm(s)gT\tiny thm% (t)⟩ = −(2kBT)Lδ(t−s). (23)

In the notation . A more detailed discussion and derivation of the equations in this regime in given in Section 3.3.

### 2.4 Regime IV

The description of the fluid-structure system can be further simplified by considering for the fluid the viscous limit in which . In this regime only the structure dynamics remain and can be shown to be given by

 dXdt = H\tiny SELM[−∇XΦ(X)]+(∇X⋅H\tiny SELM)kBT+h% \tiny thm (24) H\tiny SELM = Γ(−℘L)−1Λ (25) ⟨h\tiny thm(s)hT\tiny thm% (t)⟩ = (2kBT)H\tiny SELMδ(t−s). (26)

The denotes a projection operator imposing constraints, such as incompressibility. The adjoint property and symmetry of yields an operator which is symmetric. A more detailed discussion and derivation of the equations in this regime is given in Section 3.4.

### 2.5 Summary

This gives an overview of the SELM formalism and the associated stochastic differential equations. We remark that each of these regimes were motivated by a rather specific limit. Non-dimensional analysis of the equations can also be carried out and other limits considered to motivate working with such reduced equations. We discuss in more detail the derivation of the reduced equations in each regime in Section 3. We discuss how specific stochastic numerical methods can be developed for the SELM formalism in Section 4. We discuss applications and how the SELM formalism can be used in practice in Section 6.

## 3 Derivations for the Stochastic Eulerian Lagrangian Method

We now discuss formal derivations to motivate the stochastic differential equations used in each of the physical regimes. For this purpose, we do not present the most general derivation of the equations. For brevity, we make simplifying assumptions when convenient.

In the initial formulation of SELM, the fluid-structure system is described by

 ρdudt = Lu+Λ[Υ(v−Γu)]+λ+f\tiny thm (27) mdvdt = (28) dXdt = v. (29)

The notation and operators appearing in these equations has been discussed in detail in Section 2. For these equations, we focus primarily on the motivation for the stochastic driving fields used for the fluid-structure system.

For the thermal fluctuations of the system, we assume Gaussian random fields with mean zero and -correlated in time. For such stochastic fields, the central challenge is to determine an appropriate covariance structure. For this purpose, we use the fluctuation-dissipation principle of statistical mechanics [54, 36]. For linear stochastic differential equations of the form

 dZt=LZtdt+QdBt (30)

the fluctuation-dissipation principle can be expressed as

 G=QQT=−(LC)−(LC)T. (31)

This relates the equilibrium covariance structure of the system to the covariance structure of the stochastic driving field. The operator accounts for the dissipative dynamics of the system. For the equations 27 –  29, the dissipative operators only appear in the momentum equations. This can be shown to have the consequence that there is no thermal forcing in the equation for , this will also be confirmed in Section 5.1. To simplify the presentation, we do not represent explicitly the stochastic dynamics of the structure configuration .

For the fluid-structure system it is convenient to work with the stochastic driving fields by defining

 q = [ρ−1f\tiny thm,m−1F% \tiny thm]T. (32)

The field formally is given by and determined by the covariance structure . This covariance structure is determined by the fluctuation-dissipation principle expressed in equation 31 with

 L = [ρ−1(L−ΛΥΓ)ρ−1ΛΥm−1ΥΓ−m−1Υ] (33) C = [ρ−1kBTI00m−1kBTI]. (34)

The denotes the identity operator. The covariance was obtained by considering the fluctuations at equilibrium. The covariance is easily found since the Gibbs-Boltzmann distribution is a Gaussian with formal density . The is the normalization constant for . The energy is given by equation 4. For this purpose, we need only consider the energy in the case when . This gives the covariance structure

 G = (2kBT)[−ρ−2(L−ΛΥΓ)−m−1ρ−1ΛΥ−m−1ρ−1ΥΓm−2Υ]. (35)

To obtain this result we use that and . From the definition of , it is found the covariance of the stochastic driving fields of SELM are given by equations 6– 8. This provides a description of the thermal fluctuations in the fluid-structure system.

### 3.1 Regime I

It is convenient to reformulate the description of the fluid-structure system in terms of a field for the total momentum of the system associated with spatial location . For this purpose we define

 p(x,t)=ρu(x,t)+Λ[mv(t)](x). (36)

The operator is used to give the distribution in space of the momentum associated with the structures. Using this approach, the fluid-structure dynamics are described by

 dpdt = Lu+Λ[−∇XΦ(X)]+(∇XΛ[mv])⋅v+λ+g\tiny thm (37) mdvdt = (38) dXdt = v (39)

where and . The third term in the first equation arises from the dependence of on the configuration of the structures, .

The thermal fluctuations are taken into account by two stochastic fields and . The covariance of is obtained from

 ⟨g\tiny thmgT\tiny thm⟩ = = (2kBT)(−L+ΛΥΓ−ΛΥΛT−ΛΥΛT+ΛΥΛT) = −(2kBT)L.

This makes use of the adjoint property of the coupling operators .

One particularly convenient feature of this reformulation is that the stochastic driving field and become independent. This can be seen as follows

 ⟨g\tiny thmFT\tiny thm⟩ = ⟨f\tiny thmFT\tiny thm⟩+⟨ΛF\tiny thmFT% \tiny thm⟩ = (2kBT)(−ΛΥ+ΛΥ)=0.

This decoupling of the stochastic driving fields greatly reduces the computational effort to generate the fields with the required covariance structure. This shows the covariance structure of the stochastic driving fields of SELM are given by equations 13– 15.

### 3.2 Regime II

In many situations, inertial effects often play a relatively minor role in the structure dynamics as a consequence of the small mass of the structure relative to the displaced fluid or as a consequence of the large viscosity of the solvent fluid [31, 42]. We consider the regime in which , as discussed in Section 2.2. We shall derive formally reduced stochastic equations in the limit .

For this purpose, we focus primarily on the dynamics of the velocity of the structures . This can be expressed using the notation of Ito Stochastic Differential Equations [48] as

 dVt = −m−1Υ(Vt−Γu+Υ−1∇XΦ(X))dt+m−1(2kBTΥ)1/2dBt. (42)

The denotes throughout the standard Brownian motion on  [48]. To simplify the presentation, we consider the case when . We expect similar results to hold more generally. Treating the other degrees of freedom as fixed, we can solve equation 42 using Ito’s Lemma [48]. The stationary behavior of this stochastic process can be expressed as

 Vt = μ0+∫t−∞e−(t−s)m−1Υm−1(2kBTΥ)1/2dBs (43) μ0 = Γu−Υ−1∇XΦ(X). (44)

As a result of the integrand being deterministic in the Ito Integral, the is a Gaussian process. This has the consequence that the statistics of the process are completely determined by its mean and covariance functions. The mean of the process is given at each time by

 μ(t)=⟨Vt⟩=μ0. (45)

The covariance function can be computed using the Ito Isometry [48] to obtain

 ϕ(|τ|) = ⟨(Vt+τ−μ0)(Vt−μ0)T⟩=kBTm−1e−|τ|m−1ΥI. (46)

In the limit this can be expressed as

 ϕ(|t−s|)=2kBTΥ−1[12m−1Υe−|t−s|m−1Υ]→2kBTΥ−1δ(t−s). (47)

We have used formally as . This suggests the following approximation for the velocity of the structures in equations 11 and 12.

 v(t) → Γu−Υ−1∇XΦ(X)+(2kBTΥ−1)1/2dBtdt. (48)

To approximate the term in the limit appearing in equation 10, a different approach is required. For this purpose, we consider the process . By Ito’s Lemma this satisfies the stochastic differential equation

 dRt = ∇f(Vt)dVt+12dVTt∇2f(Vt)dVt (49)

with the formal substitutions , . To simplify the discussion we consider the case when . We expect similar results can be obtained more generally. In this case

 dRt = −2m−1γ(Rt−kBTI−12(mVtμT0+μ0(mVt)T))dt + m−1(2kBTγ)1/2(mVtdBTt+dBt(mVt)T).

From the form of the energy given by equation 4, the structure momentum has equilibrium distribution . This gives a Gaussian with mean and variance

 ⟨mVt⟩ = 0 (51) ⟨(mVt)(mVt)T⟩ = mkBTI. (52)

By reasoning similar to the arguments used to approximate , this suggests the terms involving do not make a contribution in the limit. This suggests to leading order we have

 dRt = −2m−1γ(Rt−kBTI)dt. (53)

This suggests the approximation in the limit

 Rt=mVtVTt→kBTI. (54)

By substituting this result for in equation 10, we have

 (∇XΛ[mv])⋅v→Tr[∇XΛ](kBT)=(∇X⋅Λ)(kBT). (55)

This establishes formally the reduced SELM description when the mass of the structures becomes negligible. It should be mentioned these approximations can be established more rigorously using a perturbation analysis of the Kolomogorov Equations associated with the stochastic processes [35, 25]. It should also be mentioned that other limits can be considered in which additional drift and stochastic terms arise in the reduced momentum equations. This will be the focus of another paper.

### 3.3 Regime III

The SELM stochastic equations can be further reduced if the effective viscous interactions between the structure and fluid are assumed to become very strong. This corresponds to approximating the reduced stochastic equations of Regime II in the formal limit . By this notation, we mean that all eigenvalues of the symmetric operator uniformly tend to infinity. In this formal limit the terms involving are expected to no longer make a contribution to the dynamics. This motivates the reduced stochastic equations 2123.

### 3.4 Regime IV

The description of the fluid-structure system can be further simplified by considering for the fluid the viscous limit in which . In this regime the fluid adopts a quasi-steady-state behavior with respect to the configuration of the structures and the forces acting on the fluid. In this limit only the structure dynamics remain. By approximating using arguments similar to those used in Regime II for approximating , we can derive the reduced stochastic equations 24– 26.

## 4 Computational Methodology

We now discuss briefly numerical methods for the SELM formalism. For concreteness we consider the specific case in which the fluid is Newtonian and incompressible. For now, the other operators of the SELM formalism will be treated rather generally. This case corresponds to the dissipative operator for the fluid

 Lu=μΔu. (56)

The denotes the Laplacian . The incompressibility of the fluid corresponds to the constraint

 ∇⋅u=0. (57)

This is imposed by the Lagrange multiplier . By the Hodge Decomposition, is given by the gradient of a function with . The can be interpreted as the local pressure of the fluid.

A variety of methods could be used in practice to discretize the SELM formalism, such as Finite Difference Methods, Spectral Methods, and Finite Element Methods [29, 60, 59]. We present here discretizations based on Finite Difference Methods.

### 4.1 Numerical Semi-Discretizations for Incompressible Newtonian Fluid

The Laplacian will be approximated by central differences on a uniform periodic lattice by

 (58)

The denotes the index of the lattice site. The denotes the standard basis vector in three dimensions. The incompressibility of the fluid will be approximated by imposing the constraint

 [D⋅u]m=3∑j=1ujm+ej−ujm−ej2Δx. (59)

The superscripts denote the vector component. In practice, this will be imposed by computing the projection of a vector to the sub-space , where is the total number of lattice sites. We denote this projection operation by

 u=℘u∗. (60)

The semi-discretized equations for SELM to be used in practice are

 dpdt = Lu+Λ[−∇XΦ]+(∇XΛ[mv])⋅v+λ+g\tiny thm (61) dvdt = −Υ[v−Γu]+F\tiny thm (62) dXdt = v. (63)

The component . Each of the operators now appearing are understood to be discretized. We discuss specific discretizations for and in Section 6. To obtain the Lagrange multiplier which imposes incompressibility we use the projection operator and

 λ=−(I−℘)(Lu+Υ[v−Γu]+f\tiny thm) (64)

In this expression, we let for the particular realized values of the fields and .

We remark that in fact the semi-discretized equations of the SELM formalism in this regime can also be given in terms of directly, which may provide a simpler approach in practice. The identity could be used to efficiently generate the required stochastic driving fields in the equations for . We present the reformulation here, since it more directly suggests the semi-discretized equations to be used for the reduced stochastic equations.

For this semi-discretization, we consider a total energy for the system given by

 E[u,v,X]=ρ2∑m|u(xm)|2Δx3m+m2|v|2+Φ[X]. (65)

This is useful in formulating an adjoint condition 5 for the semi-discretized system. This can be derived by considering the requirements on the coupling operators and which ensure the energy is conserved when in the inviscid and zero temperature limit.

To obtain appropriate behaviors for the thermal fluctuations, it is important to develop stochastic driving fields which are tailored to the specific semi-discretizations used in the numerical methods. Once the stochastic driving fields are determined, which is the subject of the next section, the equations can be integrated in time using traditional methods for SDEs, such as the Euler-Maruyama Method or a Stochastic Runge-Kutta Method [34]. More sophisticated integrators in time can also be developed to cope with sources of stiffness, but are beyond the scope of this paper [6]. For each of the reduced equations, similar semi-discretizations can be developed as the one presented above.

### 4.2 Stochastic Driving Fields for Semi-Discretizations

To obtain behaviors consistent with statistical mechanics, it is important stochastic driving fields be used which are tailored to the specific numerical discretization employed [6, 21, 5]. To ensure consistency with statistical mechanics, we will again use the fluctuation-dissipation principle but now apply it to the semi-discretized equations. For each regime, we then discuss the important issues arising in practice concerning the efficient generation of these stochastic driving fields.

### 4.3 Regime I

To obtain the covariance structure for this regime, we apply the fluctuation-dissipation principle as expressed in equation 31 to the semi-discretized equations 61– 63. This gives the covariance

 G=−2LC=(2kBT)⎡⎢⎣−ρ−2Δx−3L000m−2Υ0000⎤⎥⎦. (66)

The factor of arises from the form of the energy for the discretized system which gives covariance for the equilibrium fluctuations of the total momentum , see equation 65. In practice, achieving the covariance associated with the dissipative operator of the fluid is typically the most challenging to generate efficiently. This arises from the large number of lattice sites in the discretization.

One approach is to determine a factor such that the block , subscripts indicate block entry of the matrix. The required random field with covariance is then given by , where is the uncorrelated Gaussian field with the covariance structure . For the discretization used on the uniform periodic mesh, the matrices and are cyclic [58]. This has the important consequence that they are both diagonalizable in the discrete Fourier basis of the lattice. As a result, the field can be generated using the Fast Fourier Transform (FFT) with at most computational steps. In fact, in this special case of the discretization, “random fluxes” at the cell faces can be used to generate the field in computational steps [5]. Other approaches can be used to generate the random fields on non-periodic meshes and on multi-level meshes, see [4, 5].

### 4.4 Regime II

The covariance structure can be found using an approach similar to the one presented in Section 4.3. This gives

 G=(2kBT)[−ρ−2Δx−3L00Υ−1]. (67)

By factoring the covariance matrix in the Fourier basis, the field can be generated using FFTs in at most computational steps.

### 4.5 Regime III

The covariance structure can be found using an approach similar to the one presented in Section 4.3. This gives

 G=(2kBT)[−ρ−2Δx−3L000]. (68)

By factoring the covariance matrix in the Fourier basis, the field can be generated using FFT in at most computational steps.

### 4.6 Regime IV

This regime differs from the others since the fluid momentum and structure momentum are both no longer represented explicitly. Spontaneous changes in the momentum of the system were the primary source of fluctuations in the configuration of the structures in the other regimes. While the momentum is no longer represented explicitly, we can none-the-less use a discretization of the momentum equations to generate efficiently the random fields required in the over-damped dynamics. This is done by expressing the covariance of the stochastic driving field as

 G=(2kBT)H\tiny SELM=(2kBT)(Γ℘(−L)−1℘TΓT). (69)

This makes use of and properties of the specific discretized operators and . In particular, commutativity and the projection operator properties , . Let be a factor so that . Using this factor we can express the covariance as

 G=(√2kBTΓ℘U)(√2kBTΓ℘U)T. (70)

From this expression a matrix square-root of is readily obtained, .

We remark this is different than the Cholesky factor obtained for which is required to be lower triangular [61, 58]. Obtaining such a factor by Cholesky factorization would cost , where is the number of structure degrees of freedom. For the current discretization considered, the operators and are diagonalizable in Fourier space. This has the consequence that the action of the operators and can be computed using FFTs with a cost of . The is the number of lattice sites used to discretize . The stochastic driving field is computed from . This allows for the stochastic driving field to be generated in computational steps, assuming the action can be compute in steps. This is in contrast to using the often non-sparse matrix arising from Cholesky factorization which generates the stochastic field with a cost of . We remark that this approach shares some similarities with the method proposed in [7, 56]. Other methods based on splittings or multigrid can also be utilized to efficiently generate stochastic fields with this required covariance structure, see [4, 5].

## 5 Equilibrium Statistical Mechanics of SELM Dynamics

We now discuss how the SELM formalism and the presented numerical methods capture the equilibrium statistical mechanics of the fluid-structure system. This is done through an analysis of the invariant probability distribution of the stochastic dynamics. For the fluid-structure systems considered, the appropriate probability distribution is given by the Gibbs-Boltzmann distribution

 Ψ\tiny GB(z)=1Zexp[−E(z)/kBT]. (71)

The is the state of the system, is the energy, is Boltzmann’s constant, is the system temperature, and is a normalization constant for the distribution [54]. We show this Gibbs-Boltzmann distribution is the equilibrium distribution of both the full stochastic dynamics and the reduced stochastic dynamics in each physical regime.

We present here both a verification of the invariance of the Gibbs-Boltzmann distribution for the general formalism and for numerical discretizations of the formalism. The verification is rather formal for the undiscretized formalism given technical issues which would need to be addressed for such an infinite dimensional dynamical system. However, the verification is rigorous for the semi-discretization of the formalism, which yields a finite dimensional dynamical system. The latter is likely the most relevant case in practice. Given the nearly identical calculations involved in the verification for the general formalism and its semi-discretizations, we use a notation in which the key differences between the two cases primarily arise in the definition of the energy. In particular, the energy is understood to be given by equation 4 when considering the general SELM formalism and equation 65 when considering semi-discretizations.

### 5.1 Regime I

The stochastic dynamics given by equations 10– 12 is a change-of-variable of the full stochastic dynamics of the SELM formalism given by equations 1– 3. Thus verifying the invariance using the reformulated description is also applicable to equations  1– 3 and vice versa. To verify the invariance in the other regimes, it is convenient to work with the reformulated description given for Regime I. The energy associated with the reformulated description is given by

 E[p,v,X]=12ρ∫Ω|p(y)−Λ[mv](y)|2dy+m2|v|2+Φ[X]. (72)

The energy associated with the semi-discretization is

 (73)

The probability density for the current state of the system under the SELM dynamics is governed by the Fokker-Planck equation

 ∂Ψ∂t=−∇⋅J (74)

with probability flux

 J = ⎡⎢⎣L+Λ+∇XΛ⋅v+λ−Υ−∇XΦ+ζv⎤⎥⎦Ψ−12(∇⋅G)Ψ−12G∇Ψ. (75)

The covariance operator is associated with the Gaussian field by . In this regime, is given by equation 13 or 66. In the notation with the summation convention for repeated indices. To simplify the notation we have suppressed denoting the specific functions on which each of the operators act, see equations 10– 12 for these details.

The requirement that the Gibbs-Boltzmann distribution given by equation 71 be invariant under the stochastic dynamics is equivalent to the distribution yielding . We find it convenient to group terms and express this condition as

 ∇⋅J = A1+A2+∇⋅A3+∇⋅A4=0 (76)

where

 (77) A1 = [(Λ+∇XΛ⋅v+λ1)⋅∇pE+(−∇XΦ+ζ1)⋅∇vE+(v)⋅∇XE](−kBT)−1Ψ\tiny GB A2 = [∇p⋅(Λ+∇XΛ⋅v+λ1)+∇v⋅(−∇XΦ+ζ2)+∇X⋅(v)]Ψ\tiny GB A3 = −12(∇⋅G)Ψ\tiny GB A4 = ⎡⎢ ⎢⎣Lu+λ2+[Gpp∇pE+Gpv∇vE+GpX∇XE](2kBT)−1−Υ+ζ2+[Gvp∇pE+Gvv∇vE+GvX∇XE](2kBT)−1[GXp∇pE+GXv∇vE+GXX∇XE](2kBT)−1⎤⎥ ⎥⎦Ψ\tiny GB.

We assume here that the Lagrange multipliers can be split and to impose the constraints by considering in isolation different terms contributing to the dynamics, see equation 77. This is always possible for linear constraints. The block entries of the covariance operator are denoted by with . For the energy of the discretized system given by equation 4 we have

 ∇pnE = u(xn)Δx3n (78) ∇vqE = ∑mu(xm)⋅(−∇vqΛ[mv]m)Δx3m+mvq (79) ∇XqE = ∑mu(xm)⋅(−∇XqΛ[mv]m)Δx3m+∇XqΦ. (80)

where . Similar expressions for the energy of the undiscretized formalism can be obtained by using the calculus of variations [26].

We now consider and each term . The term can be shown to be the time derivative of the energy when considering only a subset of the contributions to the dynamics. Thus, conservation of the energy under this restricted dynamics would result in being zero. For the SELM formalism, we find by direct substitution of the gradients of given by equations 78– 80 into equation 77 that . When there are constraints, it is important to consider only admissible states . This shows in the inviscid and zero temperature limit of SELM the resulting dynamics are non-dissipative. This property imposes constraints on the coupling operators and can be viewed as a further motivation for the adjoint conditions imposed in equation 5.

The term gives the compressibility of the phase-space flow generated by the non-dissipative dynamics of the SELM formalism. The flow is generated by the vector field on the phase-space . When this term is non-zero there are important implications for the Liouville Theorem and statistical mechanics of the system [62]. For the current regime, we have since in the divergence each component of the vector field is seen to be independent of the variable on which the derivative is computed. This shows in the inviscid and zero temperature limit of SELM, the phase-space flow is incompressible. For the reduced SELM descriptions, we shall see this is not always the case.

The term corresponds to fluxes arising from multiplicative features of the stochastic driving fields. When the covariance has a dependence on the current state of the system, this can result in possible changes in the amplitude and correlations in the fluctuations. These changes can yield asymmetries in the stochastic dynamics which manifest as a net probability flux. In the SELM formalism it is found that in the divergence of each contributing entry is independent of the variable on which the derivative is being computed. This shows for the SELM dynamics there is no such probability fluxes, .

The last term accounts for the fluxes arising from the primarily dissipative dynamics and the stochastic driving fields. This term is calculated by substituting the gradients of the energy given by equation 78– 80 and using the choice of covariance structure given by equations 13 or 66. By direct substitution this term is found to be zero, .

This shows the invariance of the Gibbs-Boltzmann distribution under the SELM dynamics. This provides a rather strong validation of the stochastic driving fields introduced for the SELM formalism. This shows the SELM stochastic dynamics are consist with equilibrium statistical mechanics [54].

### 5.2 Regime II

For the reduced stochastic dynamics given by equations 10– 12, the probability density satisfies the Fokker-Planck equation with the probability flux

 J = [ρ−1L+Λ+(∇X⋅Λ)kBT+λρ−1Γ+Υ−1+ζ]Ψ\tiny GB−12(∇⋅G)Ψ\tiny GB−12G∇Ψ%GB. (81)

The denotes the covariance operator for the stochastic driving fields given by equation 67. The invariance of the Gibbs-Boltzmann distribution requires

 (82) ∇⋅J = A1+A2+∇⋅A3+∇⋅A4=0 A1 = [(Λ+(∇X⋅Λ)kBT+λ1)⋅∇pE+(ρ−1Γ+ζ1)⋅∇XE](−kBT)−1Ψ\tiny GB A2 = [∇p⋅(Λ+(∇X⋅Λ)kBT+λ1)+∇X⋅(ρ−1Γ+ζ1)]Ψ\tiny GB A3 = −12(∇⋅G)Ψ\tiny GB A4 = [(ρ−1L+λ2)+[Gpp∇pE+GpX∇XE](2kBT)−1(Υ−1+ζ2)+[GXp∇pE+GXX∇XE](2kBT)−1]Ψ\tiny GB.

To simplify the notation we have suppressed explicitly denoting the functions on which the operators act, which can be inferred from equation 16– 17. In the current regime and the energy given by equation 65 has gradients given by

 ∇pnE = unΔx3n (83) ∇XqE = ∇XqΦ. (84)

Similar expressions can be obtained for the undiscretized formalism using the calculus of variations [26].

We now consider and . The terms have a similar interpretation as in Section 5.1. The term can be interpreted as the time derivative of the energy when considering only a subset of the contributions to the dynamics. By direct substitution of the gradients given by equation 83– 84, we find . This differs from the non-reduced equations in which this term was zero, see Section 5.1.

The term gives the compressibility of the flow generated by the vector field on the phase-space . For the reduced equations, the phase-space flow has compressibility given by , which in general is no longer zero. However, we have that . This follows from the form of the gradients given by equation 83– 84 and from the properties of and . In particular, that the operators are linear and that they are adjoints in the sense of equation 5.

The term accounts for probability fluxes driven by multiplicative features of the stochastic driving fields. It is found this term is zero . This follows from the divergence in which each entry of is independent of the variable on which the derivative is applied. The term accounts for fluxes arising from the dissipative contributions to the dynamics and the stochastic driving fields. By direct substitution of the gradients given in equation 83– 84, and the choice made for given in equation 67, we find this term is zero, . This establishes and that the Gibbs-Boltzmann distribution is invariant for the SELM dynamics.

### 5.3 Regime III

We now discuss briefly the reduced stochastic dynamics given by equations 21– 23, in which . In this regime, the probability flux is almost identical to equation 81, with any terms involving set to zero. With this substitution, it immediately follows that the Gibbs-Boltzmann distribution is invariant under the SELM dynamics.

### 5.4 Regime IV

In the over-damped regime in which the fluid is no longer explicitly represented, we have the reduced stochastic dynamics given by equations 24– 26. The probability density for the current state of the system is governed by the Fokker-Planck equation with the probability flux

 (85) J = [H\tiny SELM[−∇XΦ]+kBT(∇X⋅H\tiny SELM)]Ψ\tiny GB−12(∇X⋅G)Ψ\tiny GB−12G∇XΨ\tiny GB.

In this regime, the covariance is given by , see Section 3.4. This gives and . Substituting these expressions in equation 85, we find . This establishes for the over-damped regime of the SELM formalism that the Gibbs-Boltzmann distribution is invariant and satisfies detailed balance.

### 5.5 Summary

For the SELM formalism, we have demonstrated in each regime that the Gibbs-Boltzmann distribution is invariant. This shows the SELM formalism yields appropriate behaviors with respect to equilibrium statistical mechanics.

## 6 Applications

To demonstrate how the SELM formalism can be used in practice, we consider spherical particles which have translation and rotational degrees of freedom. We give specific operators representing the coupling of the particles and fluid. We compare this SELM formalism with classical results from fluid mechanics. We should mention that similar approaches for the SELM formalism can be applied much more generally to represent spatially extended structures, such as filaments, membranes, or even deformable bodies. The development of representations and specific coupling operators for these structures will be the focus of future work.

### 6.1 Particles with Rotational and Translational Degrees of Freedom

To describe particles which can exhibit translational and rotational motions, we use the degrees of freedom for the center of mass and for the rotational configuration. To describe the full configuration of a particle, we define the composite vector