twalls.c - simple_DEM - a simple 2D Discrete Element Method code for educational purposes
 (HTM) git clone git://src.adamsgaard.dk/simple_DEM
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       twalls.c (2507B)
       ---
            1 #include <math.h>
            2 #include "header.h"
            3 #include "global_properties.h"
            4 
            5 /* Compute force between i and upper wall */
            6 double compute_force_upper_wall(int i, grain* g, double wup)
            7 {
            8   double dn = wup - (g[i].y + g[i].R); /* Overlap */
            9   
           10   if(dn<0.0) {
           11     double vn,fn;
           12     /* velocity (wall velocity = 0) */
           13     vn = g[i].vy;
           14     /* force */
           15     fn = -kn * dn - nu * vn;
           16     /* Update sum of forces on grains */
           17     g[i].fy -= fn;
           18     /* Add fn to pressure on grains i */
           19     g[i].p += fn;
           20     /* Update stress to the wall */
           21     return fn;
           22   }
           23   return 0.0;
           24 }
           25 
           26 double compute_force_lower_wall(int i, grain* g, double wdown)
           27 {
           28 
           29   double dn = g[i].y - g[i].R - wdown; /* Overlap */
           30   
           31   if(dn<0.0) {
           32     double vn,fn;
           33     /* velocity (wall velocity = 0) */
           34     vn = g[i].vy;
           35     /* force */
           36     fn = - kn * dn - nu * vn;
           37     /* Update sum of forces on grains */
           38     g[i].fy += fn;
           39     /* Add fn to pressure on grains i */
           40     g[i].p += fn;
           41     /* Update stress to the wall */
           42     return fn;
           43   }
           44   return 0.0;
           45 }
           46 
           47 double compute_force_left_wall(int i, grain* g, double wleft)
           48 {
           49   double dn = wleft + g[i].x - g[i].R; /* Overlap */
           50   
           51   if(dn<0.0) {
           52     double vn,fn;
           53     /* velocity (wall velocity = 0) */
           54     vn = g[i].vx;
           55     /* force */
           56     fn = -kn * dn - nu * vn;
           57     /* Update sum of forces on grains */
           58     g[i].fx += fn;
           59     /* Add fn to pressure on grains i */
           60     g[i].p += fn;
           61     /* Update stress to the wall */
           62     return fn;
           63   }
           64   return 0.0;
           65 }
           66 
           67 
           68 double compute_force_right_wall(int i, grain* g, double wright)
           69 {
           70   double dn = wright - (g[i].x + g[i].R); /* Overlap */
           71   
           72   if(dn<0.0) {
           73     double vn,fn;
           74     /* velocity (wall velocity = 0) */
           75     vn = -g[i].vx;
           76     /* force */
           77     fn = -kn * dn - nu * vn;
           78     /* Update sum of forces on grains */
           79     g[i].fx -= fn;
           80     /* Add fn to pressure on grains i */
           81     g[i].p += fn;
           82     /* Update stress to the wall */
           83     return fn;
           84   }
           85   return 0.0;
           86 }
           87 
           88 void interact_walls(grain* g, double wleft, double wright, double wup, double wdown, 
           89                         double* wp_up, double* wp_down, double* wp_left, double* wp_right)
           90 {
           91   double u = 0.0;
           92   double d = 0.0;
           93   double r = 0.0;
           94   double l = 0.0;
           95 
           96   int i;
           97 
           98   #pragma omp parallel for shared(g, u, d, r, l) private (i)
           99   for (i = 0; i < np; i++) {
          100     u += compute_force_upper_wall(i, g, wup);
          101     d += compute_force_lower_wall(i, g, wdown);
          102     r += compute_force_right_wall(i, g, wright);
          103     l += compute_force_left_wall(i, g, wleft);
          104   }
          105   *wp_up         = u;
          106   *wp_down         = d;
          107   *wp_left        = l;
          108   *wp_right        = r;
          109 }