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

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!

• Are you thinking about analytical solutions for estimates? Is the pressure known or you need one more equation to determine? Aug 28 '15 at 2:45
• Yes, but I suppose the transient solution of the above doesn't depend on the pressure gradient. So, $h^2/\nu$ would still be a good estimate. Aug 29 '15 at 0:57

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.