Estimation of time taken to reach steady state in an MD simulation of Poiseuille flow

Question

I am trying to do a Molecular Dynamics Simulation of a complex fluid, confined between solid surfaces. I would like to find the flow rate as a function of fluid film thickness, $h$, for a plane poiseuille flow. (The flow rate is not expected to show an $h^3$ dependence, in general, if $h$ is small, in the order of a few nanometers).

The question I have is, how I can get an estimate for the simulation time required for the system to attain steady state, so the flow rate measurements can be made after such time. If we assumed continuum, we would have this equation:

$ \begin{align*} \frac{\partial M}{\partial t} = \nu \frac{\partial^2 M}{\partial y^2} - \frac{1}{\rho} \frac{\partial P}{ \partial x} \end{align*}$

to describe this unsteady Poiseuille flow($M$ being the Momentum, $\nu$ the kinematic viscosity and $P$ the pressure), with the applied pressure gradient being very high (since it is applied via adding high elementary body forces to the particles) to start with, when compared to the viscous diffusion term. The pressure gradient term can be assumed to be constant along the direction of flow $x$ and in time, $t$.

Now, if we didn't have the pressure term, we would have a momentum diffusion equation from where we could estimate the time required for the unsteadiness to decay. But, from the above equation (or from any other starting point), it is not obvious to me how one can get an estimate for the simulation time required for steady state to be reached.

Thank you very much for your patience and help!

Answer 1

I am not sure i understand your setup exactly, but why can't you just use the analytic expression of the laminar solution between parallel plates to get the flow rate? Of course assuming your Reynolds number is small enough so that the laminar solution is globally stable.

A different approach, which you can apply also to fluid flows in more complex geometries, could be to estimate the decay rate of the least stable eigenmode of the linear stability problem, (again assuming the system is linearly and globally stable) to get an estimate of the time it takes to reach a certain decay of the amplitude of the initial disturbance.

Answer 2

A good estimate can be obtained by looking at: $ \begin{align*} \frac{\partial M}{\partial t} = \nu \frac{\partial^2 M}{\partial y^2} \end{align*}$

which is a diffusion equation, that has a solution of the form: $M(t) = \frac{M_0}{\sqrt{4\pi \nu t}} \exp{(-\frac{y^2}{4 \nu t})}$. In theory, therefore, $h^2/\nu $, where $h$ is the fluid film thickness is a fair estimate for the time required to reach steady state. But, in practice, for complicated fluids, the MD simulation seems to be required to be run for much longer,for the fluctuations in the binned averages of the measured quantities to die down.