New question #290129 on Yade: https://answers.launchpad.net/yade/+question/290129
Hi all, I try to calculate the external work in triaxial compression test. In fact, i found that in the script TriaxialStressController.cpp the external work is calculated when the loading is done by force with the function "controlExternalStress" . But when the loading is done with the strain rate, the external work is not calculated. How can I calculate the external work when i load with the strain rate. Best regards. Jabrane. This is the script TriaxialStressController.cpp #include"TriaxialStressController.hpp" #include<pkg/common/Sphere.hpp> #include<pkg/common/Box.hpp> #include<pkg/dem/ScGeom.hpp> #include<pkg/dem/FrictPhys.hpp> #include<core/State.hpp> #include<assert.h> #include<core/Scene.hpp> #include<pkg/dem/Shop.hpp> #include<core/Clump.hpp> #ifdef FLOW_ENGINE //#include<pkg/pfv/FlowEngine.hpp> #include "FlowEngine_FlowEngineT.hpp" #endif CREATE_LOGGER(TriaxialStressController); YADE_PLUGIN((TriaxialStressController)); TriaxialStressController::~TriaxialStressController(){} Vector3r TriaxialStressController::getStress(int boundId) {assert (boundId>=0 && boundId<=5); return stress[boundId];} Vector3r TriaxialStressController::getStrainRate() { return Vector3r ( (Body::byId(wall_right_id,scene)->state->vel[0]-Body::byId(wall_left_id,scene)->state->vel[0])/width, (Body::byId(wall_top_id,scene)->state->vel[1]-Body::byId(wall_bottom_id,scene)->state->vel[1])/height, (Body::byId(wall_front_id,scene)->state->vel[2]-Body::byId(wall_back_id,scene)->state->vel[2])/depth ); } void TriaxialStressController::updateStiffness() { Real fluidStiffness = 0.; #ifdef FLOW_ENGINE FOREACH(const shared_ptr<Engine> e, Omega::instance().getScene()->engines) { if (e->getClassName() == "FlowEngine") { TemplateFlowEngine_FlowEngineT<FlowCellInfo_FlowEngineT,FlowVertexInfo_FlowEngineT>* flow = dynamic_cast<TemplateFlowEngine_FlowEngineT<FlowCellInfo_FlowEngineT,FlowVertexInfo_FlowEngineT>*>(e.get()); if ( (flow->fluidBulkModulus > 0) && (!(flow->dead)) ) fluidStiffness = flow->fluidBulkModulus/porosity; } } #endif for (int i=0; i<6; ++i) stiffness[i] = 0; InteractionContainer::iterator ii = scene->interactions->begin(); InteractionContainer::iterator iiEnd = scene->interactions->end(); for( ; ii!=iiEnd ; ++ii ) if ((*ii)->isReal()) { const shared_ptr<Interaction>& contact = *ii; Real fn = (static_cast<FrictPhys*> (contact->phys.get()))->normalForce.norm(); if (fn!=0) { int id1 = contact->getId1(), id2 = contact->getId2(); for (int index=0; index<6; ++index) if ( wall_id[index]==id1 || wall_id[index]==id2 ) { FrictPhys* currentContactPhysics = static_cast<FrictPhys*> ( contact->phys.get() ); stiffness[index] += currentContactPhysics->kn; } } } if (fluidStiffness > 0) { stiffness[0] += fluidStiffness*width*depth/height; stiffness[1] += fluidStiffness*width*depth/height; stiffness[2] += fluidStiffness*height*depth/width; stiffness[3] += fluidStiffness*height*depth/width; stiffness[4] += fluidStiffness*width*height/depth; stiffnessWallFront = stiffness[5] += fluidStiffness*width*height/depth; } } void TriaxialStressController::controlExternalStress(int wall, Vector3r resultantForce, State* p, Real wall_max_vel) // controls walls such that Sum Forces from Sample on Wall = resultantForce { scene->forces.sync(); Real translation=normal[wall].dot(getForce(scene,wall_id[wall])-resultantForce); const bool log=false; if(log) LOG_DEBUG("wall="<<wall<<" actualForce="<<getForce(scene,wall_id[wall])<<", resultantForce="<<resultantForce<<", translation="<<translation); if (translation!=0) { if (stiffness[wall]!=0) { translation /= stiffness[wall]; stiffnessWallLeft = stiffness[0]; stiffnessWallRight = stiffness[1]; stiffnessWallBottom = stiffness[2]; stiffnessWallTop = stiffness[3]; stiffnessWallBack = stiffness[4]; stiffnessWallFront = stiffness[5]; if(log) TRVAR2(translation,wall_max_vel*scene->dt) translation = std::min( std::abs(translation), wall_max_vel*scene->dt ) * Mathr::Sign(translation); getForceLeft=normal[wall_left].dot(getForce(scene,wall_id[wall_left]));//JH getForceRight=normal[wall_right].dot(getForce(scene,wall_id[wall_right]));//JH getForceTop=normal[wall_top].dot(getForce(scene,wall_id[wall_top]));//JH getForceBottom=normal[wall_bottom].dot(getForce(scene,wall_id[wall_bottom]));//JH getForceFront=normal[wall_front].dot(getForce(scene,wall_id[wall_front]));//JH getForceBack=normal[wall_back].dot(getForce(scene,wall_id[wall_back])); //JH ResultantForce= normal[wall].dot(resultantForce);//JH deltaForceLeft=normal[wall_left].dot(getForce(scene,wall_id[wall_left])-resultantForce);//JH deltaForceRight=normal[wall_right].dot(getForce(scene,wall_id[wall_right])-resultantForce);//JH deltaForceBottom=normal[wall_bottom].dot(getForce(scene,wall_id[wall_bottom])-resultantForce);//JH deltaForceTop=normal[wall_top].dot(getForce(scene,wall_id[wall_top])-resultantForce);//JH deltaForceBack=normal[wall_back].dot(getForce(scene,wall_id[wall_back])-resultantForce);//JH deltaForceFront=normal[wall_front].dot(getForce(scene,wall_id[wall_front])-resultantForce);//JH deltaForce0=getForceLeft-getForceRight;//JH deltaForce1=getForceTop-getForceBottom;//JH deltaForce2=getForceBack-getForceFront;//JH } else translation = wall_max_vel * Mathr::Sign(translation)*scene->dt; } previousTranslation[wall] = (1-stressDamping)*translation*normal[wall] + 0.8*previousTranslation[wall];// formula for "steady-flow" evolution with fluctuations dampingWorkWall+=stressDamping*(translation*normal[wall]-previousTranslation[wall]).dot(getForce(scene,wall_id[wall]));// JH translationLeft=(previousTranslation[wall]).dot(normal[wall_left]);//JH translationRight=(previousTranslation[wall]).dot(normal[wall_right]);//JH translationTop=(previousTranslation[wall]).dot(normal[wall_top]);//JH translationBottom=(previousTranslation[wall]).dot(normal[wall_bottom]);//JH translationFront=(previousTranslation[wall]).dot(normal[wall_front]);//JH translationBack=(previousTranslation[wall]).dot(normal[wall_back]);//JH //Don't update position since Newton is doing that starting from bzr2612 // p->se3.position += previousTranslation[wall]; externalWork += previousTranslation[wall].dot(getForce(scene,wall_id[wall])); // this is important is using VelocityBins. Otherwise the motion is never detected. Related to https://bugs.launchpad.net/yade/+bug/398089 p->vel=previousTranslation[wall]/scene->dt; //if(log)TRVAR2(previousTranslation,p->se3.position); } void TriaxialStressController::action() { // sync thread storage of ForceContainer scene->forces.sync(); if (first) {// sync boundaries ids in the table wall_id[wall_bottom] = wall_bottom_id; wall_id[wall_top] = wall_top_id; wall_id[wall_left] = wall_left_id; wall_id[wall_right] = wall_right_id; wall_id[wall_front] = wall_front_id; wall_id[wall_back] = wall_back_id;} if(thickness<0) thickness=2.0*YADE_PTR_CAST<Box>(Body::byId(wall_bottom_id,scene)->shape)->extents.y(); State* p_bottom=Body::byId(wall_bottom_id,scene)->state.get(); State* p_top=Body::byId(wall_top_id,scene)->state.get(); State* p_left=Body::byId(wall_left_id,scene)->state.get(); State* p_right=Body::byId(wall_right_id,scene)->state.get(); State* p_front=Body::byId(wall_front_id,scene)->state.get(); State* p_back=Body::byId(wall_back_id,scene)->state.get(); height = p_top->se3.position.y() - p_bottom->se3.position.y() - thickness; width = p_right->se3.position.x() - p_left->se3.position.x() - thickness; depth = p_front->se3.position.z() - p_back->se3.position.z() - thickness; boxVolume = height * width * depth; if ( (first) || (updatePorosity) ) { BodyContainer::iterator bi = scene->bodies->begin(); BodyContainer::iterator biEnd = scene->bodies->end(); particlesVolume = 0; for ( ; bi!=biEnd; ++bi ) { const shared_ptr<Body>& b = *bi; if (b->isClump()) { const shared_ptr<Clump>& clump = YADE_PTR_CAST<Clump>(b->shape); const shared_ptr<Body>& member = Body::byId(clump->members.begin()->first,scene); particlesVolume += b->state->mass / member->material->density; } else if (b->isDynamic() && !b->isClumpMember()) { const shared_ptr<Sphere>& sphere = YADE_PTR_CAST<Sphere> ( b->shape ); particlesVolume += 1.3333333*Mathr::PI*pow ( sphere->radius, 3 ); } } first = false; updatePorosity = false; } max_vel1=3 * width /(height+width+depth)*max_vel; max_vel2=3 * height /(height+width+depth)*max_vel; max_vel3 =3 * depth /(height+width+depth)*max_vel; porosity = ( boxVolume - particlesVolume ) /boxVolume; position_top = p_top->se3.position.y(); position_bottom = p_bottom->se3.position.y(); position_right = p_right->se3.position.x(); position_left = p_left->se3.position.x(); position_front = p_front->se3.position.z(); position_back = p_back->se3.position.z(); // must be done _after_ height, width, depth have been calculated //Update stiffness only if it has been computed by StiffnessCounter (see "stiffnessUpdateInterval") if (scene->iter % stiffnessUpdateInterval == 0 || scene->iter<100) updateStiffness(); bool isARadiusControlIteration = (scene->iter % radiusControlInterval == 0); if (scene->iter % computeStressStrainInterval == 0 || (internalCompaction && isARadiusControlIteration) ) computeStressStrain(); if (!internalCompaction) { Vector3r wallForce (0, goal2*width*depth, 0); if (wall_bottom_activated) { if (stressMask & 2) controlExternalStress(wall_bottom, wallForce, p_bottom, max_vel2); else p_bottom->vel[1] += (-normal[wall_bottom][1]*0.5*goal2*height -p_bottom->vel[1])*(1-strainDamping); /// JH we have to verify //else { p_bottom->vel[1] += (-normal[wall_bottom][1]*0.5*goal2*height -p_bottom->vel[1])*(1-strainDamping); //externalWork += (p_bottom->vel*scene->dt).dot(getForce(scene,wall_id[wall]));} } else p_bottom->vel=Vector3r::Zero(); if (wall_top_activated) { if (stressMask & 2) controlExternalStress(wall_top, -wallForce, p_top, max_vel2); else p_top->vel[1] += (-normal[wall_top][1]*0.5*goal2*height -p_top->vel[1])*(1-strainDamping); } else p_top->vel=Vector3r::Zero(); wallForce = Vector3r(goal1*height*depth, 0, 0); if (wall_left_activated) { if (stressMask & 1) controlExternalStress(wall_left, wallForce, p_left, max_vel1); else p_left->vel[0] += (-normal[wall_left][0]*0.5*goal1*width -p_left->vel[0])*(1-strainDamping); } else p_left->vel=Vector3r::Zero(); if (wall_right_activated) { if (stressMask & 1) controlExternalStress(wall_right, -wallForce, p_right, max_vel1); else p_right->vel[0] += (-normal[wall_right][0]*0.5*goal1*width -p_right->vel[0])*(1-strainDamping); } else p_right->vel=Vector3r::Zero(); wallForce = Vector3r(0, 0, goal3*height*width); if (wall_back_activated) { if (stressMask & 4) controlExternalStress(wall_back, wallForce, p_back, max_vel3); else p_back->vel[2] += (-normal[wall_back][2]*0.5*goal3*depth -p_back->vel[2])*(1-strainDamping); } else p_back->vel=Vector3r::Zero(); if (wall_front_activated) { if (stressMask & 4) controlExternalStress(wall_front, -wallForce, p_front, max_vel3); else p_front->vel[2] += (-normal[wall_front][2]*0.5*goal3*depth -p_front->vel[2])*(1-strainDamping); } else p_front->vel=Vector3r::Zero(); } else //if internal compaction { p_bottom->vel=Vector3r::Zero(); p_top->vel=Vector3r::Zero(); p_left->vel=Vector3r::Zero(); p_right->vel=Vector3r::Zero(); p_back->vel=Vector3r::Zero(); p_front->vel=Vector3r::Zero(); if (isARadiusControlIteration) { Real sigma_iso_ = bool(stressMask & 1)*goal1 + bool(stressMask & 2)*goal2 + bool(stressMask & 4)*goal3; sigma_iso_ /= bool(stressMask & 1) + bool(stressMask & 2) + bool(stressMask & 4); if (std::abs(sigma_iso_)<=std::abs(meanStress)) maxMultiplier = finalMaxMultiplier; if (meanStress==0) previousMultiplier = maxMultiplier; else { // previousMultiplier = 1+0.7*(sigma_iso-s)*(previousMultiplier-1.f)/(s-previousStress); // = (Dsigma/apparentModulus)*0.7 // previousMultiplier = std::max(2-maxMultiplier, std::min(previousMultiplier, maxMultiplier)); if (sigma_iso_ < 0) // compressive case: we have to increase radii if meanStress > sigma_iso_, considering that sigma_iso_ < 0. We end with the same expression as before sign change previousMultiplier = 1.+(sigma_iso_-meanStress)/sigma_iso_*(maxMultiplier-1.); // = (Dsigma/apparentModulus)*0.7 else // tensile case: we have to increase radii if meanStress > sigma_iso_ too. But here sigma_iso_ > 0 => another expression previousMultiplier = 1.+(meanStress-sigma_iso_)/sigma_iso_*(maxMultiplier-1.); // = (Dsigma/apparentModulus)*0.7 } previousStress = meanStress; //Real apparentModulus = (s-previousStress)/(previousMultiplier-1.f); controlInternalStress(previousMultiplier); } } } void TriaxialStressController::computeStressStrain() { scene->forces.sync(); State* p_bottom=Body::byId(wall_bottom_id,scene)->state.get(); State* p_top=Body::byId(wall_top_id,scene)->state.get(); State* p_left=Body::byId(wall_left_id,scene)->state.get(); State* p_right=Body::byId(wall_right_id,scene)->state.get(); State* p_front=Body::byId(wall_front_id,scene)->state.get(); State* p_back=Body::byId(wall_back_id,scene)->state.get(); height = p_top->se3.position.y() - p_bottom->se3.position.y() - thickness; width = p_right->se3.position.x() - p_left->se3.position.x() - thickness; depth = p_front->se3.position.z() - p_back->se3.position.z() - thickness; meanStress = 0; if (height0 == 0) height0 = height; if (width0 == 0) width0 = width; if (depth0 == 0) depth0 = depth; strain[0] = log(width/width0); // all strain values are positiv for extension strain[1] = log(height/height0); strain[2] = log(depth/depth0); volumetricStrain=strain[0]+strain[1]+strain[2]; Real invXSurface = 1.f/(height*depth); Real invYSurface = 1.f/(width*depth); Real invZSurface = 1.f/(width*height); force[wall_bottom]=getForce(scene,wall_id[wall_bottom]); stress[wall_bottom]=force[wall_bottom]*invYSurface; // all stress values are positiv for tension force[wall_top]= getForce(scene,wall_id[wall_top]); stress[wall_top]=-force[wall_top]*invYSurface; force[wall_left]= getForce(scene,wall_id[wall_left]); stress[wall_left]=force[wall_left]*invXSurface; force[wall_right]= getForce(scene,wall_id[wall_right]); stress[wall_right]= -force[wall_right]*invXSurface; force[wall_front]= getForce(scene,wall_id[wall_front]); stress[wall_front]=-force[wall_front]*invZSurface; force[wall_back]= getForce(scene,wall_id[wall_back]); stress[wall_back]= force[wall_back]*invZSurface; for (int i=0; i<6; i++) meanStress+=stress[i].dot(pow(-1.0,i)*normal[i]); // normal[i] is always inwards meanStress/=6.; // ( sXX(xLeft) + sXX(xRight) + sYY(yBottom) + sYY(yTop) + sZZ(zBack) + sZZ(zFront) ) / 6 } void TriaxialStressController::controlInternalStress ( Real multiplier ) { particlesVolume *= pow ( multiplier,3 ); Shop::growParticles(multiplier,true,true); } /*! \fn TriaxialStressController::ComputeUnbalancedForce( bool maxUnbalanced) */ Real TriaxialStressController::ComputeUnbalancedForce( bool maxUnbalanced) {return Shop::unbalancedForce(maxUnbalanced,scene);} -- You received this question notification because your team yade-users is an answer contact for Yade. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp