yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2956
author: Christian Jakob
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
branch nick: yade
timestamp: Thu 2011-11-10 12:04:14 +0100
message:
  Add forgotten files
added:
  examples/CapillaryPhys-example.py
  pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.cpp
  pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp


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

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

=== added file 'examples/CapillaryPhys-example.py'
--- examples/CapillaryPhys-example.py	1970-01-01 00:00:00 +0000
+++ examples/CapillaryPhys-example.py	2011-11-10 11:04:14 +0000
@@ -0,0 +1,72 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+"""
+To run this script you need to have all 10 text files from https://yade-dem.org/wiki/CapillaryTriaxialTest
+in the same folder as this script!
+
+This script shows how to use Law2_ScGeom_CapillaryPhys_Capillarity. The user can switch between hertz and 
+linear model by setting "model_type"
+"""
+
+model_type		= 1	#1=Hertz model with capillary forces, else linear model with capillary model
+
+#some parameters:
+shear_modulus 	= 1e5
+poisson_ratio 	= 0.3
+young_modulus	= 2*shear_modulus*(1+poisson_ratio)
+friction		= 0.5
+angle			= atan(friction)
+local_damping 	= 0.01
+viscous_normal	= 0.021
+viscous_shear	= 0.8*viscous_normal
+lowercorner		= Vector3(0,0,0)
+uppercorner		= Vector3(0.002,0.002,0.004)
+
+#creating a material (FrictMat):
+id_SphereMat=O.materials.append(FrictMat(young=young_modulus,poisson=poisson_ratio,density=2500,frictionAngle=angle))
+SphereMat=O.materials[id_SphereMat]
+
+#generate particles:
+sp=pack.SpherePack()
+sp.makeCloud(lowercorner,uppercorner,.0002,rRelFuzz=.3)
+O.bodies.append([utils.sphere(c,r,material=SphereMat) for c,r in sp])
+
+#generate boundary:
+O.bodies.append(geom.facetBox(uppercorner/2,uppercorner/2,wire=True,fixed=True,material=SphereMat))
+
+#define engines:
+if model_type == 1:#hertz model with capillary forces
+	O.engines=[
+		ForceResetter(),
+		InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
+		InteractionLoop(
+			[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
+			[Ip2_FrictMat_FrictMat_MindlinCapillaryPhys(label='ContactModel')],#for hertz model only
+			[Law2_ScGeom_MindlinPhys_Mindlin()]#for hertz model only
+		),
+		Law2_ScGeom_CapillaryPhys_Capillarity(CapillaryPressure=10000,hertzOn=True),#for hertz model only
+		GravityEngine(gravity=(0,0,-9.81)),
+		NewtonIntegrator(damping=local_damping),
+	]
+	ContactModel.betan=viscous_normal
+	ContactModel.betas=viscous_shear
+	ContactModel.useDamping=True
+else:
+	O.engines=[
+		ForceResetter(),
+		InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
+		InteractionLoop(
+			[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
+			[Ip2_FrictMat_FrictMat_CapillaryPhys()],	#for linear model only
+			[Law2_ScGeom_FrictPhys_CundallStrack()],	#for linear model only
+		),
+		Law2_ScGeom_CapillaryPhys_Capillarity(CapillaryPressure=10000,hertzOn=False),#for linear model only
+		GravityEngine(gravity=(0,0,-9.81)),
+		NewtonIntegrator(damping=local_damping),
+	]
+
+#set time step and run simulation:
+O.dt=0.5*utils.PWaveTimeStep()
+O.run(10000,True)
+

=== added file 'pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.cpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.cpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.cpp	2011-11-10 11:04:14 +0000
@@ -0,0 +1,95 @@
+/*************************************************************************
+*  Copyright (C) 2007 by Bruno CHAREYRE                                 *
+*  bruno.chareyre@xxxxxxxxxxx                                        *
+*                                                                        *
+*  This program is free software; it is licensed under the terms of the  *
+*  GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#include"Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp"
+#include<yade/pkg/dem/ScGeom.hpp>
+#include<yade/pkg/dem/HertzMindlin.hpp>
+#include<yade/pkg/common/ElastMat.hpp>
+#include<yade/pkg/common/Dispatching.hpp>
+
+#include<yade/core/Omega.hpp>
+#include<yade/core/Scene.hpp>
+
+YADE_PLUGIN(
+            (MindlinCapillaryPhys)
+            (Ip2_FrictMat_FrictMat_MindlinCapillaryPhys)
+);
+
+MindlinCapillaryPhys::~MindlinCapillaryPhys(){};// destructor
+
+void Ip2_FrictMat_FrictMat_MindlinCapillaryPhys::go( const shared_ptr<Material>& b1 //FrictMat
+					, const shared_ptr<Material>& b2 // FrictMat
+					, const shared_ptr<Interaction>& interaction)
+{
+	if(interaction->phys) return; // no updates of an already existing contact necessary
+
+	shared_ptr<MindlinCapillaryPhys> contactPhysics(new MindlinCapillaryPhys());
+	interaction->phys = contactPhysics;
+
+	FrictMat* mat1 = YADE_CAST<FrictMat*>(b1.get());
+	FrictMat* mat2 = YADE_CAST<FrictMat*>(b2.get());
+	
+	/* from interaction physics */
+	Real Ea = mat1->young;
+	Real Eb = mat2->young;
+	Real Va = mat1->poisson;
+	Real Vb = mat2->poisson;
+	Real fa = mat1->frictionAngle;
+	Real fb = mat2->frictionAngle;
+
+	/* from interaction geometry */
+	GenericSpheresContact* scg = YADE_CAST<GenericSpheresContact*>(interaction->geom.get());		
+	Real Da = scg->refR1>0 ? scg->refR1 : scg->refR2; 
+	Real Db = scg->refR2; 
+	Vector3r normal=scg->normal; 
+
+	/* calculate stiffness coefficients */
+	Real Ga = Ea/(2*(1+Va));
+	Real Gb = Eb/(2*(1+Vb));
+	Real G = (Ga+Gb)/2; // average of shear modulus
+	Real V = (Va+Vb)/2; // average of poisson's ratio
+	Real E = Ea*Eb/((1.-std::pow(Va,2))*Eb+(1.-std::pow(Vb,2))*Ea); // Young modulus
+	Real R = Da*Db/(Da+Db); // equivalent radius
+	Real Rmean = (Da+Db)/2.; // mean radius
+	Real Kno = 4./3.*E*sqrt(R); // coefficient for normal stiffness
+	Real Kso = 2*sqrt(4*R)*G/(2-V); // coefficient for shear stiffness
+	Real frictionAngle = std::min(fa,fb);
+
+	Real Adhesion = 4.*Mathr::PI*R*gamma; // calculate adhesion force as predicted by DMT theory
+
+	/* pass values calculated from above to MindlinCapillaryPhys */
+	contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle); 
+	//mindlinPhys->prevNormal = scg->normal; // used to compute relative rotation
+	contactPhysics->kno = Kno; // this is just a coeff
+	contactPhysics->kso = Kso; // this is just a coeff
+	contactPhysics->adhesionForce = Adhesion;
+	
+	contactPhysics->kr = krot;
+	contactPhysics->ktw = ktwist;
+	contactPhysics->maxBendPl = eta*Rmean; // does this make sense? why do we take Rmean?
+
+	/* compute viscous coefficients */
+	if(en && betan) throw std::invalid_argument("Ip2_FrictMat_FrictMat_MindlinCapillaryPhys: only one of en, betan can be specified.");
+	if(es && betas) throw std::invalid_argument("Ip2_FrictMat_FrictMat_MindlinCapillaryPhys: only one of es, betas can be specified.");
+
+	// en or es specified, just compute alpha, otherwise alpha remains 0
+	if(en || es){
+		Real logE = log((*en)(mat1->id,mat2->id));
+		contactPhysics->alpha = -sqrt(5/6.)*2*logE/sqrt(pow(logE,2)+pow(Mathr::PI,2))*sqrt(2*E*sqrt(R)); // (see Tsuji, 1992)
+	}
+	
+	// betan specified, use that value directly; otherwise give zero
+	else{	
+		contactPhysics->betan=betan ? (*betan)(mat1->id,mat2->id) : 0; 
+		contactPhysics->betas=betas ? (*betas)(mat1->id,mat2->id) : contactPhysics->betan;
+	}
+};
+
+
+
+

=== added file 'pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp	2011-11-10 11:04:14 +0000
@@ -0,0 +1,60 @@
+/*************************************************************************
+*  Copyright (C) 2007 by Bruno CHAREYRE                                  *
+*  bruno.chareyre@xxxxxxxxxxx                                            *
+*                                                                        *
+*  This program is free software; it is licensed under the terms of the  *
+*  GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#pragma once
+
+#include<yade/pkg/common/Dispatching.hpp>
+#include<yade/pkg/common/ElastMat.hpp>
+#include<yade/pkg/dem/HertzMindlin.hpp>
+
+class MindlinCapillaryPhys : public MindlinPhys
+{
+	public :
+		int currentIndexes [4]; // used for faster interpolation (stores previous positions in tables)
+		
+		virtual ~MindlinCapillaryPhys();
+
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR(MindlinCapillaryPhys,MindlinPhys,"Physics (of interaction) for Law2_ScGeom_CapillaryPhys_Capillarity.",
+				((bool,meniscus,false,,"Presence of a meniscus if true"))
+				((Real,CapillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
+				((Real,Vmeniscus,0.,,"Volume of the menicus"))
+				((Real,Delta1,0.,,"Defines the surface area wetted by the meniscus on the smallest grains of radius R1 (R1<R2)"))
+				((Real,Delta2,0.,,"Defines the surface area wetted by the meniscus on the biggest grains of radius R2 (R1<R2)"))
+				((Vector3r,Fcap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
+				((short int,fusionNumber,0.,,"Indicates the number of meniscii that overlap with this one"))
+				,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
+				);
+	REGISTER_CLASS_INDEX(MindlinCapillaryPhys,MindlinPhys);
+};
+REGISTER_SERIALIZABLE(MindlinCapillaryPhys);
+
+
+class Ip2_FrictMat_FrictMat_MindlinCapillaryPhys : public IPhysFunctor
+{
+	public :
+		virtual void go(	const shared_ptr<Material>& b1,
+					const shared_ptr<Material>& b2,
+					const shared_ptr<Interaction>& interaction);
+
+	FUNCTOR2D(FrictMat,FrictMat);
+	YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_MindlinCapillaryPhys,IPhysFunctor, "RelationShips to use with Law2_ScGeom_CapillaryPhys_Capillarity\n\n In these RelationShips all the interaction attributes are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",
+	            ((Real,gamma,0.0,,"Surface energy parameter [J/m^2] per each unit contact surface, to derive DMT formulation from HM"))
+				((Real,eta,0.0,,"Coefficient to determine the plastic bending moment"))
+				((Real,krot,0.0,,"Rotational stiffness for moment contact law"))
+				((Real,ktwist,0.0,,"Torsional stiffness for moment contact law"))
+				((shared_ptr<MatchMaker>,en,,,"Normal coefficient of restitution $e_n$."))
+				((shared_ptr<MatchMaker>,es,,,"Shear coefficient of restitution $e_s$."))
+				((shared_ptr<MatchMaker>,betan,,,"Normal viscous damping coefficient $\\beta_n$."))
+				((shared_ptr<MatchMaker>,betas,,,"Shear viscous damping coefficient $\\beta_s$."))
+		);
+		DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_MindlinCapillaryPhys);
+
+
+