其仿真过程的核心算法参考了 Robert Bridson 的 Fluid Simulation for Computer Graphics。它使用 Navier-Stokes 方程的数值解来预测 Qt 的 OpenGL 小部件中显示的每一帧中粒子的密度和速度分布。有两种分辨率可供用户选择：64 × 64 和 128 × 128。用户还可以在 7 种颜色的密度和速度之间切换显示模式。
Its core algorithm of the simulation process refers to Robert Bridson's Fluid Simulation for Computer Graphics. It uses the numerical solution of Navier-Stokes equations to predict the density and velocity distribution of particles in each frame displayed in the OpenGL widget of Qt. There are two resolution options for users: 64 × 64 and 128 × 128. Users can also change display mode between density and velocity in 7 colors.
 |
 |
|---|---|
| 密度场 (Density field) | 速度场 (Velocity field) |
git clone https://github.com/MuGemSt/FluidS.git
cd FluidS此部分内容请统一参考此博文
Please refer to this blog post for this section.
The core problems solved are as follows:
Since the equations to be processed are made in two lines in real-time, one is the velocity field, and the other is the density field. And what we're going to end up visualizing is the latter. However, real-time updates of the former are also essential because the latter's real-time updates need to be supported by the former.
The solution flow for the speed field and the density field is as follows:
First, the source force item is handled, and the processing of it is also in processing concentration supplement of the density field. If we assume that the external force added throughout the process does not change dramatically,
It is an excellent approximation within the Δt time interval and is highly compatible with interactive systems. The effect of force occurs only at the beginning of each moment.
Next, the convection terms are calculated. Since the velocity field and density fields are solved similarly, the relatively complex velocity field is directly explained. It is precise because of the existence of −(u · ∇)u that the Navier-Stokes equation becomes a nonlinear equation. Foster and Metaxas use the finite difference method when processing this term. Their algorithm can only guarantee stability when the time step is tiny, and the condition is ∆t < ∆τ/|u|, where delta t is the scale of the current computing grid. Here uses Stam's semi-Lagrangian characteristic line solution method. Convection item questions can be abstracted into the following form:
, where a is a scalar field, v is a stable vector field, and a0 is the first field at time t = 0. Use P(x0, t) to represent the characteristic line of the vector field v at time t = 0 passing through the point x0:
Now let ̅a (x0, t) = a(P(x0, t)) be the value of the field along the characteristic line passing through the point x0 at time t = 0. The following difference law chain can calculate the rate of change of this quantity with time:
It indicates that the value of this scalar is not consistent along the curve. In particular, we have ̅a (x0, t) = ̅a (x0, 0) = a0(x0). Therefore, the first field and the characteristic line can completely define the solution to the convection problem. The field at a given time t and position x can be calculated by x back to the previous position along the characteristic line in real-time to get the point x0, and then approximate the first field at that point:
We use this method to solve the fluid's convection equation in the time interval [t, t + ∆t]. Where v = u(x, t) and a 0 is the velocity component of any fluid at time t. Use the above method to advance u1 to u2 as follows:
Next, the viscous term is calculated. Since the relationship between the two-dimensional viscous term and the one-dimensional diffusion term is merely one-dimensional, the solution of the viscous item is the solution to each component's diffuse issue. In the previous introduction, it was also mentioned that it could be processed by Gauss-Seidel iteration, namely:
, where the I is the unit operator. When the diffusion operator is discretized, solving this format requires solving a sparse linear system to get u3. The last thing that needs to be done is to introduce a projection operator to make the obtained field a passive field so that u3 can be used for u0 at the next moment. To solve the above problem, we need to start with the Helmholtz-Hodge decomposition. For any vector field w can be decomposed into the following form:
, where u is a passive field, and q is a scalar field. According to this conclusion, we define a projection operator P so that any vector field w can be projected to the part with a divergence of 0, that is, u = Pw. The construction of this operator can be obtained by multiplying the gradient operator on both sides of the above formula at the same time:
At this time, a Poisson equation about the scalar field q is obtained with Neumann boundary conditions ∂q/∂n = 0 on ∂D. The process of solving is the process of calculating u projection:
Applying this projection operator to both sides of the original Navier-Stokes equation gives a simple comparison of the velocity field:
According to the above description, the projection process for this step is transformed into a problem to solve the Poisson equation, namely:
After the Poisson problem's discretization, a sparse linear system is obtained with a solution complexity of O(N). Since the entire process is calculated in the ∆t period, the algorithm's operation efficiency is also critical. For the operation processing of the operators above, we use the FFT algorithm. Fourier transform is often used to deal with periodic boundary conditions. The problem, our current challenge, is a particular simple situation. Periodicity allows us to map the velocity field to Fourier space:
The gradient operator ∇ in Fourier space can be equivalent to multiplying by ik, where i^2 = −1. In this way, the calculation of the previous diffusion operator and projection operator will become very simple:
, where k = |k|. The operator ̂P maps the vector ̂w(k) to a relatively regular region for wavelet k. The Fourier transform of the velocity in the passive field therefore, always remains orthogonal to the wavelet. The diffusion operator can be regarded as a filter whose attenuation is proportional to the time step and the viscosity coefficient. Therefore, we can completely convert the solution operator into two routes. Both are optimized with the following O(NlogN) complexity FFT algorithm:
The last task is to embed the adaptive grid algorithm into each solution plate. As mentioned earlier, the density gradient is used as a criterion. If the density gradient at the center of a block exceeds the threshold, the mesh subdivides and solves the problem. According to the density gradient level, you can determine how many degrees the area needs to be divided so that it is determined to iterate on until the highest density gradient is completed.