twocylinders_v4.fe

// Angle between cylinder normals, should be between 0 and pi/2
parameter cyl_angle_degree = 30
#define cyl_angle (cyl_angle_degree*pi/180)
// Length of cylinders in 3d space
parameter cyl_len = 5
// Radius of the upper, lower and droplet cylinders
parameter rad_upper = 0.45
parameter rad_lower = 0.5
// This radius must be smaller than rad_lower and rad_upper. A workable value is either the maximum or the average of these.
parameter rad_center = 0.3
// Volume of the droplet
parameter VOL = 0.08
// Step size for finding detachment from decreasing volume. Smaller is better, but takes considerably longer. Recommend lowering VOL until detachment, and then using this for finer approximation.
parameter volstep = 0.0005
// Contact angle, in degrees because why not
parameter cont_angle_up = 60
parameter cont_angle_low = 120
// Some reference values
PARAMETER min_area = 0.002
PARAMETER min_lent = 0.05
PARAMETER ref_length = 10*min_lent
PARAMETER ref_area = 25*min_area
// Strength of gravity
PARAMETER GZ = -980
// The height should be greater than rad_lower+rad_upper
parameter height = rad_upper + rad_lower + 0.5
// density of the fluid
#define DENS 1
// the surface tension of the liquid-air interface
PARAMETER TENS = 72
// surface tension of liquid-solid interface, uses contact angle
#define tension_val_up (cos(cont_angle_up*pi/180)*TENS)
#define tension_val_low (cos(cont_angle_low*pi/180)*TENS)

// Lower Cylinder - display constraint, lies on cylinder along x-axis with radius rad_lower
constraint low_disp_con
formula: sqrt(y^2 + (height+z)^2) = rad_lower

// Upper Cylinder - display constraint, lies on cylinder along (cos(cyl_angle),sin(cyl_angle)) with radius rad_upper
constraint up_disp_con
formula: sqrt(z^2 + (y*cos(cyl_angle) - x*sin(cyl_angle))^2) = rad_upper

// Keep upper cylinder ends flat, lies on two planes in space along cylinder normal
constraint up_disp_con_plane_back
formula: x*cos(cyl_angle) + y*sin(cyl_angle) = cyl_len

constraint up_disp_con_plane_front
formula: x*cos(cyl_angle) + y*sin(cyl_angle) = -cyl_len

// Define common/long terms for the upper energy
#define denom_u (z^2 + (y*cos(cyl_angle) - x*sin(cyl_angle))^2) // upper cylinder denominator
#define common_u ((tension_val_up*rad_upper/2)/denom_u) // common front of w for upper

#define uproot (sqrt(rad_upper^2 - (y*cos(cyl_angle) - x*sin(cyl_angle))^2))
// Upper cylinder constraint with energy. This defines the upper cylinder's liquid solid interface.
constraint up_energy
formula: sqrt(z^2 + (y*cos(cyl_angle) - x*sin(cyl_angle))^2) = rad_upper
// Describes contact angle term (first term) and gravitational (second term, if there is one). Each component e1,e2,e3 forms a vector integrand, integrated over the boundary.
energy:
e1: -common_u*(-2*y*z*(sin(cyl_angle))^2 - x*z*sin(2*cyl_angle)) -(GZ*DENS*y/2)*(y^2 * (cos(cyl_angle))^2 /3 + (x*sin(cyl_angle))^2)
e2: -common_u*(2*x*z*(cos(cyl_angle))^2 + y*z*sin(2*cyl_angle)) +(GZ*DENS*x/2)*(rad_upper^2 + x*y*sin(2*cyl_angle)/2)
e3: -common_u*(-2*x*y*cos(2*cyl_angle)+ (x^2 - y^2)*sin(2*cyl_angle))
// volume content integral, doesn't work with cyl_angle=90. Can simply use 89.9, or 90.1.
content:
c1: -0.5*((rad_upper^2 / cos(cyl_angle)) * atan((y*cos(cyl_angle) - x*sin(cyl_angle))/uproot) + (y-x*tan(cyl_angle))*uproot)
c2: 0
c3: 0


// Lower cylinder with energy. This defines the lower cylinder's liquid solid interface.
constraint low_energy
formula: sqrt(y^2 + (height+z)^2)= rad_lower
// Describes contact angle term (first term) and gravitational (second term, if there is one). Each component e1,e2,e3 forms a vector integrand, integrated over the boundary.
energy:
e1: 0
e2: -tension_val_low*rad_lower*x*(z+height) / (y^2+(z+height)^2) -(GZ*DENS*x/2)*(-height + sqrt(rad_lower^2 - y^2))
e3: tension_val_low*rad_lower*x*y / (y^2+(z+height)^2)
// volume content integral
content:
c1: 0
c2: height*x - x*sqrt(rad_lower^2 - y^2)
c3: 0

// named quantity for arbitrary direction gravity on facets, calculates the gravitational effects on the liquid-air interface. Refer to Surface Evolver manual for this wizardry.
quantity arb_grav energy method facet_vector_integral global
vector_integrand:
q1: 0
q2: 0
q3: GZ*DENS*(z^2) / 2
//
// The following is a list of vertices, edges and faces. Refer to the vertex images on GitHub. The relevant constraints are applied where needed.
// useful parameter, bad to calculate multiple times
PARAMETER tan_num =tan(cyl_angle/2)
vertices
// Lower cylinder display - defines geometry
1 -cyl_len 0 -rad_lower-height constraint low_disp_con FIXED
2 -cyl_len rad_lower -height constraint low_disp_con FIXED
3 -cyl_len 0 rad_lower-height constraint low_disp_con FIXED
4 -cyl_len -rad_lower -height constraint low_disp_con FIXED
5  cyl_len -rad_lower -height constraint low_disp_con FIXED
6  cyl_len 0 rad_lower-height constraint low_disp_con FIXED
7  cyl_len rad_lower -height constraint low_disp_con FIXED
8  cyl_len 0 -rad_lower-height constraint low_disp_con FIXED

// Upper cylinder display - defines geometry with horrendously long formulae.
9 -cyl_len*cos(cyl_angle) -cyl_len*sin(cyl_angle) -rad_upper constraint up_disp_con up_disp_con_plane_back FIXED
10 (-(-cyl_len)*tan_num^2 - 2*(rad_upper)*tan_num + (-cyl_len))/(tan_num^2 + 1) (2*(-cyl_len)*tan_num - (rad_upper)*tan_num^2 + (rad_upper))/(tan_num^2 + 1) 0 constraint up_disp_con up_disp_con_plane_back FIXED
11 -cyl_len*cos(cyl_angle) -cyl_len*sin(cyl_angle) rad_upper constraint up_disp_con up_disp_con_plane_back FIXED
12 (-(-cyl_len)*tan_num^2 - 2*(-rad_upper)*tan_num + (-cyl_len))/(tan_num^2 + 1) (2*(-cyl_len)*tan_num - (-rad_upper)*tan_num^2 + (-rad_upper))/(tan_num^2 + 1) 0 constraint up_disp_con up_disp_con_plane_back FIXED
13 (-(cyl_len)*tan_num^2 - 2*(-rad_upper)*tan_num + (cyl_len))/(tan_num^2 + 1) (2*(cyl_len)*tan_num - (-rad_upper)*tan_num^2 + (-rad_upper))/(tan_num^2 + 1) 0 constraint up_disp_con up_disp_con_plane_front FIXED
14 cyl_len*cos(cyl_angle) cyl_len*sin(cyl_angle) rad_upper constraint up_disp_con up_disp_con_plane_front FIXED
15 (-(cyl_len)*tan_num^2 - 2*(rad_upper)*tan_num + (cyl_len))/(tan_num^2 + 1) (2*(cyl_len)*tan_num - (rad_upper)*tan_num^2 + (rad_upper))/(tan_num^2 + 1) 0 constraint up_disp_con up_disp_con_plane_front FIXED
16 cyl_len*cos(cyl_angle) cyl_len*sin(cyl_angle) -rad_upper constraint up_disp_con up_disp_con_plane_front FIXED

// several useful parameters to reduce the length of the formulae
#define RC (rad_center)
#define RU (rad_upper)
#define RL (rad_lower)
#define h (height)
// two parameters that show up repeatedly, first is a corner, second describes the y component
#define z_low_corn (-h+sqrt(RL^2 - (RC^2) / 2))
#define z_low_y (-h+sqrt(RL^2 - RC^2))
// Droplet initial geometry
17 -RC 0 -h+RL constraint low_energy
18 -RC/sqrt(2) RC/sqrt(2) z_low_corn constraint low_energy
19 0 RC z_low_y constraint low_energy
20 RC/sqrt(2) RC/sqrt(2) z_low_corn constraint low_energy
21 RC 0 -h+RL constraint low_energy
22 RC/sqrt(2) -RC/sqrt(2) z_low_corn constraint low_energy
23 0 -RC z_low_y constraint low_energy
24 -RC/sqrt(2) -RC/sqrt(2) z_low_corn constraint low_energy
25 -RC 0 -sqrt(RU^2 - ((0)*cos(cyl_angle) -(-RC)*sin(cyl_angle))^2) constraint up_energy
26 -RC/sqrt(2) RC/sqrt(2) -sqrt(RU^2 - ((RC/sqrt(2))*cos(cyl_angle) -(-RC/sqrt(2))*sin(cyl_angle))^2) constraint up_energy
27 0 RC -sqrt(RU^2 - ((RC)*cos(cyl_angle) -(0)*sin(cyl_angle))^2) constraint up_energy
28 RC/sqrt(2) RC/sqrt(2) -sqrt(RU^2 - ((RC/sqrt(2))*cos(cyl_angle) -(RC/sqrt(2))*sin(cyl_angle))^2) constraint up_energy
29 RC 0 -sqrt(RU^2 - ((0)*cos(cyl_angle) -(RC)*sin(cyl_angle))^2) constraint up_energy
30 RC/sqrt(2) -RC/sqrt(2) -sqrt(RU^2 - ((-RC/sqrt(2))*cos(cyl_angle) -(RC/sqrt(2))*sin(cyl_angle))^2) constraint up_energy
31 0 -RC -sqrt(RU^2 - ((-RC)*cos(cyl_angle) -(0)*sin(cyl_angle))^2) constraint up_energy
32 -RC/sqrt(2) -RC/sqrt(2) -sqrt(RU^2 - ((-RC/sqrt(2))*cos(cyl_angle) -(-RC/sqrt(2))*sin(cyl_angle))^2) constraint up_energy

edges
// Lower cylinder display - fixed, as not required to evolve
1  1 2 constraint low_disp_con FIXED
2  2 3 constraint low_disp_con FIXED
3  3 4 constraint low_disp_con FIXED
4  4 1 constraint low_disp_con FIXED
5  1 8 constraint low_disp_con FIXED
6  2 7 constraint low_disp_con FIXED
7  3 6 constraint low_disp_con FIXED
8  4 5 constraint low_disp_con FIXED
9  8 7 constraint low_disp_con FIXED
10 7 6 constraint low_disp_con FIXED
11 6 5 constraint low_disp_con FIXED
12 5 8 constraint low_disp_con FIXED

// Upper cylinder display - fixed, as not required to evolve
13 9 10 constraint up_disp_con up_disp_con_plane_back FIXED
14 10 11 constraint up_disp_con up_disp_con_plane_back FIXED
15 11 12 constraint up_disp_con up_disp_con_plane_back FIXED
16 12 9 constraint up_disp_con up_disp_con_plane_back FIXED
17 9 16 constraint up_disp_con FIXED
18 16 15 constraint up_disp_con FIXED
19 15 10 constraint up_disp_con FIXED
20 15 14 constraint up_disp_con FIXED
21 13 14 constraint up_disp_con up_disp_con_plane_front FIXED
22 16 13 constraint up_disp_con up_disp_con_plane_front FIXED
23 12 13 constraint up_disp_con up_disp_con_plane_front FIXED
24 14 11 constraint up_disp_con up_disp_con_plane_front FIXED

// Droplet - edges on the cylinder boundaries have energy constraints
25 17 24 constraint low_energy
26 24 23 constraint low_energy
27 23 22 constraint low_energy
28 22 21 constraint low_energy
29 21 20 constraint low_energy
30 20 19 constraint low_energy
31 19 18 constraint low_energy
32 18 17 constraint low_energy
33 25 17
34 32 25 constraint up_energy
35 24 32
36 31 32 constraint up_energy
37 31 23
38 30 31 constraint up_energy
39 22 30
40 29 30 constraint up_energy
41 29 21
42 28 29 constraint up_energy
43 20 28
44 27 28 constraint up_energy
45 27 19
46 26 27 constraint up_energy
47 18 26
48 25 26 constraint up_energy

faces
// Lower cylinder display
1 9 -6 -1 5 constraint low_disp_con density 0 FIXED// right bottom
2 10 -7 -2 6 constraint low_disp_con density 0 FIXED// right top
3 -3 7 11 -8 constraint low_disp_con density 0 FIXED// left top
4 8 12 -5 -4 constraint low_disp_con density 0 FIXED// left bottom

// Upper cylinder display
5 20 24 -14 -19 constraint up_disp_con density 0 FIXED// right top
6 -15 -24 -21 -23 constraint up_disp_con density 0 FIXED// left top
7 23 -22 -17 -16 constraint up_disp_con density 0 FIXED// left bot
8 19 -13 17 18 constraint up_disp_con density 0 FIXED// right bot

// Droplet - has density TENS to deal with liquid-air surface energy
9 47 -48 33 -32 density TENS
10 -45 -46 -47 -31 density TENS
11 43 -44 45 -30 density TENS
12 -41 -42 -43 -29 density TENS
13 41 -28 39 -40 density TENS
14 -39 -27 -37 -38 density TENS
15 35 -36 37 -26 density TENS
16 -33 -34 -35 -25 density TENS

// Lets SE know this is a fluid.
bodies
1 -9 -10 -11 -12 -13 -14 -15 -16 volume VOL density DENS

// COMMANDS:
// Other script commands to use during runtime
// Type these into the surface evolver console
read

// change from linear mode to quadratic

// Making nicely refined cylinders and then freezing them so they don't refine more.
ref_cylinders := {
  set edge no_refine where not fixed;
  r; t 0.1; u; r;
  set facet no_refine where fixed;
  unset edge no_refine;
  set edge no_refine where fixed;
}
// Adds colours to cylinders
do_cols := {
foreach facet ff where on_constraint low_disp_con do set ff color red;
foreach facet ff where on_constraint up_disp_con do set ff color red;	
set body[1].facets color cyan;
}
// run "calc" to calculate the minimal surface - You should pretty much always run this first
calc := {
G 0;
ref_cylinders;
do_cols;
r; {g 50; V; u; } 5; M 2;{g 10; V; u; } 5; r; {g 10; V; u; } 5; u; {g 10; V; u; } 5; U; {g 10; V; u; } 5; U; {g 10; V; u; } 5;
}
// dead command
gogo := {
//grow droplet:
{
g 400;
};
};
// useful command to evolve mesh without human input
gogoconjugate := {{g 200; U; g 200; U; V 5} 3}
// Shrinks droplet by 0.01
shrink := {body[1].target += -volstep; print body[1].target};
// combination of shrink and gogoconjugate, convenience command
shrinkgoconj := {shrink; gogoconjugate}
// combination of calc and gogoconjugate, convenience command
calcgoconj := {calc; gogoconjugate}

//Commands to deal with separation of the droplet from one of the cylinders.
// Use "separateupper" to cut off from the upper cylinder, and "separatelower" to cut off from the lower cylinder. 
deleteupper := {delete edges where on_constraint up_energy}
deletelower := {delete edges where on_constraint low_energy}
// unsets vertices on boundary
unsetupper := {unset vertices constraint up_energy; unset edges constraint up_energy}
unsetlower := {unset vertices constraint low_energy; unset edges constraint low_energy}
// USE THESE TO SEPARATE FROM UPPER OR LOWER
separateupper := {deleteupper; unsetupper}
separatelower := {deletelower; unsetlower}

// Automated shrinking- checks if boundary is small enough and then disconnects automatically, then evolves
auto_shrink :={
do {shrinkgoconj} while sum(edges where on_constraint low_energy,length) > 4*min_lent and sum(edges where on_constraint up_energy,length) > 4*min_lent;
if sum(edges where on_constraint low_energy,length) < 4*min_lent then {separatelower;separatelower;g 100; U; g 100;V 5;g 100; U; g 100};
if sum(edges where on_constraint up_energy,length) < 4*min_lent then {separateupper;separateupper;g 100; U; g 100;V 5;g 100; U; g 100};
print body[1].target;
}
//print angle and energy, change file location to your own
giv_ang := {printf " %f,",cyl_angle_degree >> "C:/Users/Horep/Documents/Extra Study PDF Files/Surface Evolver Project/Surface Evolver Files/garbagefiles/angles.txt"}
giv_energy := {printf "%10.20g,",total_energy >> "C:/Users/Horep/Documents/Extra Study PDF Files/Surface Evolver Project/Surface Evolver Files/garbagefiles/energies.txt"}
giv_area := {printf "%10.20g,",total_area >> "C:/Users/Horep/Documents/Extra Study PDF Files/Surface Evolver Project/Surface Evolver Files/garbagefiles/area.txt"}
// changes the angle by some amount, and then evolves, repeats until 90degrees
doangleiteration := {
while cyl_angle_degree < 90 do { gogoconjugate;
giv_ang;
giv_area;
cyl_angle_degree:=cyl_angle_degree+1}}
// outputs the width along the cylinder normal for upper and lower cylinder
check_width:=
{
cyl_angle:=cyl_angle_degree*pi/180;
width_upper := max(vertices vv where on_constraint up_energy,vv.x*cos(cyl_angle)+vv.y*sin(cyl_angle))-min(vertices vv where on_constraint up_energy,vv.x*cos(cyl_angle)+vv.y*sin(cyl_angle));
width_lower := max(vertices vv where on_constraint low_energy,vv.x)-min(vertices vv where on_constraint low_energy,vv.x);
print "lower width";
print width_lower;
print "upper width";
print width_upper
}
// Removes the display cylinders, leaving only the droplet. Useful for visualisation.
remove_cylinders :=
{
dissolve facets where fixed;
dissolve edges where valence ==0;
dissolve vertices where valence ==0;
print total_area
};

set background white; showq;