The Exact Solution of All of Continuum Mechanics

By Richard D. James, Alessia Nota, Gunjan Pahlani, and Juan J. L. Velázquez

Let \(A\) be any \(3 \times 3\) matrix, assume \(t\) is small enough that \(\det(I + tA)>0\), and consider the function

\[v(x, t) = A(I + tA)^{-1} x, \quad x \in \mathbb{R}^3.\tag1\]

Make one of the many choices of \(A\) such that \(\textrm{div} \, v = 0\) and substitute into the Navier-Stokes equations of fluid mechanics:

\[\frac{\partial v}{\partial t} +  (v \cdot \nabla)v  = -\nabla p + \nu \Delta v.\tag2\]

Choosing for example \(p = p(t)\), a quick calculation shows that \(v(x,t)\) is an exact solution. In fact, both sides of the Navier-Stokes equations vanish. Equation \((2)\) is for the incompressible case, but \(v(x,t)\) also solves the Navier-Stokes equations for a compressible fluid; for the latter, choose any \(A\) and the corresponding density given by \(\rho(x,t) = \rho_0 \det (I + tA)^{-1}\), \(\rho_0 = const. > 0\). 

Now try the equations of elasticity (linear or nonlinear). Again, they yield an exact solution for every \(A\). To witness this phenomenon in the general nonlinear case, we must first convert from the Eulerian to the Lagrangian description of motion by solving the system of ordinary differential equations (ODEs) \(\dot{y}(z,t) = v(y(z,t),t)\), \(y(z,0) = z\); the solution is \(y(z,t) = (I + tA)z\). We then immediately see that the following holds for every elastic material:

\[\rho_0 \frac{\partial^2 y}{\partial t^2} = {\rm div} \,T, \quad T = \tilde{T} (\nabla y).\tag3\]

At first glance, this is quite surprising. Aren’t solutions of partial differential equations supposed to depend on the choice of coefficients? Why should a particular solution for air also work for water, rubber, and steel? Very different physics—and very different kinds of atomic forces—underly the flow of a gas or liquid when compared to the forces that bind a crystalline solid. Why is the exact same nine-parameter family of functions possible in both scenarios? The motions that \(v(x,t)\) describe are not trivial; they may be volume preserving, they can have vorticity that changes in time, they have a strong singularity when \(\det(I + tA) \to 0\), and they can most certainly exist far from equilibrium.

How about something more exotic? Consider any one of the accepted but highly nonlinear models of non-Newtonian fluids, the theory of liquid crystals, or nonlinear viscoelasticity. Or think about one of the models of phase transformations with a change of type (regularized or not). These examples again produce an exact solution for every \(A\). 

What about thermodynamics? Here we finally find a departure from universality. In all of the aforementioned cases—augmented by the first and second laws of thermodynamics—\(v(x,t)\) still satisfies the laws of conservation of mass and momentum, with temperature \(\theta(t)\) assumed to be a function of time only. But the energy equation becomes an ODE for \(\theta(t)\). The coefficients of this ODE depend on the material, and the corresponding evolution of temperature is material dependent.

Symmetry often accounts for coincidences like this, as is the case here. So, where is the symmetry group?

It is best to descend to atomistic theory in order to answer this question. Let us consider the molecular dynamics (MD) of atoms with positions \(y_k(t)\) and masses \(m_k\). These satisfy the equations of MD:

\[m_k \ddot{y}_k(t) = -\frac{\partial \varphi}{\partial y_k} (\dots,  y_{i-1}(t), y_i(t), y_{i+1}(t),  \dots).\tag4\]

The well-known symmetries of quantum mechanics—the frame indifference and permutation invariance of atomic forces—are here. We wish to examine cases wherein infinitely many atoms fill space, analogous to the function \(v(x,t), \, x \in \mathbb{R}^3\). The potential energy \(\varphi\) is typically infinite in these cases, but the infinite system of ODEs of MD can still make perfect sense. We would therefore like to express these symmetries in terms of the forces \(-{\partial \varphi}/{\partial y_k}\). The fundamental symmetries are

\[\frac{\partial \varphi}{\partial y_k} (\dots, Qy_{i-1} + c, Qy_i + c, Qy_{i+1} + c, \dots) \,  = \, Q \frac{\partial \varphi}{\partial y_k} (\dots, y_{i-1} , y_i , y_{i+1} , \dots) \,\, \textrm{(frame indifference)}\]

\[\frac{\partial \varphi}{\partial y_k} (\dots,  y_{\Pi(i-1)}(t), y_{\Pi(i)}(t), y_{\Pi(i+1)}(t),  \dots) \, = \, \frac{\partial \varphi}{\partial y_{\Pi(k)}} (\dots, y_{i-1} , y_i , y_{i+1} , \dots)  \, \, \textrm {(permutation invariance)},\]

where \(Q \in \textrm{O}(3)\), \(c \in \mathbb{R}^3\), and \(\Pi(\cdot)\) is a permutation that preserves the species of atom.

Both \(\textrm{O}(3)\) and the translation group are continuous groups, but we have discrete atomic positions. This fact suggests that we should examine discrete groups of isometries, i.e., elements \((Q|c)\)—often written in this notation—with \(Q \in \textrm{O}(3)\), \(c \in \mathbb{R}^3\). These can generate many groups of the form \(\mathcal{G} = \{ (Q_1|c_1), (Q_2|c_2),\) \(\dots, (Q_N|c_N)\}\), some of which are listed in the International Tables for Crystallography. \(N\) can be infinite, and the group product is \((Q_1|c_1) (Q_2|c_2) = (Q_1Q_2 | c_1 + Q_1 c_2)\). One last crucial assumption is that we will allow \(c_1, c_2, \dots\) to depend on time, but we only allow affine time dependence: \(c_1 = a_1 t + b_1\), \(c_2 = a_2 t + b_2\), \(c_3 = a_3 t + b_3, \dots\). The \(a_i\) and \(b_i\) must be chosen so that \(\mathcal{G}\) remains a group for, say, all \(t>0\). Although this vaguely resembles a Galilean transformation, it is not because the \(a_i, b_i\) vary between group elements.

Consider any (discrete) isometry group \(\mathcal{G}\)—possibly time dependent as described—and designate a subset of atoms with positions \(y_1(t), \dots y_M(t)\) as the simulated atoms, with initial positions \(y_1^{\circ}, \dots, y_M^{\circ}\) and initial velocities \(v_1^{\circ}, \dots, v_M^{\circ}\). We obtain the other atom positions by applying \(\mathcal{G}\) to the simulated atoms using the obvious rule: \(Q_1 y_1(t) + c_1, Q_1y_2(t) + c_1, \dots,  Q_2 y_1(t) + c_2, Q_2y_2(t) + c_2, \dots, Q_N y_1(t) + c_N, Q_Ny_2(t) + c_N,\dots\). Now we write the MD equations from above, but only for \(M\) simulated atoms. On the right-hand side of \((4)\), substitute these formulas for all other nonsimulated atoms. Because the nonsimulated atoms are given by formulas in terms of the simulated atoms, the equations for the simulated atoms become a system of (nonatonomous) ODEs in standard form. The necessary initial conditions are \(y_1^{\circ}, \dots, y_M^{\circ}\), \(v_1^{\circ}, \dots, v_M^{\circ}\). We then solve these equations. A straightforward theorem reveals that even though the nonsimulated atoms are given by formulas, these atoms exactly satisfy the MD equations for their forces; we utilize both frame indifference and permutation invariance to show this. Animation 1 depicts a simulation.

Now we are ready to explain \(v(x,t)\). We select the simplest time-dependent translation group, with elements that can be written as \((I \, | \, (I + t A)  (\nu^1 e_1 + \nu^2 e_2 + \nu^3 e_3))\). Here, \(\nu^1, \nu^2, \nu^3\) are integers and \(e_1, e_2, e_3\) are linearly independent. This is clearly a group, and it has the affine time dependence. We can also determine the macroscopic velocity field of this solution with minimal difficulty: it is exactly \(v(x,t)\). A clear victory for continuum mechanics!

By examining this argument more closely, we can see that it really pertains to a largely unstudied, explicit, and time-dependent invariant manifold of the MD equations. In fact, lots of groups, time dependencies, and these types of invariant manifolds are present. Their forms are independent of the material; the same manifolds apply to air, water, and steel.

We now arrive at the article’s main purpose. Recall that the molecular density function \(f(t, x, v)\) of the kinetic theory of gases represents the probability density of finding an atom with velocity \(v\) in a small neighborhood of \(x\) at time \(t\), in Eulerian form. Let us stay with the time-dependent translation group and examine the statistics of the MD solutions. Draw a ball \(\mathcal{B}_0\) of any diameter that is centered at the origin. Now choose integers \(\nu^1, \nu^2, \nu^3\) and draw a ball \(\mathcal{B}_\nu\) of the same diameter that is centered at \(x = (I + t A)  (\nu^1 e_1 + \nu^2 e_2 + \nu^3 e_3)\). Since the simulated atoms quickly diffuse into the nonsimulated atoms during a simulation (see Animation 1), it is not unusual for \(\mathcal{B}_0\) and \(\mathcal{B}_\nu\) to contain both simulated and nonsimulated atoms at any particular time. The velocities of atoms in \(\mathcal{B}_0\) and \(\mathcal{B}_\nu\) are different. Nevertheless, if we know the velocities of atoms in \(\mathcal{B}_0\), we can then use explicit formulas to calculate the velocities of atoms in \(\mathcal{B}_\nu\). But based on its interpretation, this calculation implies an ansatz for the molecular density function \(f:\) \(f(t, x, v) = g(t, v - A(I + tA)^{-1} x)\). Ignore the fact that \(x\) was a special point and substitute the ansatz into the general form of the Boltzmann equation. It yields an exact reduction.

Notice that we lose the periodicity regardless of whether we pass from MD to the reduced Boltzmann equation or to continuum mechanics. This is comforting.

We have been studying this reduced equation for \(g\) and simulating the invariant manifold of MD itself. We do not believe that anything like a statistical description of nonequilibrium exists in general. After all, the Boltzmann equation—which treats only the simplest kind of material—gives an initial value problem for a nonlinear integro-differential equation in seven variables. There is no chance that a simple reduction to an explicit form of \(f\) exists far from equilibrium. However, on the manifold given by \(g\)—or perhaps on one of the invariant manifolds of the aforementioned general equations of MD—a kind of explicit non-equilibrium statistical mechanics might exist that approaches the simplicity of equilibrium statistical mechanics, which also lies on an invariant manifold \(H = const\). Finding this, even in special cases, would be a tremendous achievement. Our work on the Boltzmann equation perhaps gives some hints, but more ideas are necessary.

Richard D. James teaches continuum mechanics at the University of Minnesota. Comments can be addressed to him at [email protected]. Alessia Nota is a tenure-track assistant professor at the University of L’Aquila (DISIM). She was previously a postdoctoral researcher at the University of Bonn and the University of Helsinki. Her research area is non-equilibrium statistical mechanics, with a particular focus on the kinetic theory of gases and plasmas. Gunjan Pahlani is a Ph.D. student in the Department of Aerospace Engineering and Mechanics at the University of Minnesota. Her research focuses on the formulation of computational methods at the molecular level that exploit long-range symmetries for the simulation of hypersonic gas flows and plastic deformation in materials. Juan J. L. Velázquez is a professor in the Institute for Applied Mathematics at the University of Bonn, prior to which he was a professor in the Department of Applied Mathematics at the Complutense University of Madrid and a postdoctoral Fulbright Fellow at the University of Minnesota. His research focuses on several problems of partial differential equations, particularly in the study of kinetic equations.