# Multiphase flow modelling – 2: Lagrangian particle tracking

#### Parameters

• Maximum number of particles: npmax
• Particle diameter: dp
• Particle density: rhop
• Fluid density: rhof
• Particles flow rate: PFR (number of particles entering per second)
• Particle inlet velocity: up_in
• Two-way coupling(Boolean): use_twoway
##### Define a derived data type for particle with following detail
• double x, y : position
• int i,j : mesh location
• double u,v : velocity
• double dt : Time step size
• int or boolean stop: active or not
##### Create an array of particles part, ptmp

Identify the current mesh location(i,j) of the particle based on particle position

This is rather easy in a structured grid compared to an unstructured one. The idea is particle would be somewhere near to the location it was at the previous time step. Hence we shall search the immediate four(in 2D) neighbours and if not found, search neighbours’ neighbours.

#### Algorithm

• grab the old location: iold, jold
• grab the current position: xp, yp
• define now location: inew, jnew
• Search to the right: do while (xp >=x(inew+1) && inew+1 <= nx); inew=inew+1;
• Search to the left: do while (xp <= x(inew) && inew-1>= 1 ); inew=inew-1;
• Either of the above search will only execute.
• Do the same for yp

#### Bilinear interpolation

The new location we arrive is marked in the figure below. Linearly interpolate for U & V.

U_{P} = \frac{U1*Green+U2*Purple+U3*Yellow+U4*Red}{Green+Purple+Yellow+Red}
wx1=(x0-x (i))/(x (i+1)-x (i)); wx2=1.0_WP-wx1
wy1=(y0-ym(j))/(ym(j+1)-ym(j)); wy2=1.0_WP-wy1
// Rectangular area in a matrix
ww(1,1)=wx1*wy1
ww(2,1)=wx2*wy1
ww(1,2)=wx1*wy2
ww(2,2)=wx2*wy2
//Bilinear interpolation of U
myu=sum(ww*U(i:i+1,j:j+1))

Particle initialisation

• give x position that of inlet
• y position random between y(1) and y(ny)
• locate i,j using previously defined plocate algorithm
• If any particle are in wall (mask ==1), remove that particle (nadded++, np++)
• Initialise u velocity to the particle

#### ODE solver for a single particle: Second order Runge-Kutta

Step 1: Predictor (Forward Euler) : \phi^{n+\frac{1}{2}} = \phi^{n} + \frac{\Delta t}{2}f(\phi^{n})

Step 2: Corrector (Centered in time) : \phi^{n+1} = \phi^{n} +\Delta tf(\phi^{n+\frac{1}{2}})

Algorithm

for each particle p

Predictor – Euler prediction

• copy p to pn
• Calculate drag(dragx and dragy) by calling lpt_drag(dragx,dragy).
• p.x = pn.x +0.5*p.dt *p.u
• p.y = pn.y +0.5*p.dt *p.v
• p.u = pn.u +0.5*p.dt *(dragx+gravity_x+collision_x)
• p.v = pn.v +0.5*p.dt *(dragy+gravity_y+collision_y)

Corrector – Midpoint rule

• Recalculate drag(dragx and dragy) by calling lpt_drag(dragx,dragy).
• p.x = pn.x +p.dt *p.u
• p.y = pn.y +p.dt *p.v
• p.u = pn.u +p.dt *(dragx+gravity_x+collision_x)
• p.v = pn.v +p.dt *(dragy+gravity_y+collision_y)

call plocate to relocalize particle on mesh

Update optimal timestep size p.dt=taup

#### Drag Model (lpt_drag(dragx,dragy).)

Step 1: Obtain particle velocity by interpolation. Input: i,j,x,y Output: u,v

ug=get_u(p.i,p.j,p.x,p.y) and vg=get_v(p.i,p.j,p.x,p.y)

Step 2: Calculate Stokes response time

taup=dp^2/(18.0*knu)

Step 3: Calculate particle Reynolds number

Re_{p} = \sqrt{[p.u - u_{g}]^{2} +[p.v - v_{g}]^{2} } . \frac{dp}{knu}

Step 4: Shiller & Naumann drag correlation

taup=taup/(1.0+0.15*Rep*(0.687))

Step 5: Acceleration due to drag
accx=(ug – p.u)/taup
accy=(vg – p.v)/taup

accx, accy are dragx and dragy

#### two way coupling

For two-way coupling, add momentum source term calculated from particle to the u* calculation in the velocity step.

 // Increment exchange term
end if

// Transfer source to flow solver
if (use_twoway) then
do j=1,ny
do i=2,nx
end do
end do
do j=2,ny
do i=1,nx

//In velocity_step
V=V+HV1+ABcoeff*(HV1-HV2)+srcV