yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 4187
committer: Christian Jakob <jakob@xxxxxxxxxxxxxxxxxxx>
timestamp: Wed 2014-10-08 07:24:38 +0200
message:
  include fluid stiffness for TSC for undrained flow condition
modified:
  pkg/dem/TriaxialStressController.cpp


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'pkg/dem/TriaxialStressController.cpp'
--- pkg/dem/TriaxialStressController.cpp	2014-08-13 13:31:41 +0000
+++ pkg/dem/TriaxialStressController.cpp	2014-10-08 05:24:38 +0000
@@ -17,6 +17,11 @@
 #include<yade/pkg/dem/Shop.hpp>
 #include<yade/core/Clump.hpp>
 
+#ifdef FLOW_ENGINE
+//#include<yade/pkg/pfv/FlowEngine.hpp>
+#include "FlowEngine_FlowEngineT.hpp"
+#endif
+
 CREATE_LOGGER(TriaxialStressController);
 YADE_PLUGIN((TriaxialStressController));
 
@@ -32,8 +37,17 @@
 	);
 }
 
-void TriaxialStressController::updateStiffness ()
-{
+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) 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();
@@ -46,12 +60,19 @@
 			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;
+				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;
+		stiffness[5] += fluidStiffness*width*height/depth;
+	}
 }
 
 void TriaxialStressController::controlExternalStress(int wall, Vector3r resultantForce, State* p, Real wall_max_vel)