sandbox/Antoonvh/cap.c

    Why do plants transpire?

    To pump nutrients through their capillaries:

    Without transpiration, a balance at rest emerges. With evaporation: the capillary pumps fluid in a quasi-steady balance.

    #include "navier-stokes/centered.h"
    #include "contact.h"
    #include "two-phase.h"
    #include "tension.h"
    #include "view.h"
    #define BVIEW 1
    #include "particles.h"

    The contact angle is setup.

    bid capillary;
    vector h[];
    double theta0 = 50.;
    h.t[left] = contact_angle (theta0*pi/180.);
    h.t[capillary] = contact_angle (x < 2 ? theta0*pi/180.: pi/2.);
    
    face vector av[];
    
    int main() {
      L0 = 11.;
      mu1 = mu2 = 0.1;
      a = av;
      rho1 = 2.;
      rho2 = 1.;
      f.sigma = 0.0015;
      f.height = h;
      N = 32;
      run();
    }
    
    event init (t = 0) {
      foreach()
        f[]  = (y < 4);

    Contant angles do work on straight masked cells.

      mask (y > 2 && y < 8 && x > 1 && x < 3 ? capillary : none);

    Nutrients are suspended in the water.

      foreach() 
        if (f[] >= 0.9)
          n_part++;
      loc = (coord*) malloc (sizeof(coord)*n_part);
      n_part = 0;
      foreach()
        if (f[] >= 0.9) {
          coord new = {x, y, z};
          loc[n_part++] = new;
        }
    }

    Gravity is an important part of the column’s force balance.

    event acceleration (i++) {
      foreach_face(y)
        av.y[] = -0.001;
    }

    Evaporation of the fluid in the cappilary is hacked in via the event below. The aim is to implement a fixed evaportation at the surface of the column.

    bool evaporate = false;
    event turn_on_evaporation (t = 2000)
      evaporate = true;
    
    event evaporation (i++) {
      if (evaporate) {
        double evap_rate = 0.01;
        foreach() {
          if (f[] > 1e-5 && f[] < 1. - 1e-5 && x < 2.) {
    	coord n = interface_normal (point, f), p;
    	double alpha = plane_alpha (f[], n);
    	double l = plane_area_center (n, alpha, &p);
    	double flux = dt*l*evap_rate;
    	f[] = max (f[] - flux, 0);
          }
        }
      }
    }
    
    event movie (t += 10; t <= 12500) {
      char str[99];
      scalar white[];
      foreach()
        white[] = 0;
      view (tx = -0.5, ty = -0.5, bg = {0.9, 0.9, 0.9});
      if (evaporate)
        sprintf (str, "With evaporation");
      else
        sprintf (str, "Without evaporation");
      draw_string (str, size = 30);
      draw_vof ("f");
      draw_vof ("f", filled = 1, fc = {0.2, 0.4, 0.9});
      scatter (loc);
      translate (z = -0.1)
        squares ("white", map = cool_warm, min = -1., max = 1);
      save ("movie.mp4");
    }
     
    event adapt (i++) 
      adapt_wavelet ({f}, (double[]){0.001}, 8, 6);