yade-users team mailing list archive

[Luc Scholtes <scholtes@xxxxxxxxxxxxxxx>]

Hi,

I would like to introduce an interaction force which can occur without
any contact between 2 discrete elements in order to simulate capillary
phenomenon in unsaturated granular media. I thought increasing the size 
of bouding boxes a little bit would cause neighbooring (but without 
contact) particles to be included in the list of volatile interctions, 
so that I could iterate over this "augmented" list of interaction and 
retrieve particles on which cappilary forces apply. I have created the
corresponding files, but there is an aspect that I can't overcome: the
resultant force is not effective if the particules are not in contact,
even if their bounding volume interact. This is certainly due to my
definition of the interactions...

I join the files. If you have any suggestion...

hpp :

   #ifndef CAPILLARY_COHESIVE_LAW_HPP
   #define CAPILLARY_COHESIVE_LAW_HPP

   #include <yade/yade-core/InteractionSolver.hpp>

   #include <set>
   #include <boost/tuple/tuple.hpp>

   class PhysicalAction;

   class CapillaryCohesiveLaw : public InteractionSolver
   {
       private :
           shared_ptr<PhysicalAction> actionForce;
           shared_ptr<PhysicalAction> actionMomentum;

       public :
           int sdecGroupMask;
           bool momentRotationLaw;

           CapillaryCohesiveLaw();
           void action(Body* body);

       protected :
           void registerAttributes();
       REGISTER_CLASS_NAME(CapillaryCohesiveLaw);
       REGISTER_BASE_CLASS_NAME(InteractionSolver);

   };

   REGISTER_SERIALIZABLE(CapillaryCohesiveLaw,false);

   #endif // CAPILLARY_COHESIVE_LAW_HPP


cpp :

   #include "CapillaryCohesiveLaw.hpp"
   #include "BodyMacroParameters.hpp"
   #include "SpheresContactGeometry.hpp"
   #include "ElasticContactParameters.hpp"
   #include <yade/yade-core/Omega.hpp>
   #include <yade/yade-core/MetaBody.hpp>
   #include <yade/yade-package-common/Force.hpp>
   #include <yade/yade-package-common/Momentum.hpp>
   #include <yade/yade-core/PhysicalAction.hpp>
   #include <yade/yade-lib-wm3-math/Math.hpp>


   //using namespace Math;

   CapillaryCohesiveLaw::CapillaryCohesiveLaw() : InteractionSolver() ,
   actionForce(new Force) , actionMomentum(new Momentum)
   {
       sdecGroupMask=1;
       momentRotationLaw = true;
   }


   void CapillaryCohesiveLaw::registerAttributes()
   {
       InteractionSolver::registerAttributes();
       REGISTER_ATTRIBUTE(sdecGroupMask);
       REGISTER_ATTRIBUTE(momentRotationLaw);
   }

   void CapillaryCohesiveLaw::action(Body* body)
   {
       MetaBody * ncb = dynamic_cast<MetaBody*>(body);
       shared_ptr<BodyContainer>& bodies = ncb->bodies;


       InteractionContainer::iterator ii =
   ncb->volatileInteractions->begin();
       InteractionContainer::iterator iiEnd =
   ncb->volatileInteractions->end();
       // persistent interactions?

       for(  ; ii!=iiEnd ; ++ii )
       {

           if ((*ii)->isReal)   // pairs of BoundingBox in interaction?
           {
               const shared_ptr<Interaction>& interaction = *ii;
               unsigned int id1 = interaction->getId1();
               unsigned int id2 = interaction->getId2();

               if( !( (*bodies)[id1]->getGroupMask() &
   (*bodies)[id2]->getGroupMask() & sdecGroupMask)  )
                   continue;


               BodyMacroParameters* de1         =
dynamic_cast<BodyMacroParameters*>((*bodies)[id1]->physicalParameters.get());

               BodyMacroParameters* de2         =
dynamic_cast<BodyMacroParameters*>((*bodies)[id2]->physicalParameters.get());

               SpheresContactGeometry* currentContactGeometry     =

dynamic_cast<SpheresContactGeometry*>(interaction->interactionGeometry.get());
                // SDECLinkGeometry?

               ElasticContactParameters* currentContactPhysics       =

dynamic_cast<ElasticContactParameters*>(interaction->interactionPhysics.get());
               // SDECLinkPhysics?

               // intergranular distance
               Real un =
(de2->se3.position-de1->se3.position).length()-(currentContactGeometry->radius1+currentContactGeometry->radius2);



               if (un <
2*std::max(currentContactGeometry->radius2,currentContactGeometry->radius1))

               {
                   // capillary force (test)
                   Vector3r Fcap=100000*currentContactGeometry->normal;

                   static_cast<Force*>   ( ncb->actionParameters->find(
   id1 , actionForce   ->getClassIndex() ).get() )->force    += Fcap;
                   static_cast<Force*>   ( ncb->actionParameters->find(
   id2 , actionForce   ->getClassIndex() ).get() )->force    -= Fcap;

                   currentContactPhysics->prevNormal =
   currentContactGeometry->normal;
               }
           }
       }

thanks in advance, so long.



_______________________________________________
Yade-users mailing list
Yade-users@xxxxxxxxxxxxxxxx
http://lists.berlios.de/mailman/listinfo/yade-users