How can I convert a Brusselator type model of reaction-diffusion to include vortical motion? A system called Model G, which is similar to the Brusselator, describes reaction-diffusion processes taking place among hypothetical sub-quantum units. We are investigating whether it would be possible to add vortical motion to this system?
The proposed sub-quantum units differ from quarks in that they exhibit only standard classical properties such as reaction and diffusion. As a result of their processes, quantum structures are found to emerge which have properties and behaviors similar to those exhibited by subatomic
particles.
This would be to mathematically describe the reaction-diffusion system in terms of fluid dynamics which includes vortices. One idea is to include Navier-Stokes equations or something similar for the convection process, i.e. reaction-diffusion-convection (non-linear) or vortical convection of the hypothetical sub-particle reactants. How would one go about this?
The equations for the Brusselator type system (Model G) in their present form, without vortical motion, are shown below. The figure shows the Brusselator type Model G, as a schematic. The 2nd image shows the five kinetic equations, followed by their PDEs.


\begin{equation}
\frac{\partial G(x,y,z,t)}{\partial t}=D_G \nabla^2 G - (k_{-1} + k_2)G + k_{-2} X + k_1 A \quad \textrm{(4a)}
\end{equation}
\begin{equation}
\frac{\partial X(x,y,z,t)}{\partial t}=D_X \nabla^2 X + k_2 G - (k_{-2} + k_3 B + k_5) X + k_{-3}ZY - k_{-4}X^3 + k_4 X^2Y +k_{-5} \Omega \quad \textrm{(4b)}
\end{equation}
\begin{equation}
\frac{\partial Y(x,y,z,t)}{\partial t}=D_Y \nabla^2 Y + k_3 B X - k_{-3} Z Y + k_{-4} X^3 - k_4 X^2 Y \quad \textrm{(4c)}
\end{equation}
Letter symbols A and B denote the concentrations of the initial particle reactants; G, X, and Y denote the concentrations of the intermediate reactants; and Z and Ω denote the concentrations of the final reaction products. $D_i$ is a constant of proportionality called the diffusion coefficient. 
The forward reaction rate constants, $k_{i}$ , above each arrow, and the reverse reaction rate constants, $k_{-i}$ , below each arrow specify the rate at which reactants (back of arrow) transform into products (front of arrow).
 A: Here is a quote from NASA Navier-Stokes site. It mentions CFD as a way to solve approximations to Navier-Stokes, as these equations are too complex to solve analytically.
https://www.grc.nasa.gov/WWW/k-12/airplane/nseqs.html
"The equations are a set of coupled differential equations and could, in theory, be solved for a given flow problem by using methods from calculus. But, in practice, these equations are too difficult to solve analytically. In the past, engineers made further approximations and simplifications to the equation set until they had a group of equations that they could solve. Recently, high speed computers have been used to solve approximations to the equations using a variety of techniques like finite difference, finite volume, finite element, and spectral methods. This area of study is called Computational Fluid Dynamics or CFD".
One could use finite element method to solve the Model G reaction-diffusion equations for X, Y and G concentrations. And additonally an approximation to Navier-Stokes could 
be run separately within finite elements to calculate the vortical motion. Therefore having 2 simulations running along side each other, evolving together and complementing each other.
In the Navier-Stokes solution, one could start with the supposition that the reaction process is the flux inlet or source and the diffusion is the outlet or sink. One could then test from the simulation whether this gave the correct outcome of subatomic particle spin, by comparing it with current particle data sets.
