yade-users team mailing list archive

[Yor1 <question290129@xxxxxxxxxxxxxxxxxxxxx>]

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.