zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Hacking in a sympletic, leapfrog integrator for interparticle forces.

Browse Source 
Conserves kinetic energy correctly.
This commit is contained in:
graham
2009-09-08 19:29:37 +01:00
parent e2c0466b21
commit 2d017d6cc4
2 changed files with 16 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
src/lagrangian/intermediate/clouds/Templates/InteractingKinematicCloud/InteractingKinematicCloud.C
View File

@ -254,10 +254,24 @@ void Foam::InteractingKinematicCloud<ParcelType>::evolve()
resetSourceTerms(); resetSourceTerms();
} }
this->collision().collide(); const scalar deltaT = mesh_.time().deltaT().value();
forAllIter(typename InteractingKinematicCloud<ParcelType>, *this, iter)
{
ParcelType& p = iter();
p.U() += 0.5*deltaT*p.f()/p.mass();
}
Cloud<ParcelType>::move(td); Cloud<ParcelType>::move(td);
this->collision().collide();
forAllIter(typename InteractingKinematicCloud<ParcelType>, *this, iter)
{
ParcelType& p = iter();
p.U() += 0.5*deltaT*p.f()/p.mass();
}
postEvolve(); postEvolve();
} }

3
src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.C
View File

@ -167,14 +167,13 @@ const Foam::vector Foam::InteractingKinematicParcel<ParcelType>::calcVelocity
mass*td.cloud().forces().calcCoupled(cellI, dt, rhoc_, rho, Uc_, U); mass*td.cloud().forces().calcCoupled(cellI, dt, rhoc_, rho, Uc_, U);
const vector FNonCoupled = const vector FNonCoupled =
mass*td.cloud().forces().calcNonCoupled(cellI, dt, rhoc_, rho, Uc_, U); mass*td.cloud().forces().calcNonCoupled(cellI, dt, rhoc_, rho, Uc_, U);
const vector& FCollision = f();
// New particle velocity // New particle velocity
//~~~~~~~~~~~~~~~~~~~~~~ //~~~~~~~~~~~~~~~~~~~~~~
// Update velocity - treat as 3-D // Update velocity - treat as 3-D
const scalar As = this->areaS(d); const scalar As = this->areaS(d);
const vector ap = Uc_ + (FCoupled + FNonCoupled + FCollision + Su)/(utc*As); const vector ap = Uc_ + (FCoupled + FNonCoupled + Su)/(utc*As);
const scalar bp = 6.0*utc/(rho*d); const scalar bp = 6.0*utc/(rho*d);
IntegrationScheme<vector>::integrationResult Ures = IntegrationScheme<vector>::integrationResult Ures =