diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index 0248e4375d..583f3ed303 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -106,9 +106,6 @@ void Foam::KinematicParcel::calc vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx, Cud, dUTrans); - // Constrain the new velocity for reduced -D cases - meshTools::constrainDirection(mesh, mesh.solutionD(), U1); - // Accumulate carrier phase source terms // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -144,6 +141,8 @@ const Foam::vector Foam::KinematicParcel::calcVelocity vector& dUTrans ) const { + const polyMesh& mesh = this->cloud().pMesh(); + // Return linearised term from drag model Cud = td.cloud().drag().Cu(U - Uc_, d, rhoc_, rho, muc_); @@ -160,6 +159,9 @@ const Foam::vector Foam::KinematicParcel::calcVelocity vector Unew = td.cloud().UIntegrator().integrate(U, dt, ap, bp).value(); + // Apply correction to velocity for reduced-D cases + meshTools::constrainDirection(mesh, mesh.solutionD(), Unew); + // Calculate the momentum transfer to the continuous phase // - do not include gravity impulse dUTrans = -mass*(Unew - U - dt*td.g()); @@ -188,6 +190,9 @@ bool Foam::KinematicParcel::move(TrackData& td) while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL) { + // Apply correction to position for reduced-D cases + meshTools::constrainToMeshCentre(mesh, p.position()); + // Set the Lagrangian time-step scalar dt = min(dtMax, tEnd);