← Back to branch summary

~ubuntu-branches/ubuntu/trusty/yade/trusty

« back to all changes in this revision

Viewing changes to pkg/dem/ViscoelasticCapillarPM.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky, cf3f8d9
  • Date: 2013-10-30 20:56:33 UTC
  • mfrom: (20.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20131030205633-1f01r7hjce17d723
Tags: 1.05.0-2
http://bugs.debian.org/726009
[cf3f8d9] Pass -ftrack-macro-expansion=0 only if gcc>=4.8. (Closes: #726009)

Show diffs side-by-side

added added

removed removed

Lines of Context:
pkg/dem/ViscoelasticCapillarPM.cpp
 
1
#include"ViscoelasticPM.hpp"
 
2
#include<yade/core/State.hpp>
 
3
#include<yade/pkg/dem/ScGeom.hpp>
 
4
#include<yade/core/Omega.hpp>
 
5
#include<yade/core/Scene.hpp>
 
6
#include<yade/pkg/common/Sphere.hpp>
 
7
 
 
8
Real Law2_ScGeom_ViscElPhys_Basic::calculateCapillarForce(const ScGeom& geom, ViscElPhys& phys) {
 
9
  Real fC = 0.0;
 
10
  
 
11
  /* Capillar
 
12
    * Some equations have constants, which can be calculated only once per contact. 
 
13
    * No need to recalculate them at each step. 
 
14
    * It needs to be fixed.
 
15
    * 
 
16
    */
 
17
     
 
18
  if (phys.CapillarType  == "Weigert") {
 
19
      /* Capillar model from [Weigert1999]
 
20
       */
 
21
        Real R = phys.R;
 
22
        Real a = -geom.penetrationDepth;
 
23
        Real Ca = (1.0 + 6.0*a/(R*2.0));                                                          // [Weigert1999], equation (16)
 
24
        Real Ct = (1.0 + 1.1*sin(phys.theta));                                                    // [Weigert1999], equation (17)
 
25
        
 
26
        /*
 
27
        Real Eps = 0.36;                                                                          // Porosity
 
28
        Real fi = phys.Vb/(2.0*M_PI/6.0*pow(R*2.0,3.));                                           // [Weigert1999], equation (13)
 
29
        Real S = M_PI*(1-Eps)/(pow(Eps, 2.0))*fi;                                                 // [Weigert1999], equation (14)
 
30
        Real beta = asin(pow(((S/0.36)*(pow(Eps, 2.0)/(1-Eps))*(1.0/Ca)*(1.0/Ct)), 1.0/4.0));     // [Weigert1999], equation (19)
 
31
        */
 
32
        
 
33
        Real beta = asin(pow(phys.Vb/(0.12*Ca*Ct*pow(2.0*R, 3.0)), 1.0/4.0));                     // [Weigert1999], equation (15), against Vb
 
34
        
 
35
        Real r1 = (2.0*R*(1-cos(beta)) + a)/(2.0*cos(beta+phys.theta));                           // [Weigert1999], equation (5)
 
36
        Real r2 = R*sin(beta) + r1*(sin(beta+phys.theta)-1);                                      // [Weigert1999], equation (6)
 
37
        Real Pk = phys.gamma*(1/r1 - 1/r2);                                                       /* [Weigert1999], equation (22),
 
38
                                                                                                   * see also a sentence over the equation
 
39
                                                                                                   * "R1 was taken as positive and R2 was taken as negative"
 
40
                                                                                                   */ 
 
41
 
 
42
        //fC = M_PI*2.0*R*phys.gamma/(1+tan(0.5*beta));                                           // [Weigert1999], equation (23), [Fisher]
 
43
        
 
44
        fC = M_PI/4.0*pow((2.0*R),2.0)*pow(sin(beta),2.0)*Pk +
 
45
             phys.gamma*M_PI*2.0*R*sin(beta)*sin(beta+phys.theta);                                // [Weigert1999], equation (21)
 
46
        
 
47
      } else if (phys.CapillarType  == "Willett_numeric") {
 
48
      
 
49
        /* Capillar model from [Willett2000]
 
50
         */ 
 
51
        
 
52
        Real R = phys.R;
 
53
        Real s = -geom.penetrationDepth;
 
54
        Real Vb = phys.Vb;
 
55
        
 
56
        Real VbS = Vb/(R*R*R);
 
57
        Real Th1 = phys.theta;
 
58
        Real Th2 = phys.theta*phys.theta;
 
59
        Real Gamma = phys.gamma;
 
60
        
 
61
        /*
 
62
         * [Willett2000], equations in Anhang
 
63
        */
 
64
        Real f1 = (-0.44507 + 0.050832*Th1 - 1.1466*Th2) + 
 
65
                  (-0.1119 - 0.000411*Th1 - 0.1490*Th2) * log(VbS) +
 
66
                  (-0.012101 - 0.0036456*Th1 - 0.01255*Th2) *log(VbS)*log(VbS) +
 
67
                  (-0.0005 - 0.0003505*Th1 - 0.00029076*Th2) *log(VbS)*log(VbS)*log(VbS);
 
68
        
 
69
        Real f2 = (1.9222 - 0.57473*Th1 - 1.2918*Th2) +
 
70
                  (-0.0668 - 0.1201*Th1 - 0.22574*Th2) * log(VbS) +
 
71
                  (-0.0013375 - 0.0068988*Th1 - 0.01137*Th2) *log(VbS)*log(VbS);
 
72
                  
 
73
        Real f3 = (1.268 - 0.01396*Th1 - 0.23566*Th2) +
 
74
                  (0.198 + 0.092*Th1 - 0.06418*Th2) * log(VbS) +
 
75
                  (0.02232 + 0.02238*Th1 - 0.009853*Th2) *log(VbS)*log(VbS) +
 
76
                  (0.0008585 + 0.001318*Th1 - 0.00053*Th2) *log(VbS)*log(VbS)*log(VbS);
 
77
        
 
78
        Real f4 = (-0.010703 + 0.073776*Th1 - 0.34742*Th2) +
 
79
                  (0.03345 + 0.04543*Th1 - 0.09056*Th2) * log(VbS) +
 
80
                  (0.0018574 + 0.004456*Th1 - 0.006257*Th2) *log(VbS)*log(VbS);
 
81
  
 
82
        Real sPl = (s/2.0)/sqrt(Vb/R);
 
83
        
 
84
        Real lnFS = f1 - f2*exp(f3*log(sPl) + f4*log(sPl)*log(sPl));
 
85
        Real FS = exp(lnFS);
 
86
        
 
87
        fC = FS * 2.0 * M_PI* R * Gamma;
 
88
      } else if (phys.CapillarType  == "Willett_analytic") {
 
89
        /* Capillar model from Willet [Willett2000] (analytical solution), but 
 
90
         * used also in the work of Herminghaus [Herminghaus2005]
 
91
         */
 
92
         
 
93
        Real R = phys.R;
 
94
        Real Gamma = phys.gamma;
 
95
        Real s = -geom.penetrationDepth;
 
96
        Real Vb = phys.Vb;
 
97
                
 
98
        /*
 
99
        
 
100
        Real sPl = s/sqrt(Vb/R);                                                            // [Herminghaus2005], equation (sentence between (7) and (8))
 
101
        fC = 2.0 * M_PI* R * Gamma * cos(phys.theta)/(1 + 1.05*sPl + 2.5 *sPl * sPl);       // [Herminghaus2005], equation (7)
 
102
        
 
103
        */ 
 
104
        
 
105
        Real sPl = (s/2.0)/sqrt(Vb/R);                                                      // [Willett2000], equation (sentence after (11)), s - half-separation, so s*2.0
 
106
        Real f_star = cos(phys.theta)/(1 + 2.1*sPl + 10.0 * pow(sPl, 2.0));                 // [Willett2000], equation (12)
 
107
        fC = f_star * (2*M_PI*R*Gamma);                                                     // [Willett2000], equation (13), against F
 
108
        
 
109
      } else if (phys.CapillarType  == "Rabinovich") {
 
110
        /* Capillar model from Rabinovich [Rabinov2005]
 
111
         */
 
112
         
 
113
        Real R = phys.R;
 
114
        Real Gamma = phys.gamma;
 
115
        Real H = -geom.penetrationDepth;
 
116
        Real V = phys.Vb;
 
117
        
 
118
        Real alpha = 0.0;
 
119
        Real dsp = 0.0;
 
120
        if (H!=0.0) {
 
121
          alpha = sqrt(H/R*(-1+ sqrt(1 + 2.0*V/(M_PI*R*H*H))));                          // [Rabinov2005], equation (A3)
 
122
        
 
123
          dsp = H/2.0*(-1.0 + sqrt(1.0 + 2.0*V/(M_PI*R*H*H)));                           // [Rabinov2005], equation (20)
 
124
        
 
125
          fC = -(2*M_PI*R*Gamma*cos(phys.theta))/(1+(H/(2*dsp))) - 
 
126
              2*M_PI*R*Gamma*sin(alpha)*sin(phys.theta + alpha);                         // [Rabinov2005], equation (19)
 
127
        } else {
 
128
          fC = -(2*M_PI*R*Gamma*cos(phys.theta)) - 
 
129
              2*M_PI*R*Gamma*sin(alpha)*sin(phys.theta + alpha);                         // [Rabinov2005], equation (19)
 
130
        }
 
131
        
 
132
        fC *=-1;
 
133
        
 
134
      } else {
 
135
        throw runtime_error("CapillarType is unknown, please, use only Willett_numeric, Willett_analytic, Weigert or Rabinovich");
 
136
      }
 
137
  return fC;
 
138
}