Dynamics fixed - using notation from Dullweber et al, storing body local angular momentum and torque. Care required when measuring global angular momentum.

This commit is contained in:
graham
2008-11-25 18:14:52 +00:00
parent edfcceb31c
commit 5037380a83
6 changed files with 245 additions and 119 deletions

View File

@ -48,9 +48,29 @@ scalar singleStepTotalPE = 0.0;
scalar singleStepTotalrDotf = 0.0;
vector singleStepCentreOfMass(vector::zero);
{
IDLList<molecule>::iterator mol(molecules.begin());
for
(
mol = molecules.begin();
mol != molecules.end();
++mol
)
{
label molId = mol().id();
scalar molMass(molecules.constProps(molId).mass());
singleStepTotalMass += molMass;
singleStepCentreOfMass += mol().position()*molMass;
}
singleStepCentreOfMass /= singleStepTotalMass;
for
(
mol = molecules.begin();
@ -66,13 +86,19 @@ scalar singleStepTotalrDotf = 0.0;
const vector& molV(mol().v());
const vector& molOmega(mol().omega());
const vector& molOmega(inv(molMoI) & mol().pi());
// Info<< mol().pi()
// << nl << molOmega
// << nl << (molOmega & mol().pi())/(mag(molOmega) * mag(mol().pi()))
// << endl;
vector molPiGlobal = mol().Q() & mol().pi();
singleStepTotalLinearMomentum += molV * molMass;
singleStepTotalAngularMomentum += (molMoI & molOmega);
singleStepTotalMass += molMass;
singleStepTotalAngularMomentum += molPiGlobal
+((mol().position() - singleStepCentreOfMass) ^ (molV * molMass));
if(mag(molV) > singleStepMaxVelocityMag)
{
@ -123,7 +149,8 @@ if (singleStepNMols)
<< "Average linear momentum per mol = "
<< singleStepTotalLinearMomentum/singleStepNMols << nl
<< "Average angular momentum per mol = "
<< singleStepTotalAngularMomentum/singleStepNMols << nl
<< singleStepTotalAngularMomentum << ' '
<< mag(singleStepTotalAngularMomentum)/singleStepNMols << nl
<< "Maximum |velocity| = "
<< singleStepMaxVelocityMag << nl
<< "Average linear KE per mol = "

View File

@ -31,45 +31,36 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::tensor Foam::molecule::rotationTensor(scalar deltaT) const
Foam::tensor Foam::molecule::rotationTensorX(scalar phi) const
{
scalar phi1 = 0.5*deltaT*omega_.x();
tensor U1
return tensor
(
tensor
(
1, 0, 0,
0, Foam::cos(phi1), -Foam::sin(phi1),
0, Foam::sin(phi1), Foam::cos(phi1)
)
1, 0, 0,
0, Foam::cos(phi), -Foam::sin(phi),
0, Foam::sin(phi), Foam::cos(phi)
);
}
scalar phi2 = 0.5*deltaT*omega_.y();
tensor U2
Foam::tensor Foam::molecule::rotationTensorY(scalar phi) const
{
return tensor
(
tensor
(
Foam::cos(phi2), 0, Foam::sin(phi2),
0, 1, 0,
-Foam::sin(phi2), 0, Foam::cos(phi2)
)
Foam::cos(phi), 0, Foam::sin(phi),
0, 1, 0,
-Foam::sin(phi), 0, Foam::cos(phi)
);
}
scalar phi3 = deltaT*omega_.z();
tensor U3
Foam::tensor Foam::molecule::rotationTensorZ(scalar phi) const
{
return tensor
(
tensor
(
Foam::cos(phi3), -Foam::sin(phi3), 0,
Foam::sin(phi3), Foam::cos(phi3), 0,
0, 0, 1
)
Foam::cos(phi), -Foam::sin(phi), 0,
Foam::sin(phi), Foam::cos(phi), 0,
0, 0, 1
);
return (U1 & U2 & U3 & U2 & U1);
}
@ -89,14 +80,17 @@ Foam::molecule::trackData::trackData
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::molecule::move(molecule::trackData& td)
{
td.switchProcessor = false;
td.keepParticle = true;
const constantProperties& constProps(td.molCloud().constProps(id_));
const diagTensor& momentOfInertia(constProps.momentOfInertia());
scalar deltaT = cloud().pMesh().time().deltaT().value();
scalar tEnd = (1.0 - stepFraction())*deltaT;
scalar dtMax = tEnd;
if (td.part() == 0)
{
@ -105,12 +99,15 @@ bool Foam::molecule::move(molecule::trackData& td)
v_ += 0.5*deltaT*a_;
omega_ += 0.5*deltaT*alpha_;
pi_ += 0.5*deltaT*tau_;
}
else if (td.part() == 1)
{
// Leapfrog tracking part
scalar tEnd = (1.0 - stepFraction())*deltaT;
scalar dtMax = tEnd;
while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
{
// set the lagrangian time-step
@ -128,22 +125,40 @@ bool Foam::molecule::move(molecule::trackData& td)
// but after tracking stage, i.e. rotation carried once linear motion
// complete.
R_ = R_ & rotationTensor(deltaT);
tensor R;
setSitePositions(td.molCloud().constProps(id_));
R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
pi_ = R & pi_;
Q_ = Q_ & R.T();
R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
pi_ = R & pi_;
Q_ = Q_ & R.T();
R = rotationTensorZ(deltaT*pi_.z()/momentOfInertia.zz());
pi_ = R & pi_;
Q_ = Q_ & R.T();
R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
pi_ = R & pi_;
Q_ = Q_ & R.T();
R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
pi_ = R & pi_;
Q_ = Q_ & R.T();
setSitePositions(constProps);
}
else if (td.part() == 3)
{
// Second leapfrog velocity adjust part, required after tracking+force
// part
const diagTensor& I(td.molCloud().constProps(id_).momentOfInertia());
scalar m = td.molCloud().constProps(id_).mass();
scalar m = constProps.mass();
a_ = vector::zero;
vector tau(vector::zero);
tau_ = vector::zero;
forAll(siteForces_, s)
{
@ -151,14 +166,12 @@ bool Foam::molecule::move(molecule::trackData& td)
a_ += f/m;
tau += (sitePositions_[s] - position_) ^ f;
tau_ += ((Q_.T() & f) ^ constProps.siteReferencePositions()[s]);
}
alpha_ = R_ & inv(I) & R_.T() & tau;
v_ += 0.5*deltaT*a_;
omega_ += 0.5*deltaT*alpha_;
pi_ += 0.5*deltaT*tau_;
}
else
{
@ -171,6 +184,88 @@ bool Foam::molecule::move(molecule::trackData& td)
return td.keepParticle;
}
// bool Foam::molecule::move(molecule::trackData& td)
// {
// td.switchProcessor = false;
// td.keepParticle = true;
// scalar deltaT = cloud().pMesh().time().deltaT().value();
// scalar tEnd = (1.0 - stepFraction())*deltaT;
// scalar dtMax = tEnd;
// if (td.part() == 0)
// {
// // First leapfrog velocity adjust part, required before tracking+force
// // part
// v_ += 0.5*deltaT*a_;
// omega_ += 0.5*deltaT*alpha_;
// }
// else if (td.part() == 1)
// {
// // Leapfrog tracking part
// while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
// {
// // set the lagrangian time-step
// scalar dt = min(dtMax, tEnd);
// dt *= trackToFace(position() + dt*v_, td);
// tEnd -= dt;
// stepFraction() = 1.0 - tEnd/deltaT;
// }
// }
// else if (td.part() == 2)
// {
// // Leapfrog orientation adjustment, carried out before force calculation
// // but after tracking stage, i.e. rotation carried once linear motion
// // complete.
// R_ = R_ & rotationTensor(deltaT);
// setSitePositions(td.molCloud().constProps(id_));
// }
// else if (td.part() == 3)
// {
// // Second leapfrog velocity adjust part, required after tracking+force
// // part
// const diagTensor& I(td.molCloud().constProps(id_).momentOfInertia());
// scalar m = td.molCloud().constProps(id_).mass();
// a_ = vector::zero;
// vector tau(vector::zero);
// forAll(siteForces_, s)
// {
// const vector& f = siteForces_[s];
// a_ += f/m;
// tau += (sitePositions_[s] - position_) ^ f;
// }
// alpha_ = R_ & inv(I) & R_.T() & tau;
// v_ += 0.5*deltaT*a_;
// omega_ += 0.5*deltaT*alpha_;
// }
// else
// {
// FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
// << td.part()
// << " is an invalid part of the integration method."
// << abort(FatalError);
// }
// return td.keepParticle;
// }
void Foam::molecule::transformProperties(const tensor& T)
{}
@ -187,7 +282,7 @@ void Foam::molecule::transformProperties(const vector& separation)
void Foam::molecule::setSitePositions(const constantProperties& constProps)
{
sitePositions_ = position_ + (R_ & constProps.siteReferencePositions());
sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions());
}

View File

@ -177,15 +177,15 @@ private:
// -# Don't go scalar-label, scalar-label, because in 64bit mode,
// the labels will be padded by 4bytes.
tensor R_;
tensor Q_;
vector v_;
vector a_;
vector omega_;
vector pi_;
vector alpha_;
vector tau_;
List<vector> siteForces_;
@ -204,7 +204,11 @@ private:
// Private Member Functions
tensor rotationTensor(scalar deltaT) const;
tensor rotationTensorX(scalar deltaT) const;
tensor rotationTensorY(scalar deltaT) const;
tensor rotationTensorZ(scalar deltaT) const;
public:
@ -218,11 +222,11 @@ public:
const Cloud<molecule>& c,
const vector& position,
const label celli,
const tensor& R,
const tensor& Q,
const vector& v,
const vector& a,
const vector& omega,
const vector& alpha,
const vector& pi,
const vector& tau,
const vector& specialPosition,
const constantProperties& constProps,
const label special,
@ -257,8 +261,8 @@ public:
// Access
inline const tensor& R() const;
inline tensor& R();
inline const tensor& Q() const;
inline tensor& Q();
inline const vector& v() const;
inline vector& v();
@ -266,11 +270,11 @@ public:
inline const vector& a() const;
inline vector& a();
inline const vector& omega() const;
inline vector& omega();
inline const vector& pi() const;
inline vector& pi();
inline const vector& alpha() const;
inline vector& alpha();
inline const vector& tau() const;
inline vector& tau();
inline const List<vector>& siteForces() const;
inline List<vector>& siteForces();

View File

@ -66,11 +66,11 @@ inline Foam::molecule::molecule
const Cloud<molecule>& c,
const vector& position,
const label celli,
const tensor& R,
const tensor& Q,
const vector& v,
const vector& a,
const vector& omega,
const vector& alpha,
const vector& pi,
const vector& tau,
const vector& specialPosition,
const constantProperties& constProps,
const label special,
@ -79,11 +79,11 @@ inline Foam::molecule::molecule
)
:
Particle<molecule>(c, position, celli),
R_(R),
Q_(Q),
v_(v),
a_(a),
omega_(omega),
alpha_(alpha),
pi_(pi),
tau_(tau),
siteForces_(constProps.nSites(), vector::zero),
sitePositions_(constProps.nSites()),
specialPosition_(specialPosition),
@ -250,15 +250,15 @@ inline Foam::label Foam::molecule::trackData::part() const
// * * * * * * * * * * * * molecule Member Functions * * * * * * * * * * * * //
inline const Foam::tensor& Foam::molecule::R() const
inline const Foam::tensor& Foam::molecule::Q() const
{
return R_;
return Q_;
}
inline Foam::tensor& Foam::molecule::R()
inline Foam::tensor& Foam::molecule::Q()
{
return R_;
return Q_;
}
@ -286,27 +286,27 @@ inline Foam::vector& Foam::molecule::a()
}
inline const Foam::vector& Foam::molecule::omega() const
inline const Foam::vector& Foam::molecule::pi() const
{
return omega_;
return pi_;
}
inline Foam::vector& Foam::molecule::omega()
inline Foam::vector& Foam::molecule::pi()
{
return omega_;
return pi_;
}
inline const Foam::vector& Foam::molecule::alpha() const
inline const Foam::vector& Foam::molecule::tau() const
{
return alpha_;
return tau_;
}
inline Foam::vector& Foam::molecule::alpha()
inline Foam::vector& Foam::molecule::tau()
{
return alpha_;
return tau_;
}

View File

@ -38,11 +38,11 @@ Foam::molecule::molecule
)
:
Particle<molecule>(cloud, is),
R_(tensor::zero),
Q_(tensor::zero),
v_(vector::zero),
a_(vector::zero),
omega_(vector::zero),
alpha_(vector::zero),
pi_(vector::zero),
tau_(vector::zero),
siteForces_(List<vector>(0,vector::zero)),
sitePositions_(List<vector>(0,vector::zero)),
specialPosition_(vector::zero),
@ -58,11 +58,11 @@ Foam::molecule::molecule
{
if (is.format() == IOstream::ASCII)
{
is >> R_;
is >> Q_;
is >> v_;
is >> a_;
is >> omega_;
is >> alpha_;
is >> pi_;
is >> tau_;
is >> siteForces_;
is >> sitePositions_;
is >> specialPosition_;
@ -75,12 +75,12 @@ Foam::molecule::molecule
{
is.read
(
reinterpret_cast<char*>(&R_),
sizeof(R_)
reinterpret_cast<char*>(&Q_),
sizeof(Q_)
+ sizeof(v_)
+ sizeof(a_)
+ sizeof(omega_)
+ sizeof(alpha_)
+ sizeof(pi_)
+ sizeof(tau_)
+ sizeof(siteForces_)
+ sizeof(sitePositions_)
+ sizeof(specialPosition_)
@ -107,8 +107,8 @@ void Foam::molecule::readFields(moleculeCloud& mC)
return;
}
IOField<tensor> R(mC.fieldIOobject("R", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, R);
IOField<tensor> Q(mC.fieldIOobject("Q", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, Q);
IOField<vector> v(mC.fieldIOobject("v", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, v);
@ -116,11 +116,11 @@ void Foam::molecule::readFields(moleculeCloud& mC)
IOField<vector> a(mC.fieldIOobject("a", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, a);
IOField<vector> omega(mC.fieldIOobject("omega", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, omega);
IOField<vector> pi(mC.fieldIOobject("pi", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, pi);
IOField<vector> alpha(mC.fieldIOobject("alpha", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, alpha);
IOField<vector> tau(mC.fieldIOobject("tau", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, tau);
IOField<vector> specialPosition
(
@ -139,11 +139,11 @@ void Foam::molecule::readFields(moleculeCloud& mC)
{
molecule& mol = iter();
mol.R_ = R[i];
mol.Q_ = Q[i];
mol.v_ = v[i];
mol.a_ = a[i];
mol.omega_ = omega[i];
mol.alpha_ = alpha[i];
mol.pi_ = pi[i];
mol.tau_ = tau[i];
mol.specialPosition_ = specialPosition[i];
mol.special_ = special[i];
mol.id_ = id[i];
@ -158,11 +158,11 @@ void Foam::molecule::writeFields(const moleculeCloud& mC)
label np = mC.size();
IOField<tensor> R(mC.fieldIOobject("R", IOobject::NO_READ), np);
IOField<tensor> Q(mC.fieldIOobject("Q", IOobject::NO_READ), np);
IOField<vector> v(mC.fieldIOobject("v", IOobject::NO_READ), np);
IOField<vector> a(mC.fieldIOobject("a", IOobject::NO_READ), np);
IOField<vector> omega(mC.fieldIOobject("omega", IOobject::NO_READ), np);
IOField<vector> alpha(mC.fieldIOobject("alpha", IOobject::NO_READ), np);
IOField<vector> pi(mC.fieldIOobject("pi", IOobject::NO_READ), np);
IOField<vector> tau(mC.fieldIOobject("tau", IOobject::NO_READ), np);
IOField<vector> specialPosition
(
mC.fieldIOobject("specialPosition", IOobject::NO_READ),
@ -176,22 +176,22 @@ void Foam::molecule::writeFields(const moleculeCloud& mC)
{
const molecule& mol = iter();
R[i] = mol.R_;
Q[i] = mol.Q_;
v[i] = mol.v_;
a[i] = mol.a_;
omega[i] = mol.omega_;
alpha[i] = mol.alpha_;
pi[i] = mol.pi_;
tau[i] = mol.tau_;
specialPosition[i] = mol.specialPosition_;
special[i] = mol.special_;
id[i] = mol.id_;
i++;
}
R.write();
Q.write();
v.write();
a.write();
omega.write();
alpha.write();
pi.write();
tau.write();
specialPosition.write();
special.write();
id.write();
@ -207,11 +207,11 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const molecule& mol)
os << token::SPACE << static_cast<const Particle<molecule>&>(mol)
<< token::SPACE << mol.face()
<< token::SPACE << mol.stepFraction()
<< token::SPACE << mol.R_
<< token::SPACE << mol.Q_
<< token::SPACE << mol.v_
<< token::SPACE << mol.a_
<< token::SPACE << mol.omega_
<< token::SPACE << mol.alpha_
<< token::SPACE << mol.pi_
<< token::SPACE << mol.tau_
<< token::SPACE << mol.siteForces_
<< token::SPACE << mol.sitePositions_
<< token::SPACE << mol.specialPosition_
@ -225,12 +225,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const molecule& mol)
os << static_cast<const Particle<molecule>&>(mol);
os.write
(
reinterpret_cast<const char*>(&mol.R_),
sizeof(mol.R_)
reinterpret_cast<const char*>(&mol.Q_),
sizeof(mol.Q_)
+ sizeof(mol.v_)
+ sizeof(mol.a_)
+ sizeof(mol.omega_)
+ sizeof(mol.alpha_)
+ sizeof(mol.pi_)
+ sizeof(mol.tau_)
+ sizeof(mol.siteForces_)
+ sizeof(mol.sitePositions_)
+ sizeof(mol.specialPosition_)

View File

@ -485,7 +485,7 @@ Foam::moleculeCloud::moleculeCloud
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const Cloud<molecule>& cloud = *this;
vector position1 = vector(0,0,1.5e-9);
vector position1 = vector(0,0,4e-9);
addParticle
(
@ -495,9 +495,9 @@ Foam::moleculeCloud::moleculeCloud
position1,
mesh_.findCell(position1),
tensor(1, 0, 0, 0, 1, 0, 0, 0, 1),
vector(34, 12, -65),
vector(0,0,-150),
vector::zero,
vector(0.2, -0.012, 0.32),
vector(1e-35, 2e-35, 3e-35),
vector::zero,
vector::zero,
constProps(0),
@ -506,7 +506,7 @@ Foam::moleculeCloud::moleculeCloud
)
);
vector position2 = vector(0,0,0);
vector position2 = vector(0,0,-4e-9);
addParticle
(
@ -516,9 +516,9 @@ Foam::moleculeCloud::moleculeCloud
position2,
mesh_.findCell(position2),
tensor(1, 0, 0, 0, 1, 0, 0, 0, 1),
vector(0,0,150),
vector::zero,
vector::zero,
vector(-0.4, 0.015, -0.2),
vector(1e-35, 2e-35, 3e-35),
vector::zero,
vector::zero,
constProps(1),