# Numerically Solving Fluid Dynamics

I'm trying to model the behavior of a scalar field $$\rho(\vec x)$$. Associated with $$\rho$$ there is an energy density $$\mu$$ like so: $$\mu = \tau[\rho,Q] + u[\rho,Q]$$ $$u = u_{ext}[\rho,Q] + u_{ee}[\rho]$$ The symbols $$\tau$$, $$u$$, $$u_{ext}$$ and $$u_{ee}$$ represent scalar fields that will be computed from $$\rho$$ and the knowledge of a set of "point charges" $$Q$$ that don't change over time but which influence the dynamics of $$\rho$$. The terms between square brackets represent the dependency of each field in the equations with either $$\rho$$ or $$Q$$.

I want to study two types of equations of motion, the first one is a continuum version of Newton's second law of motion:

$$\partial_t \vec v = -\nabla \mu$$ The second type is $$\partial_t \rho = -\nabla \mu$$ I understand that these two equations describe different types of motion and one of the points of my study is to determine which one is more appropriate for my purposes, so ideally I would like to program both.

In order to do that I have to know how to compute $$\rho_{t+\delta t}$$ as a function of $$\rho_{t}$$ from either of the equations above. I read somewhere that in order to determine the dynamics of a fluid I have to solve a set of differential equations called the Navier-Stokes Equations.

I understand that to describe the dynamics of a set of point particles following classical mechanics, for example, formally I have to solve Newton's equations of motion. But when you actually sit down to write a software that does that, there are actually a number of algorithms you can use, such as Verlet or Leap-frog, that generate trajectories which are solutions to Newton's equations. I was wondering whether there was any such algorithm I could similarly implement to solve my hydrodynamical problem.

So my question is: What time-stepping algorithms are available to model the dynamics of a 3D field like this one?

• I'm a little confused, you start off mentioning simulating molecules, but then discuss "cubarray" and how you can do math with these things and ask how to use them. If you are interested in molecules, hydrodynamics probably isn't what you're looking for (you'd want molecular dynamics instead). If you're wanting hydrodynamics, then you just need to solve the NS/Eulerian system as a coupled set of PDEs. For references, see this post. Commented Oct 9, 2022 at 17:59
• I'm modeling the electron density of a molecular geometry as a scalar field. I want to know how to update that field as time goes on as I am not very used to working with the dynamics of volumetric data. Do you think it would help if I erased any mention of molecules from the post and left things more mathematical? Commented Oct 9, 2022 at 19:23
• I would cut it down to what your question is. You seem to be mixing a few things together that I'm not really sure what's being asked (best I can see is "What time-stepping algorithms apply for hydrodynamics?" but I don't know how that fits in with what you've said w.r.t. fields). Commented Oct 10, 2022 at 1:02
• @KyleKanos "What time-stepping algorithms apply for hydrodynamics?" seems like a really good idea of what I want. I'll edit the question. Commented Oct 10, 2022 at 10:02

While hydrodynamic equations are represented as a coupled set of PDEs, the time stepping is generally treated as an ODE instead, $$\frac{\mathrm d\mathbf u}{\mathrm d t}=F(\mathbf u(t))$$ where $$\mathbf u=(\rho,\,\pi_x,\,\pi_y,\,\pi_z,\,E)^T$$ is the state vector ($$\rho$$ the density, $$\pi_i=\rho v_i$$ the momentum density in the direction $$i$$ and $$E$$ the total energy) and $$F(\cdot)$$ the flux. Hence, any of the methods used in solving an ODE are applicable to hydrodynamics.
• Forward Euler method update $$\rho(t+\Delta t)=\rho(t)-\nabla\mu(t)\cdot\Delta t$$ which is ubiquitous due to its simplicity.
• Backward Euler method, $$\rho(t+\Delta t)=\rho(t)-\nabla\mu(t+\Delta t)\cdot\Delta t,$$ typically involves matrix operations/reductions and are thus more costly than the forward update, but is still an option.
• Predictor-corrector method, which solves the state at $$t+\Delta t$$ using a predicted value & an interpolated value, \begin{align} \hat{\rho}(t+\Delta t)&=\rho(t)-\nabla\mu(t,\,\rho(t))\cdot\Delta t \\ \rho(t+\Delta t)&=\rho(t)-\frac{1}{2}\left(\nabla\mu(t+\Delta t,\,\hat{\rho})+\nabla\mu(t,\,\rho)\right)\cdot\Delta t \end{align}
Note that, in hydrodynamics, your time step $$\Delta t$$ is constrained by the Courant-Friedrichs-Lewy condition, which essentially argues that the rate of information transfer (e.g., speed of sound) limits how far material can travel from one cell to another in one time step, lest unphysical negative densities/pressures or superluminal velocities occur. This conditions enforces that, $$\Delta t\leq \mathcal{C}\frac{\Delta x}{u_\text{max}}$$ where $$u_\text{max}$$ is the maximum characteristic speed over the whole domain at the current point and $$\mathcal{C}\leq1$$ is a constant for the whole simulation. One can derive such time step constraints by considering the von Neumann stability analysis.