Greetings! In the previous blog, I mentioned that I applied to Parker-Sochacki to another part of the FLIP fluid simulation: the integrate forces section of the pop network. Inside the integrate forces section, is an Ordinary Differential Equation that handles air resistance serving as a type of wind drag. This updates the target velocity.
This focuses on how the wind drag was updated so there is no simulation attached.
First go to the Particle Forces Section of the Flip Solver.
Second dive into the integrate forces node
Third go inside the implicit_targetv
That is where the following sticky note is that talks about the Ordinary Differential Equation:
————————————————————-
Our ODE is
dv/dt = (- a |v-vt| – b |v-vt|^2) ((v-vt)/|v-vt|)
First reframe into the local wind speed, so v-vt becomes just v.
Multiplying through we get
dv/dt = -av – b|v|v
Ignoring all other forces active over the timestep, we can scalarize this by only considering change along v, so with v now the scalar speed,
dv/dt = -av -bv^2
We want to solve for the initial condtions v0 = incoming velocity at time timestep.
So, ignoring initial conditions, we have something of the order
v = a / (b * (exp(a * t) – 1)) [A]
Note in the case of a == 0, it is 1/tb, and in the case of b == 0, exp(-at).
So, solving for the t0 given a v0,
t0 = ln(a/(b*v0)+1) / a
if a == 0, it is
t0 = 1 / v0*b
if b == 0 it is
t0 = -ln(v0)/a
Our new speed is thus evaluating [A] at t0 + timestep
We scale (v-vt) by this ratio and add back in vt to restore our original frame.
Note we can compute the ratio directly and get a fair bit of cancellation to simplify the equations.
————————————————————-
It has an analyical solution so there is no advantage to solving it with Parker-Sochacki. This was just Parker-Sochacki/Houdini practice. We successfully implemented the following: ODE:
dv/dt = -av -bv^2
where a and b are both constants for air resistance.
This could be done with the following iteration:
I originally did this in a more complicated way where I substituted c = v^2 and had another iteration recurrence to update c but that isn’t necessary.
Here is a snippet of code that was implemented. This was done in HScript, which is very similar to C++ for this example. I set a to 1 and b to 2. I gave velocity an initial value of 2 for x,y,and z components.
————————————————————-
//PS code
int seriesSize = 70;
vector a = { 1, 1, 1 };
vector b = { 2, 2, 2 };
v = { 2, 2, 2};
float a_x = a.x;
float a_y = a.y;
float a_z = a.z;
float b_x = b.x;
float b_y = b.y;
float b_z = b.z;
float v_x[];
float v_y[];
float v_z[];
resize(v_x, seriesSize);
resize(v_y, seriesSize);
resize(v_z, seriesSize);
v_x[0] = v.x;
v_y[0] = v.y;
v_z[0] = v.z;
for (int j = 0; j < i; j++)
{
sum_x += v_x[j]*v_x[i-j-1];
sum_y += v_y[j]*v_y[i-j-1];
sum_z += v_z[j]*v_z[i-j-1];
}
for (int i = 1; i <= seriesSize; i++)
{
v_x[i] = (-a_x*v_x[i-1] – b_x*sum_x) / i;
v_y[i] = (-a_y*v_y[i-1] – b_y*sum_y) / i;
v_z[i] = (-a_z*v_z[i-1] – b_z*sum_z) / i;
float vsum_x = 0;
float vsum_y = 0;
float vsum_z = 0;
float tinctemp = .014667;
for (int seriesindex = 0; seriesindex <= seriesSize; seriesindex++)
{
vsum_x += v_x[seriesindex]*pow(tinctemp, seriesindex);
vsum_y += v_y[seriesindex]*pow(tinctemp, seriesindex);
vsum_z += v_z[seriesindex]*pow(tinctemp, seriesindex);
}
v.x = vsum_x;
v.y = vsum_y;
v.z = vsum_z;
adamwermus 