File: ringblob.htm

package info (click to toggle)
evolver 2.70+ds-4
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 17,148 kB
  • sloc: ansic: 127,395; makefile: 209; sh: 98
file content (284 lines) | stat: -rw-r--r-- 11,219 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
<!DOCTYPE HTML>
<HEAD><title>Surface Evolver Documentation - Torus Example</title>
<link rel="stylesheet" type="text/css" href="evdoc-style.css" />
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" >
</head>
<body>

<h1 class="center">
<a href="http://www.susqu.edu/brakke/evolver/evolver.htm" class="comic">
Surface Evolver</a> Documentation</h1>

<a href="evolver.htm#doc-top">Back to top of Surface Evolver documentation.</a>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<a href="index.htm">Index.</a><br>
<a href="column.htm">Next: Liquid column example.</a> <br>
<a href="tutorial.htm#tutorial">Back to top of tutorial.</a>

<a   id="spinning-ring-example"></a>
<a   id="ringblob"></a> <h1>Example:  Ring around rotating rod</h1>


This example consists of a ring of liquid forming a torus around a
rod rotating along its long axis (z axis) in weightlessness.  The liquid has
controllable contact angle with the rod.  The interesting question is
the stability of the ring as the spin increases.
<p>
The effect of the rotation is incorporated in the energy through an integral
using the Divergence Theorem:
<pre>
        ///                    
  E = - ||| (1/2) p r^2 w^2 dV 
	///B                   

        //
    = - ||       (1/2) p w^2 (x^2+y^2) z <b>k&middot;dA</b>
        //bdry B

</pre>
where B is the region of the liquid, r is radius from the axis,
p is the fluid density and w is the angular velocity.
Note the energy is negative, because spin makes the liquid want to move
outward.  This has to be countered by surface tension forces holding
the liquid on the rod.  If p is negative, then one has a toroidal
bubble in a rotating liquid, and high spin stabilizes the torus.
The spin energy is put in the datafile using the named quantity syntax
(see below).
"<code>centrip</code>" is a user-chosen name for the quantity, "<code>energy</code>"
declares that this quantity is part of the total energy, 
"<code>global_method</code>"
says that the following method is to be applied to the whole surface,
"<code>facet_vector_integral</code>"  is the pre-defined name of the method
that integrates vector fields over facets, and "<code>vector_integrand</code>"
introduces the components of the vectorfield.
<p>
The rod surface is defined to be constraint 1 with equation x^2 + y^2 = R^2,
where R is the radius of the rod.  The contact energy of the liquid with
the rod is taken care of with an edge integral over the edges where the
liquid surface meets the rod:
<pre>
      //                          /                  /
  E = || -T cos(a) dA = -T cos(a) |  z ds = T cos(a) | (z/R)(y<b>i</b> - x<b>j</b>)<b>&middot;ds</b>
      //S                         /bdry S            / bdry S
</pre>
Here S is the rod surface not included as facets in bdry B,
T is the surface tension of the free surface, and a is the
internal contact angle.  A more intuitive way to arrive at this integral
is to think in cylindrical coordinates and
 imagine the replaced surface of the rod between the contact line and
z = 0 to be divided into thin vertical strips of height z and width R dtheta.
The energy of a strip is <code>-T cos(a) z R dtheta</code>.  Converting back
to Cartesian coordinates, dtheta = (x dy - y dx)/(x^2 + y^2), so
<pre> 
                 //                                   //
  E =  -T cos(a) || z R (x dy - y dx)/R^2 = -T cos(a) || (z/R)(x dy - y dx)
                 //                                   //
</pre>
<p>
Constraint 2 is a horizontal symmetry
plane.  By assuming symmetry, we only have to do half the work.
<p>
Constraint 3 is a one-sided constraint that keeps the liquid outside the
rod.  Merely having boundary edges on the rod with constraint 1 is not
enough in case the contact angle is near 180 degrees and the liquid volume
is large.  Constraint 3 may be put on any vertices, edges, or faces
likely to try to invade the rod.  However, it should be noted that
if you put constraint 3 on only some vertices and edges, equiangulation
will be prevented between facets having different constraints.
<p>
Constraint 4 is a device to keep the vertices on the rod surface evenly
spaced.  Edges on curved constraints often tend to become very uneven,
since long edges short-cutting the curve can save energy.  Hence the
need for a way to keep the vertices evenly spread circumferentially, but
free to move vertically.  One way to do that is with another constraint
with level sets being vertical planes through the z axis 
at evenly spaced angles.  Constraint 4 uses the real modulus
function with arctangent to create a periodic constraint.  Each refinement,
the parameters need to be halved to cut the period in half.  This is done
by redefining the "<code>r</code>" refinement command at the end of the datafile.
Note that autorecalc is temporarily turned off to prevent projecting
vertices to the constraint when it is in an invalid state.  Also note the
pi/6 offset to avoid the discontinuity in the modulus function.
pi/6 was cleverly chosen so that all refinements would also avoid
the discontinuity.
<p>
One way of detecting stability is to perturb the torus and seeing if
it evolves back to equilibrium. The datafile defines a command
<pre>
    perturb := set vertex y y+.01 where not on_constraint 1
</pre>
This sets the y coordinate of each vertex to y+.01.
This command is defined in the "<code>read</code>" section
at the end of the datafile, where you can put whatever commands you
want to execute immediately after the datafile is loaded.
To detect small perturbations, and get numerical values for the size
of perturbations, the y moment of the liquid is calculated in the
named quantity "<code>ymoment</code>".  It is not part of the energy, as
indicated by the <code>info_only</code> keyword.  You can see the value
with the `<a href="single.htm#v"><code>v</code></a>' command.
<p>
A better way to check stability is to examine the eigenvalues of
the Hessian matrix.  First, evolve normally to reasonably near
an equilibrium.  Then use the <a href="commands.htm#hessian-command">hessian</a>
command several times to converge exactly to the equilibrium. Then
use the command "<code><a href="commands.htm#ritz">ritz</a>(0,5)</code>"
to display the 5 eigenvalues of
the Hessian nearest 0.  By gradually raising the spin and tracking
the lowest eigenvalue, one can detect the onset of instability, where
the lowest eigenvalue becomes 0.  Note that the datafile toggles on
<a href="toggle.htm#hessian_normal">hessian_normal</a> so that the
Hessian only considers perturbations normal to the surface.

<table>
<tr><td>
<img src="ringblobbare.gif"  alt="ringblob"></td>
<td>  The initial ringblob skeleton, with
vertices and edges numbered.  Taking advantage 
of symmetry, just the top half is represented.</td></tr>
</table>


<pre>
// ringblob.fe

// Toroidal liquid ring on a rotating rod in weightlessness.
// Half of full torus
// Using second periodic constraint surface intersecting rod to
// confine vertices on rod to vertical motion. 

// Important note to user: Use only the 'rr' command defined at
// the end of this file to do refinement.  This is due to the
// nature of constraint 4 below.

// This permits drawing both halves of the ring
view_transforms 1
1 0 0 0
0 1 0 0
0 0 -1 0 
0 0 0 1

// Basic parameters.  These may be adjusted at runtime with the
// 'A' command.  Only spin is being adjusted in these experiments.
parameter rodr = 1  // rod radius
parameter spin = 0.0 // angular velocity
parameter angle = 30 // internal contact angle with rod
parameter tens = 1   // surface tension of free surface
#define rode (-tens*cos(angle*pi/180))  // liquid-rod contact energy
parameter dens = 1  // density of liquid, negative for bubble

// spin centripetal energy
quantity centrip energy global_method facet_vector_integral
vector_integrand:
q1: 0
q2: 0
q3: -0.5*dens*spin*spin*(x^2+y^2)*z

// y moment, for detecting instability
quantity ymoment info_only global_method facet_vector_integral
vector_integrand:
q1: 0
q2: 0
q3: y*z

// Constraint for vertices and edges confined to rod surface,
// with integral for blob area on rod
constraint 1
formula: x^2 + y^2 = rodr^2
energy: 
e1: -rode*z*y/rodr
e2: rode*z*x/rodr
e3: 0

// Horizontal symmetry plane
constraint 2 
formula: z = 0

// Rod surface as one-sided constraint, to keep stuff from caving in
// Can be added to vertices, edges, facets that try to cave in
constraint 3 nonnegative 
formula: x^2 + y^2 = rodr^2

// Constraint to force vertices on rod to move only vertically.
// Expressed in periodic form, so one constraint fits arbitrarily
// many vertices. Note offset to pi/6 to avoid difficulties with
// modulus discontinuity at 0.
parameter pp = pi/2    /* to be halved each refinement */
parameter qq = pi/6    /* to be halved each refinement */
constraint 4 
formula:  (atan2(y,x)+pi/6) % pp = qq

//initial dimensions
#define ht 2
#define wd 3

vertices
1  0   -wd 0  constraints 2    // equatorial vertices
2  wd    0 0  constraints 2
3  0    wd 0  constraints 2
4  -wd   0 0  constraint 2
5  0   -wd ht                 // upper outer corners
6  wd    0 ht
7  0    wd ht
8  -wd   0 ht 
9  0 -rodr ht constraints 1,4   // vertices on rod
10 rodr  0 ht constraints 1,4
11 0  rodr ht constraints 1,4
12 -rodr  0 ht constraints 1,4

edges
1   1 2 constraint 2  // equatorial edges
2   2 3 constraint 2
3   3 4 constraint 2
4   4 1 constraint 2
5   5 6               // upper outer edges 
6   6 7
7   7 8
8   8 5
9   9  10 constraint 1,4   // edges on rod
10  10 11 constraint 1,4
11  11 12 constraint 1,4
12  12  9 constraint 1,4
13   1  5        // vertical outer edges
14   2  6
15   3  7
16   4  8
17   5  9        // cutting up top face
18   6 10 
19   7 11 
20   8 12 

faces  /* given by oriented edge loop */
1   1 14 -5 -13 tension tens // side faces
2   2 15 -6 -14 tension tens  // Remember you can't change facet tension
3   3 16 -7 -15 tension tens  // dynamically just by changing tens; you have
4   4 13 -8 -16 tension tens  // to do "tens := 2; set facet tension tens"
5   5 18 -9 -17 tension tens  // top faces
6   6 19 -10 -18 tension tens
7   7 20 -11 -19 tension tens
8   8 17 -12 -20 tension tens 

bodies  /* one body, defined by its oriented faces */
1   1 2 3 4 5 6 7 8  volume 25.28

read  // some initializations
transforms off     // just show fundamental region to start with

// special refinement command redefinition
r :::= { autorecalc off;  pp := pp/2; qq := qq % pp; 'r'; autorecalc on; }

// a slight perturbation, to check stability
perturb := set vertex y y+.01 where not on_constraint 1

hessian_normal // to make Hessian well-behaved
linear_metric  // to normalize eigenvalues 

</pre>
<hr>
<a href="column.htm">Next: Liquid column example.</a> <br>
<a href="tutorial.htm#tutorial">Back to top of tutorial.</a>
<br>
<a href="evolver.htm#doc-top">Back to top of Evolver documentation.</a>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<a href="index.htm">Index.</a>
</body>
</html>