Closure method for spatially averaged dynamics of particle chains

Closure method for spatially averaged dynamics of
particle chains

Alexander Panchenko Department of Mathematics, Washington State University, Pullman, WA 99164 Lyudmyla L. Barannyk Department of Mathematics University of Idaho, Moscow, ID 83843  and  Robert P. Gilbert Department of Mathematical Sciences, University of Delaware, Newark, DE 19716

We study the closure problem for continuum balance equations that model mesoscale dynamics of large ODE systems. The underlying microscale model consists of classical Newton equations of particle dynamics. As a mesoscale model we use the balance equations for spatial averages obtained earlier by a number of authors: Murdoch and Bedeaux, Hardy, Noll and others. The momentum balance equation contains a flux (stress), which is given by an exact function of particle positions and velocities. We propose a method for approximating this function by a sequence of operators applied to average density and momentum. The resulting approximate mesoscopic models are systems in closed form. The closed from property allows one to work directly with the mesoscale equaitons without the need to calculate underlying particle trajectories, which is useful for modeling and simulation of large particle systems. The proposed closure method utilizes the theory of ill-posed problems, in particular iterative regularization methods for solving first order linear integral equations. The closed from approximations are obtained in two steps. First, we use Landweber regularization to (approximately) reconstruct the interpolants of relevant microscale quantitites from the average density and momentum. Second, these reconstructions are substituted into the exact formulas for stress. The developed general theory is then applied to non-linear oscillator chains. We conduct a detailed study of the simplest zero-order approximation, and show numerically that it works well as long as fluctuations of velocity are nearly constant.

2000 Mathematics Subject Classification:
82D25, 35B27, 35L75, 37Kxx, 70F10, 70Hxx, 74Q10, 82C21, 82C22

Key Words: FPU chain, particle chain, oscillator chain, upscaling, model reduction, dimension reduction, closure problem,

1. Introduction

In a series of papers, [1], [2], [3], [4], Murdoch and Bedeaux studied continuum mechanical balance equations for mesoscopic space time averages of discrete systems. Earlier work of Irving and Kirkwood [14], Noll [16], and Hardy [13] on closely related topics should be also mentioned here. The fluxes in balance equations (e. g. stress) are given by exact formulas as functions of particle positions and velocities. This is useful for linking microscale dynamics with mesoscale phenomena. However, using these formulas requires a complete knowledge of underlying particle dynamics. Since many particle systems of interest have enormous size, direct simulation of particle trajectories may be intractable. Consequently, it makes sense to look for closed form approximations of fluxes in terms of other mesoscale quantities (e.g., average density and velocity), rather than microscopic variables.

In this paper we address the above closure problem for spatially averaged mesoscale dynamics of large size classical particle chains. The design of the method was influenced by the following considerations.

  1. The quantities of interest are space-time continuum averages, such as density, linear momentum, stress, energy and others. This choice of averages is natural because these quantities are experimentally measurable, and also because of their importance in coupled multiscale simulations involving both continuum and discrete models. In addition, by working directly with space-time averages instead of ensemble averages one can bypass a difficult problem of relating probabilistic and space-time averages.

  2. It is desirable to be able to predict behavior of averages on arbitrary time intervals, no matter how short. This perspective comes from PDE problems, where observation time is often arbitrary and long time behavior is not of interest. When one tracks an ODE systems on an arbitrary time interval, transients may be all that is observed. Therefore, we do not use qualitative theory of ODEs, primarily concerned with describing long time features of dynamics. This significantly decreases the range of available tools. However, the closure problem for mesoscopic PDEs turned out to be a question that can be still answered in a satisfactory way. The methods developed in this fashion can be helpful in situations where long time features are not of interest: modeling transient and short-lived phenomena, working with metastable systems, and dealing with problems for which relaxation times can be hard to estimate.

  3. We consider particle systems with initial conditions that either known precisely, or ar least such that the possible initial positions and velocities are strongly restricted by available a priori information. This is in contrast to statistical mechanics, where uncertainty of initial conditions is a major problem. In this regard we note that our approach makes sense for discrete models of solid-fluid continuum systems, where the smallest relevant length scale is still much larger than a typical intermolecular distance. For other particle systems, our method can be used to run deterministic simulations repeatedly, in order to accumulate statistical information about the underlying probability density.

  4. Because of widespread use of computers in physical and engineering sciences, it is useful to develop theories tailored for computation, rather than ”paper and pencil” modeling. As far as the closure problem is concerned, traditional phenomenological approach to formulating constitutive equations can be subsumed by a more general problem of finding a computational closure method. In particular, a closure method can be realized as an iterative procedure where one inputs the values of the primary variables (e.g. density and velocity) computed at the previous moment of time, and the algorithm generates the flux (e.g. stress) at the next moment. Then primary variables are updated using mesoscopic balance equations, and the process is repeated. In addition, focusing on computing one can obtain unconventional but useful continuum mechanical models. By replacing a simple, but possibly crude, Taylor series truncation with an algorithm we make it harder to obtain exact solutions. Since such solutions are rarely available even for simple classical systems, (e.g. Navier-Stokes equations), this is not a serious drawback. On the positive side, computational closure generally contains an explicit (explicitly computable) link between micro- and mesoscale properties.

  5. An important potential application of closure is development of fast numerical methods for simulating meso-scopic dynamics of particle systems. Mesoscale solvers usually employ coarse meshes with mesh size much larger than a typical interparticle distance. Then the averages would be usually given by their coarse mesh values, while interpolants of microscale quantities are discretized on a fine scale mesh. Consequently, a closure method might consist of two generic blocks: (i) reconstruction on mesoscale mesh thereby a coarse approximations of fine scale quantities are obtained from averages; and (ii) interpolation of the obtained coarse scale discretizations to fine scale.

The closure algorithm developed in the paper is based on iterative regularization methods for solving first kind integral equations. We observe that primary mesoscale averages are related to the interpolants of microscale variables via a linear convolution operator. The kernel of this operator is the ”window function” used in [1] to generate averages. Such integral operators are usually compact. A compact operator may be invertible, but the inverse operator is not continuous. Therefore, the problem of reconstructing microscale quantities from given averages is ill-posed. Such problems are well studied in the literature [12, 6, 9, 15, 18]. A particular method used in the paper for inverting convolutions is Landweber iteration [5], [10]. It is known that if the error in the data tends to zero, the Landweber method produces successive approximations converging to the exact solution. For the merely bounded data error, convergence is replaced by a stopping criterion. This criterion provides the optimal number of iterations needed to approximate the solution with the accuracy proportional to the error in the data. As a consequence, our method has desirable feature: one can improve the approximation quality at the price of increasing the algorithm complexity. This means that predictive capability of the method can be regulated depending on available computing power.

The paper is organized as follows. In Section 2 we describe a general multi-dimensional microscopic model. The equations of motion are classical Newton equations. We limit ourselves to the case of short range interaction forces that may be either conservative or dissipative. The scaling of particle masses and forces reflects a continuum mechanical perspective, that is a family of particle systems of increasing size should represent a hypothetical continuum material. As , the total mass of the system should remain fixed, and the total particle energy should be either fixed, or at least bounded independent of . Next, we recall the main points of averaging theory of Murdoch-Bedeaux and provide mesoscopic balance equations and exact formulas for the stress from [4]. In Section 3 we develop integral approximations of averages, and describe the use of Landweber iterative regularization for approximate reconstruction. Section 4 contains the formulation of the scaled ODE equations of the so-called Fermi-Pasta-Ulam (FPU) chains. In Section 5 we derive closed form mesoscopic continuum equations of chain dynamics. The complexity of these continuum models increases with the order of the iterative deconvolution approximation. Section 6 is devoted to the detailed study of the simplest closed model with , which we call zero-order closure. Essentially, zero-order closure means that the microscopic quantities are replaced by their averages. Such an approximation can work well only for systems with small fluctuations. To quantify fluctuation size we introduce upscaling temperature and the related notion of quasi-isothermal dynamics. For such dynamics, we show how to interpolate averages given by mesoscopic mesh values, in order to initialize approximate particle positions and velocities. The interpolation procedure is problem-specific: it conserves microscopic energy and preserves quasi-isothermal nature of the dynamics. Section 7 contains the results of computational tests. Here we apply our zero-closure algorithm to a Hamiltonian chain with the finite range repulsive potential , decreasing as a power of distance. The results show good agreement of zero-order approximations with the exact stress produced by direct simulations with 10000-80000 particles, provided the initial conditions have small fluctuations. In our example, the initial conditions are such that the upscaling temperature is nearly zero during the observation time. We also demonstrate that increasing fluctuations of initial velocities leads to a considerable increase in the approximation error, indicating that higher order closure algorithms should be used instead of zero-order closure. Applicability of the zero-order closure is further discussed in Section 8. Finally, conclusions are provided in Section 9.

2. Microscale equations and mesoscale spatial averages

2.1. Scaled ODE problems

The starting point is the microscale ODE problem. In this paper we shall work with classical Newton equations of point particle dynamics. The same equations may arise as discretization of the momentum balance equation for continuum systems. Consider a system containing identical particles, denoted by . The mass of each particle is , where is the total mass of the system. Suppose that during the observation time , remain inside a bounded domain in , where is the physical space dimension, usually or . The positions and velocities of particles satisfy a system of ODEs


subject to the initial conditions


Here denotes external forces, such as gravity and confining forces. The interparticle forces , where are pair interaction forces which depend on the relative positions and velocities of the respective particles.

We are interested in investigating asymptotic behavior of the system as . Thus it is convenient to introduce a small parameter


characterzing a typical distance between neighboring particles. As approaches zero, the number of particles goes to infinity, and the distances between neighbors shrink. Consequently, the forces in (2.2) should be properly scaled. The guiding principle for scaling is to make the energy of the system bounded independent of , as . In addition, the energy of the initial conditions should be bounded uniformly in .

As an example of scaling, consider forces generated by a finite range potential and assume that each particle interacts with no more than a fixed number of neighbors (this is the case, e.g., for particle chains with nearest neighbor interaction, where a particle always interacts with two neighbors). The fixed number of interacting neighbors implies that there are about interacting pairs. Assuming also that the system is sufficiently dense, and variations of particle concentrations are not large, we can suppose that a typical distance between interacting particles is on the order . The resulting scaling


makes the potential energy of an isolated system bounded independent of . Kinetic energy will be under control provided the total energy of the initial conditions is bounded independent of . If exterior forces are present, they should be scaled as well.

Remark. Superficially, the system (2.1), (2.2) looks similar to the parameter-dependent ODE systems studied in numerous works on ODE time homogenization (see e g. [17] and references therein). In the problem under study, depends on the system dimension , while in the works on time-homogenization and ODE perturbation theory, the system size is usually fixed as .

2.2. Length scales

We introduce the following length scales:
- macroscopic length scale ;
- microscopic length scale ;
- mesoscopic length scale ,
where is a parameter that characterizes spatial mesoscale resolution. This parameter is chosen based on the desired accuracy, the computational cost requirements, available information about initial conditions and behavior of ODE trajectories etc.

The computational domain is subdivided into mesoscopic cells , , with the side length on the order of . The centers of are the nodes of the meso-mesh. The number of unknowns in the mesoscopic system will be on the order of . For computational efficiency, one should have . This does not mean that is close to one. In fact, it makes sense to keep as small as possible in order to have an additional asymptotic control over the system behavior. Decreasing will in general make computations more expensive.

2.3. Averages and their evolution

To define averages we first select a fast decreasing window function satisfying There are many possible choices of the window function. In the paper we assume, unless otherwise indicated, that is a compactly supported, differentiable on the interior of its support, and non-negative. Next, define

Once the window function is chosen, we can evaluate the averages of various continuum mechanical variables, following [1], [4]. The mesoscopic average density and momentum are given by


The meaning of the above definitions becomes clear if one considers , where is a characteristic function of the unit ball in , and is the volume of the unit ball. Then

The sum in the right hand side gives the number of particles located within distance of at time . Multiplying by we get the total mass of these particles, and dividing by (the volume of -ball) gives the usual particle density.

Differentiating (2.6), (2.7) in , and using the ODEs (2.1), (2.2) one can obtain [1] exact mesoscopic balance equations for all primary variables. For example, for an isolated system with (), mass conservation and momentum balance equations take the form:


The stress [4], where


is the convective stress, and


is the interaction stress. The summation in (2.11) is over all pairs of particles that interact with each other.

Discretizing balance equations on the mesoscopic mesh yields a discrete system of equations, called the meso-system, written for mesh values of and . The dimension of the meso-system is much smaller than the dimension of the original ODE problem. However, at this stage we still have no computational savings, since the meso-system is not closed. This means that mesoscopic fluxes such as (2.10), (2.11) are expressed as functions of the microscopic positions and velocities. To find these positions and velocities, one has to solve the original microscale system (2.1), (2.2). To achieve computational savings we need to replace exact fluxes with approximations that involve only mesoscale quantities. We refer to the procedure of generating such approximations as a closure method. This closure-based approach has much in common with continuum mechanics. The important difference is that the focus is on computing, rather than continuum mechanical style modeling of constitutive equations.

3. Closure via regularized deconvolutions

3.1. Outline

Our approach is based on a simple idea: the integral approximations of primary averages (such as density and velocity) are related to the corresponding microscopic quantities via convolution with the kernel . Therefore, given primary variables we can (approximately) recover the microscopic positions and velocities by numerically inverting convolution operators. The results are inserted into equations for secondary averages (or fluxes), such as stress in the momentum balance. This yields closed form balance equations that can be simulated efficiently on the mesoscopic mesh.

3.2. Integral approximation of discrete averages

To exploit the special structure of primary averages, it is convenient to approximate sums such as


by integrals. Since particle positions are not periodically spaced, (3.1) is not in general a Riemann sum for . To interpret the sum correctly, we introduce interpolants of positions and velocities, associated with the microscopic ODE system (2.1), (2.2). At these interpolants satisfy

where , are points of -periodic rectangular lattice in . At other times,

Then we can rewrite (3.1) as


where denotes the volume (Lebesgue measure) of . Eq. (3.2) is a Riemann sum generated by partitioning into cells of volume centered at . This yields


up to discretization error. Now suppose that the map is invertible for each , that is . Changing the variables in the integral we obtain a generic integral approximation




up to discretization error.

3.3. Regularized deconvolutions

Define an operator by

To simplify exposition, suppose that is injective. For example, a Gaussian produces an injective operator, which is not difficult to check using Fourier transform and uniqueness of analytic continuation. If is injective, then there exists the single-valued inverse operator , that we call the deconvolution operator. Unfortunately, this operator is unbounded, since is compact in . This is the underlying reason for the popular belief that averaging destroys the high-frequency information contained in the microscopic quantities. In fact, this information is still there (the inverse operator exists), but it is difficult to recover in a stable manner, because of unboundedness. This does not make the situation hopeless, as has been recognized for some time. Reconstructing from the knowledge of ) is a classical example of an unstable ill-posed problem (small perturbations of the right hand side may produce large perturbations of the solution). The exact nature of ill-posedness and methods of regularizing the problem are well investigated both analytically and numerically (see, e. g. [6, 9, 15, 18, 11, 12]). Accordingly, we interpret notation as a suitable regularized approximation of the exact operator. Many regularizing techniques are currently available: Tikhonov regularization, iterative methods, reproducing kernel methods, the maximum entropy method, the dynamical system approach and others. It is very fortunate that this vast array of knowledge can be used for the ODE model reduction. On the conceptual level, our approach makes it clear that instability associated with ill-posedness is a fundamental difficulty in the process of closing the continuum mechanical equations.

A family of Landweber iterative deconvolution methods [5], [10] seems to be particularly convenient in the present context. In the simplest version, approximations to the solution of the operator equation


are generated by the formula


The number of iterations plays the role of regularization parameter. In (3.7), denotes the identity operator.

The first three low-order approximations are


4. Microscale particle chain equations

In this section, the general method outlined above is detailed in the case of one-dimensional Hamiltonian chain of oscillators that consists of identical particles with nearest neighbor interaction. The domain is an interval . Particle positions, denoted by , , satisfy

at all times, i.e. the particles cannot occupy the same position or jump over each other. Next, define a small parameter

and microscale step size


The interparticle forces


are defined by a finite range potential . We suppose that for all within the range. Note that in (4.2) can take only two values: or . Also, observe , as it should be by the third law of Newton, and also that the sign of is the same as sign of . This means that the force exerted on by say, is repulsive. The total interaction force acting on the particle is

for .

Each particle has mass , where is the total mass of the system. Particles have velocities denoted by , . Writing the second Newton’s law as a system of first order equations yields the scaled microscale ODE system


subject to the initial conditions


5. Integral approximation of stresses for particle chains. Mesoscopic continuum equations

In the one-dimensional case stress is a scalar quantity, and (2.10), (2.11) reduce to, respectively,




The sum in (5.2) is simplified compared to the general expression, since we have exactly interacting pairs of particles.

To obtain integral approximations of stresses, we define interpolants , as in Sect. 3.2. Assuming as before that is invertible and repeating the calculations we get


Remark. Many equalities in the paper, including (5.3) hold up to a discretization error. To simplify presentation, we do not mention this in the sequel when discrete sums are approximated by integrals.

The interaction stress can be rewritten as


Next we approximate by . This approximation is in fact exact, provided the interpolant is chosen to be piecewise linear. Note also that

Inserting this into (5.4), replacing Riemann sum with an integral and changing variable of integration as in Sect. 3.2, we obtain the integral approximation of the interaction stress:


Equations (5.3), (5.5) contain two microscale quantities: and . Approximating sums in the definitions of the primary averages (2.6), (2.7) by integrals we see that and are obtained by applying the convolution operator to, respectively and :


The discretization error in (5.6) can be made small by imposing suitable requirements on the microscopic interpolants. Fortunately, the theory of ill-posed problems allows for errors in the right hand side of integral equations. The size of the error determines the choice of regularization parameter. In the present case, the error determines the number of iterations needed for the optimal reconstruction, according to so-called stopping criteria. These criteria are available in the literature on ill-posed porblems (see e.g. [9]). Detailed investigation of these questions is left to future work.

Denote by the iterative Landweber regularizing operators

Applying in (5.6) yields a sequence of approximations


and a corresponding sequence of closed form mesoscopic continuum equations (written here for an isolated system with zero exterior forces)


where are given by


with given by (5.7).

6. Zero-order closure for particle chains

Let us consider zero-order approximations in detail. The mesoscopic mesh consists of points

where , presumed to be an integer satisfying . Meso-cells are intervals of length

centered at .

Suppose that the only primary variables of interest are density and linear momentum . These variables will be computed by the mesoscale solver. For simplicity, suppose that the meso-solver is explicit in time. Then the average density and average velocity will be available at the previous moment of time. Our task is to design an update step for computing density and velocity at the next time moment. To construct a closed form update step, we need to approximate stress in (2.9) in terms of . From the knowledge of we can approximately recover and . The zero-order approximation (3.8) corresponds to


In other words, the microscale quantities are approximated by their averages. The corresponding closed form approximations for stress are obtained by inserting (6.1), (6.2) into (5.10), (5.11):


For computation, a numerical quadrature should be used. In this regard, note that all average quantities are computed on the mesoscale mesh, while the formulas (2.10), (2.11) are fine scale discretizations. Therefore, one might wonder if a straightforward mesoscale quadrature of (6.3), (6.4) is too crude. A better approach is to interpolate by prescribing approximate particle positions and velocities , , compatible with the given . Once this is done, (6.3), (6.4) can be discretized on a fine scale mesh with mesh nodes .

Interpolants cannot be unique. For zero-order closure, we are choosing positions and velocities that produce the given average density and average velocity. Clearly, there are many different position-velocity configurations with the same averages. The choice made in the paper is motivated by the practical requirement of achieving low operation count, as well as by certain expectations about the nature of dynamics. From the continuum mechanical point of view, if a system can be adequately modeled by balance equations of mass and momentum, then it must have have a trivial energy balance. Most often this means that the deformation is nearly isothermal. To mimic such isothermal dynamics we suppose that at each time step, there exists a positive number (it can be called ”upscaling temperature”) such that


Here the summation is over all particles located in a meso-cell . The temperature is the same for all . We emphasize that the actual value of is not as important as the fact that its value is the same for all meso-cells. This is because (6.5) would yield constant mesoscale mesh node values of As a result, the finite difference approximation of on the mesoscale mesh is identically zero. We interpret this by saying that convective stress does not contribute to the mesoscopic dynamics in the isothermal case. Another observation is that need not be the same at different moments of time, so our assumption is somewhat more flexible than the standard isothermal deformation approximation. Also, we note that validity (6.5) depends on the choice of . For bigger , it is more likely that (6.5) holds for the same underlying microscopic dynamics. Details on this are provided below in Section 6.2. Additionally, other features of the microscopic dynamics should be taken into account. Most importantly, interpolated velocities must be such that the collection conserves microscopic energy :


where is the microscale potential energy corresponding to the positions .

6.1. Prescribing particle positions

The objective of this section is to assign approximate particle positions . We start by interpolating . The simplest interpolant is piecewise-constant: where is the characteristic function of the meso-cell . A simple choice of the position map compatible with this interpolant is a piecewise linear map having the prescribed constant value of in each meso-cell. In practical terms, this means that in each meso-cell, particles are spaced at equal intervals from each other. The local interparticle spacing


is determined by the mesh value of the average density. To explain (6.7), note that the total mass of particles contained in the meso-cell can be approximated by . Dividing by the mass of one particle, we obtain an approximate number of particles inside :

and thus . We emphasize that are chosen based only on the known mesh values of the density , and that will be different from the actual particle positions .

Now we approximate the integral in (6.4) by its Riemann sum generated by the partition :


6.2. Prescribing particle velocities

In order to approximate the convective stress in the fine scale discretization of (6.3), we need to choose approximations of the true particle velocities . The choice of must satisfy (6.6) and be compatible with the available average velocity at the mesoscale mesh nodes.

For each , we set

where is the local average velocity, and is a perturbation to be defined.

Next, we show that the energy-conserving collection of velocity always exists, provided its upscaling temperature is suitably prescribed. This prescription will be based only on the available mesoscale information. For definitiveness, in the rest of this section we suppose that satisfies the following condition:


To make the algebra simpler, we make another assumption: for each ,


where is either or . The second summation is over all such that . Assumption (6.10) holds when is small outside of .

Averaging of should produce the known average velocity . This yields an equation for :



(6.11) holds provided


Now we look for perturbations in the form


where the are to be determined. Next, narrow down the choice of by setting


where is either one or negative one. To satisfy (6.12) we need to be even (one more point can be easily inserted if the actual is odd); in addition, the number of positive and negative must be the same.

To simplify further calculations, we write conservation of energy (6.6) in the form


where are any numbers satisfying


Our goal now is to show that there is a choice of , and such that defined by (6.13), (6.14) satisfy equations (6.5), (6.15), and (6.16). Inserting (6.13) into (6.5) and (6.15) yields, respectively,


Combining these equations we get


Substituting into (6.16) yields the choice of :


The choice of all constants now should be made as follows:
1) Given , find by (6.21);
2) Determine from (6.20);
3) Determine from (6.19);
4) Choose by (6.13), (6.14).
Note that step 4 introduces non-uniqueness, but we are concerned only with existence of suitable velocity perturbations. The actual choice of will not change the mesoscopic discretization of momentum balance equation. Indeed, once are chosen, we can approximate the integral in (6.3) (for at the mesoscale mesh nodes) by a Riemann sum corresponding to the partition of :

Therefore, the mesoscale mesh values of are all equal. This implies that a finite difference approximation of on the mesoscale mesh vanishes. The conclusion is that for isothermal dynamics, any suitable choice of a velocity perturbation produces a convective stress that has zero divergence on the mesoscale.

6.3. Zero-order isothermal continuum model

Combining the approximation with (6.3), (6.4) we obtain an isothermal zero-order continuum model


where is given by an integral expression (6.4) (or by a discretization (6.8)). We can interpret interaction stress as pressure. Then (6.4) provides dependence of pressure on density, which is non-local in space and non-linear. For small (large ), it can be approximated by

which is still non-local. In the limiting case , observing that convolution with is an approximate identity, this equation can be reduced a local equation of state

If is nearly constant, this equation can be linearized to produce a classical gas dynamics linear equation of state. This shows that zero-order closure (6.4) generalizes several classical phenomenological equations of state. The connection between micro- and mesoscales is made explicit in (6.4). Using higher order closure approximations, one can obtain other non-classical continuum models worth further investigation.

7. Computational results

In this section, the method developed in the previous sections is tested for a chain of to particles interacting with a non-linear finite rage potential