Commit 24984608 authored by Igor Ostanin's avatar Igor Ostanin

A prototype for potential-based bond model

parent aade9649
Pipeline #5029 failed with stage
in 112 minutes and 25 seconds
......@@ -55,10 +55,10 @@ typedef boost::tuple<Sphere, Plane, LRCapsule> BodyTypeTuple ;
// For the units system used see PFC 5 implementation docs
// Timestep is fixed to 20 fs.
real_t dt = 20.;
real_t dt = 30.;
// vdW interaction radius w/r to inertial segment radius
real_t CutoffFactor = 100.;
real_t CutoffFactor = 4.;
// CNT lattice numbers (m,n) CNT
real_t mm = 10.;
......@@ -97,7 +97,7 @@ real_t mass_T = ro * T * M_C * 104.397;
real_t vol = (4./3.) * math::PI * (R_CNT * R_CNT * R_CNT) + math::PI * R_CNT * R_CNT * T;
// Density of a capsule
real_t dens = mass_T / vol;
real_t dens = 0.1 * mass_T / vol;
// Parallel bond parameters;
real_t knorm = (En * 0.006242) / T;
......@@ -185,29 +185,65 @@ int main( int argc, char ** argv )
//! [Gas]
//! DEPOSITION OF CNTS
uint_t Num_CNTs = 4;
uint_t Num_Segs = 20;
real_t stretch = 1.4;
T = T + margin; // bugfix for cap segfaults
uint_t numParticles = uint_c(0);
Vec3 pos(0.0,0.0,0.0);
// Create few parallel CNTs
for (int tube = 0; tube < 4; tube++)
for (int tube = 0; tube < Num_CNTs; tube++)
{
for (int segment = 0; segment < 4; segment++)
for (int segment = 0; segment < Num_Segs; segment++)
{
pos[0] = 20. * segment;
// Flat deposition - test stretching
//pos[0] = 2. * R_CNT * stretch * segment;
//pos[1] = 20 * tube;
// Circular deposition - test bending
real_t R_B = 200;
pos[0] = R_B * std::cos( ( 2 * R_CNT / R_B ) * segment);
pos[2] = R_B * std::sin( ( 2 * R_CNT / R_B ) * segment);
pos[1] = 20 * tube;
WALBERLA_LOG_DEVEL(pos[0]);
LRCapsuleID sp = createLRCapsule( *globalBodyStorage, *forest, storageID, 0, pos, outer_radius, inner_radius, length, material, false, true, false);
if (sp != NULL) sp->setSegmentID(segment);
if (sp != NULL) sp->setClusterID(tube);
if (sp != NULL) sp->rotate(0., 0., 1., math::PI/4.*segment );
if (sp != NULL) sp->rotate(0., 1., 0.,- 0.5 * math::PI - ( 2 * R_CNT / R_B ) * segment );
if (sp != NULL) numParticles++;
if (sp != NULL) WALBERLA_LOG_DEVEL(sp);
// if (sp != NULL) WALBERLA_LOG_DEVEL(sp);
}
}
WALBERLA_LOG_DEVEL(numParticles);
// for (auto it = grid_generator::SCIterator(math::AABB(0,0,0,20,20,20), Vector3<real_t>(spacing, spacing, spacing) * real_c(0.5), spacing); it != grid_generator::SCIterator(); ++it)
// {
// LRCapsuleID sp;
......@@ -246,7 +282,7 @@ int main( int argc, char ** argv )
{
WALBERLA_LOG_DEVEL_ON_ROOT( "Timestep " << i << " / " << simulationSteps );
}
//WALBERLA_LOG_DEVEL_ON_ROOT( "!!!");
cr.timestep( real_c(dt) );
syncNextNeighbors<BodyTypeTuple>(*forest, storageID);
vtkWriter->write(true);
......
......@@ -64,7 +64,7 @@ namespace pe {
// For the units system used see PFC 5 implementation docs
// Timestep is fixed to 20 fs.
//real_t dt = 20.;
real_t dt = 100.;
// vdW interaction radius w/r to inertial segment radius
real_t CutoffFactor = 4.;
......@@ -273,14 +273,12 @@ void DEM::resolveContact( ContactID c, real_t dt) const
real_t delta( -c->getDistance() );
// Contact normal
const Vec3 n( (seg1->getPosition() - seg2->getPosition()).getNormalized() );
// CM -specific calculations
if (1) // Vdw contact
if (0) // Vdw contact
{
// WALBERLA_LOG_DEVEL( "Contact model type: vdW " );
......@@ -525,9 +523,14 @@ void DEM::resolveContact( ContactID c, real_t dt) const
else if ( L<=2*R_*1.2*TH ) { // Small distance
//WALBERLA_LOG_DEVEL( "Small distance");
seg1->addForce( n);
seg2->addForce(-n);
}
else { // big distance
seg1->addForce( n*0);
seg2->addForce(-n*0);
}
//==========================================================
......@@ -539,7 +542,7 @@ void DEM::resolveContact( ContactID c, real_t dt) const
//WALBERLA_LOG_DEVEL( "Contact model type: zero " );
}
if (0) // Parallel bond contact
if (0) // Parallel bond - invalid model as we know
{
......@@ -579,7 +582,6 @@ void DEM::resolveContact( ContactID c, real_t dt) const
seg2->addTorque(-(dMn + dMs));
}
if (0)
{
......@@ -601,6 +603,56 @@ void DEM::resolveContact( ContactID c, real_t dt) const
}
if (1) // Potential bond model
{
if ((CID_1==CID_2)&&(abs(SID_1 - SID_2) == 1))
{
//WALBERLA_LOG_DEVEL( "POTENTIAL BASED MODEL...");
real_t k_n = knorm;
double L0 = 2*R_CNT;
double L = ( seg1->getPosition() - seg2->getPosition() ).length();
Vec3 force = - k_n * (L - L0) * n;
seg1->addForce( force);
seg2->addForce(-force);
// Moment interaction - bending
Vec3 b1 = seg1->getRotation() * Vec3(1.0, 0.0, 0.0);
Vec3 b2 = seg2->getRotation() * Vec3(1.0, 0.0, 0.0);
real_t cm = 1;
real_t cb = 2;
Vec3 mom = cm * (b1 % b2);
seg1->addTorque( mom);
seg2->addTorque(-mom);
// Shear forces
Vec3 r12 = seg1->getPosition() - seg2->getPosition();
Vec3 sforce = cb * (b1 * (b1 * r12) - r12);
seg1->addForce( sforce);
seg2->addForce(-sforce);
// Viscous damping
real_t damp = 100.0;
Vec3 relVelocity = seg1->getLinearVel() - seg2->getLinearVel();
Vec3 visc_force = - damp * relVelocity;
//WALBERLA_LOG_DEVEL( "Visc force: = " << visc_force );
seg1->addForce( visc_force);
seg2->addForce(-visc_force);
}
}
}
//*************************************************************************************************
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment