From fa8f316eccade3896e23cfe000d0507d0107e47d Mon Sep 17 00:00:00 2001 From: graham Date: Mon, 12 Jan 2009 17:35:57 +0000 Subject: [PATCH] Adding pressure measurement via r dot f and temperature measurement and control including rotational energy. Adding random orientation on initialisation. Modifying constantProperties construction to detect point masses. Tidy up hitWallPatch function to remove commented out stochastic wall code. --- .../solvers/molecularDynamics/mdFoam/mdFoam.C | 20 +++++ .../molecule/mdTools/temperatureAndPressure.H | 20 +++-- .../mdTools/temperatureAndPressureVariables.H | 6 +- .../molecule/molecule/molecule.C | 36 ++------- .../molecule/molecule/moleculeI.H | 81 ++++++++++--------- .../molecule/moleculeCloud/moleculeCloud.C | 23 +++++- .../molecule/moleculeCloud/moleculeCloudI.H | 20 ++--- 7 files changed, 114 insertions(+), 92 deletions(-) diff --git a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C index 4357c5c0b1..f65424b437 100644 --- a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C +++ b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C @@ -46,20 +46,40 @@ int main(int argc, char *argv[]) # include "temperatureAndPressureVariables.H" + label nAveragingSteps = 0; + Info << "\nStarting time loop\n" << endl; while (runTime.run()) { runTime++; + nAveragingSteps++; + Info << "Time = " << runTime.timeName() << endl; molecules.evolve(); # include "meanMomentumEnergyAndNMols.H" +# include "temperatureAndPressure.H" + + if (runTime.outputTime()) + { + molecules.applyConstraintsAndThermostats + ( + 485.4, + averageTemperature + ); + } + runTime.write(); + if (runTime.outputTime()) + { + nAveragingSteps = 0; + } + Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H index 24d6fb88e5..ddf99c30d5 100644 --- a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H +++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H @@ -33,11 +33,13 @@ Description \*---------------------------------------------------------------------------*/ -accumulatedTotalMomentum += singleStepTotalMomentum; +accumulatedTotalLinearMomentum += singleStepTotalLinearMomentum; accumulatedTotalMass += singleStepTotalMass; -accumulatedTotalKE += singleStepTotalKE; +accumulatedTotalLinearKE += singleStepTotalLinearKE; + +accumulatedTotalAngularKE += singleStepTotalAngularKE; accumulatedTotalPE += singleStepTotalPE; @@ -55,12 +57,12 @@ if (runTime.outputTime()) averageTemperature = ( - 2.0/(3.0 * moleculeCloud::kb * accumulatedNMols) + 2.0/(6.0 * moleculeCloud::kb * accumulatedNMols) * ( - accumulatedTotalKE + accumulatedTotalLinearKE + accumulatedTotalAngularKE - - 0.5*magSqr(accumulatedTotalMomentum)/accumulatedTotalMass + 0.5*magSqr(accumulatedTotalLinearMomentum)/accumulatedTotalMass ) ); @@ -82,7 +84,7 @@ if (runTime.outputTime()) Info << "----------------------------------------" << nl << "Averaged properties" << nl << "Average |velocity| = " - << mag(accumulatedTotalMomentum)/accumulatedTotalMass + << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass << " m/s" << nl << "Average temperature = " << averageTemperature << " K" << nl @@ -98,11 +100,13 @@ if (runTime.outputTime()) // reset counters - accumulatedTotalMomentum = vector::zero; + accumulatedTotalLinearMomentum = vector::zero; accumulatedTotalMass = 0.0; - accumulatedTotalKE = 0.0; + accumulatedTotalLinearKE = 0.0; + + accumulatedTotalAngularKE = 0.0; accumulatedTotalPE = 0.0; diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H index 283535c7a4..45b5c0d221 100644 --- a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H +++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H @@ -30,11 +30,13 @@ Description \*---------------------------------------------------------------------------*/ -vector accumulatedTotalMomentum(vector::zero); +vector accumulatedTotalLinearMomentum(vector::zero); scalar accumulatedTotalMass = 0.0; -scalar accumulatedTotalKE = 0.0; +scalar accumulatedTotalAngularKE = 0.0; + +scalar accumulatedTotalLinearKE = 0.0; scalar accumulatedTotalPE = 0.0; diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C index 2d8cd8f00b..5f90b7856e 100644 --- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C +++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C @@ -246,38 +246,12 @@ void Foam::molecule::hitWallPatch nw /= mag(nw); scalar vn = v_ & nw; -// vector vt = v_ - vn*nw; - -// Random rand(clock::getTime()); - -// scalar tmac = 0.8; - -// scalar wallTemp = 2.5; - -// if (rand.scalar01() < tmac) -// { -// // Diffuse reflection -// -// vector tw1 = vt/mag(vt); -// -// vector tw2 = nw ^ tw1; -// -// V_ = sqrt(wallTemp/mass_)*rand.GaussNormal()*tw1 -// + sqrt(wallTemp/mass_)*rand.GaussNormal()*tw2 -// - mag(sqrt(wallTemp/mass_)*rand.GaussNormal())*nw; -// } - -// else -// { - // Specular reflection - - if (vn > 0) - { - v_ -= 2*vn*nw; - } - -// } + // Specular reflection + if (vn > 0) + { + v_ -= 2*vn*nw; + } } diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H index 7d3f5c3aa7..c1d69a87c4 100644 --- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H +++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H @@ -63,46 +63,47 @@ inline Foam::molecule::constantProperties::constantProperties mass_ = sum(siteMasses); + // Calculate the centre of mass of the body and subtract it from each + // position + + vector centreOfMass(vector::zero); + + forAll(siteReferencePositions_, i) { - // Calculate the centre of mass of the body and subtract it from each - // position + centreOfMass += siteReferencePositions_[i]*siteMasses[i]; + } - vector centreOfMass(vector::zero); + centreOfMass /= mass_; - forAll(siteReferencePositions_, i) - { - centreOfMass += siteReferencePositions_[i]*siteMasses[i]; - } + siteReferencePositions_ -= centreOfMass; - centreOfMass /= mass_; + // Calculate the inertia tensor in the current orientation - siteReferencePositions_ -= centreOfMass; + tensor momOfInertia(tensor::zero); - // Calculate the inertia tensor in the current orientation + forAll(siteReferencePositions_, i) + { + const vector& p(siteReferencePositions_[i]); - tensor I(tensor::zero); - - forAll(siteReferencePositions_, i) - { - const vector& p(siteReferencePositions_[i]); - - I += siteMasses[i]*tensor - ( - p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(), - -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(), - -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y() - ); - } + momOfInertia += siteMasses[i]*tensor + ( + p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(), + -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(), + -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y() + ); + } + if(eigenValues(momOfInertia).z() > VSMALL) + { // Normalise the inertia tensor magnitude to avoid SMALL numbers in the // components causing problems - I /= eigenValues(I).z(); + momOfInertia /= eigenValues(momOfInertia).z(); - tensor e = eigenVectors(I); + tensor e = eigenVectors(momOfInertia); - // Calculate the transformation between the principle axes to the global - // axes + // Calculate the transformation between the principle axes to the + // global axes tensor Q = vector(1,0,0)*e.x() + vector(0,1,0)*e.y() + vector(0,0,1)*e.z(); @@ -110,13 +111,11 @@ inline Foam::molecule::constantProperties::constantProperties // Transform the site positions siteReferencePositions_ = (Q & siteReferencePositions_); - } - { // Recalculating the moment of inertia with the new site positions - // The rotation was around the centre of mass but remove any components - // that have crept in due to floating point errors + // The rotation was around the centre of mass but remove any + // components that have crept in due to floating point errors vector centreOfMass(vector::zero); @@ -129,16 +128,16 @@ inline Foam::molecule::constantProperties::constantProperties siteReferencePositions_ -= centreOfMass; - // Calculate the moment of inertia in the principle component reference - // frame + // Calculate the moment of inertia in the principle component + // reference frame - tensor I(tensor::zero); + tensor momOfInertia(tensor::zero); forAll(siteReferencePositions_, i) { const vector& p(siteReferencePositions_[i]); - I += siteMasses[i]*tensor + momOfInertia += siteMasses[i]*tensor ( p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(), -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(), @@ -146,9 +145,19 @@ inline Foam::molecule::constantProperties::constantProperties ); } - momentOfInertia_ = diagTensor(I.xx(), I.yy(), I.zz()); + momentOfInertia_ = diagTensor + ( + momOfInertia.xx(), + momOfInertia.yy(), + momOfInertia.zz() + ); } + else + { + Info<< nl << "Molecule " << dict.name() << " is a point mass." << endl; + momentOfInertia_ = diagTensor(0, 0, 0); + } } diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C index c15a629417..1d0e2da8d1 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C +++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C @@ -982,6 +982,25 @@ void Foam::moleculeCloud::createMolecule cP.momentOfInertia() ); + scalar phi(rndGen_.scalar01()*mathematicalConstant::twoPi); + + scalar theta(rndGen_.scalar01()*mathematicalConstant::twoPi); + + scalar psi(rndGen_.scalar01()*mathematicalConstant::twoPi); + + const tensor Q + ( + cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi), + cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi), + sin(psi)*sin(theta), + - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi), + - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi), + cos(psi)*sin(theta), + sin(theta)*sin(phi), + - sin(theta)*cos(phi), + cos(theta) + ); + addParticle ( new molecule @@ -989,7 +1008,7 @@ void Foam::moleculeCloud::createMolecule cloud, position, cell, - I, + Q, v, vector::zero, pi, @@ -1124,6 +1143,8 @@ void Foam::moleculeCloud::applyConstraintsAndThermostats for (mol = this->begin(); mol != this->end(); ++mol) { mol().v() *= temperatureCorrectionFactor; + + mol().pi() *= temperatureCorrectionFactor; } } diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H index ca095562cd..72e5617a8f 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H +++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H @@ -91,11 +91,9 @@ inline void Foam::moleculeCloud::evaluatePair molJ->potentialEnergy() += 0.5*potentialEnergy; - // What to do here? + molI->rf() += rsIsJ * fsIsJ; - // molI->rf() += rIJ * fIJ; - - // molJ->rf() += rIJ * fIJ; + molJ->rf() += rsIsJ * fsIsJ; } } @@ -128,11 +126,9 @@ inline void Foam::moleculeCloud::evaluatePair molJ->potentialEnergy() += 0.5*potentialEnergy; - // What to do here? + molI->rf() += rsIsJ * fsIsJ; - // molI->rf() += rIJ * fIJ; - - // molJ->rf() += rIJ * fIJ; + molJ->rf() += rsIsJ * fsIsJ; } } } @@ -200,9 +196,7 @@ inline void Foam::moleculeCloud::evaluatePair molReal->potentialEnergy() += 0.5*potentialEnergy; - // What to do here? - - // molReal->rf() += rRealRef * fRealRef; + molReal->rf() += rsRealsRef * fsRealsRef; } } @@ -233,9 +227,7 @@ inline void Foam::moleculeCloud::evaluatePair molReal->potentialEnergy() += 0.5*potentialEnergy; - // What to do here? - - // molReal->rf() += rRealRef * fRealRef; + molReal->rf() += rsRealsRef * fsRealsRef; } } }