BGA Example bga-4 for the Surface Evolver

This file is the same as bga-3.fe, with some additions to calculate the vertical force on the upper pad.

Notable features:

• There are seven methods of force calculation illustrated here, which can all be invoked by the one command "forces". The first six methods, "zforce1" through "zforce6", are scripts defined in the file "zforce.cmd". The seventh is a named quantity, "zforce7", defined in the datafile.
• Most of the force calculation methods use the principle of virtual work, that is, the force associated with a parameter like height is the rate of change of energy with respect to the parameter. If global constraints like volume are present, then the volume change associated with the parameter change can be compensated for by subtracting from the energy the Lagrange multiplier of the constraint times the change in the constrained quantity value. This method requires that the surface be at equilibrium for the exact form of the deformation to be irrelevant.
• The first method is a one-way difference quotient moving only the top pad, by changing the value of "height". That is, zforce = -(energy(height+dz)-energy(height))/dz. A "recalc" is done after changing "height" to move the vertices on the top pad to the new height and recalculate the energy and volume. Note that the height is put back to its original value after the calculation is done; commands like these should leave the surface as close as possible to the state they found them in.
• The second method is a central difference moving just the top pad. Numerically, this is more accurate, but here the gain in accuracy is negligible compared to other inaccuracies present.
• The third method is a one-way difference quotient, moving the top by changing "height" and moving the other surface vertices by linear interpolation in z. The reason for doing this is that at equilibrium, the form of the deformation does not matter to first order, but it does to second order. Moving just the top pad puts all the deformation into the top row of facets, making a larger second-order error term.
• The fourth method is a central difference with linear interpolation in z.
• The fifth method is to actually re-evolve the surface after changing the height. This is the most fool-proof method, but also the most time-consuming. Also, if you want the surface to be exactly back to its original state, extra work must be done, for example saving all the original vertex coordinates in extra vertex attributes.
• The sixth method is geometric: calculate the pressure times area on the top pad and subtract off the vertical force exerted by the facets around the rim of the top pad (edge length times surface tension times horizontal component of facet normal). This turns out to be a bad method, because it turns out the facets that adjoin the rim at one vertex also matter. Also, it becomes even worse in the quadratic and higher order models.
• The seventh method is done by setting up a named quantity for the variational integrals of energy and volume for the deformation vector field, here (0,0,z/height). This gives the limit of the difference methods as dz goes to zero. This method should be recognizable to people familiar with the calculus of variations.
• The magnitude of the deformation, dz = 1e-5*height, was chosen to maximize accuracy in the central difference methods. As a rule, the maximum accuracy of single difference is half the digits of machine precision (dz being the other half of digits), and for central difference it is two-thirds the digits (one-third being used in dz).
• Some results of "forces" (after deleting some iteration output):
  After "gogo":

  zforce: 24.194068567832495  (single difference, top pad only)
  zforce: 24.197468861212691  (central difference, top pad only)
  zforce: 24.583191594966060  (single difference, linear move)
  zforce: 24.583596787882314  (central difference, linear move)
  zforce: 24.583612490876774  (central difference, linear move, reconverge)
  zforce: 17.323186420157086  (geometric)
  zforce: 24.583607839182264  (variational integrands)
  

  After another refinement:

  zforce: 18.549582101456959  (single difference, top pad only)
  zforce: 18.556654456010808  (central difference, top pad only)
  zforce: 18.667719719663484  (single difference, linear move)
  zforce: 18.668116454160590  (central difference, linear move)
  zforce: 18.668112948816425  (central difference, linear move, reconverge)
  zforce: 15.328939872764476  (geometric)
  zforce: 18.668111551491222  (variational integrands)
  

  After another refinement:

  zforce: 17.150752669910691  (single difference, top pad only)
  zforce: 17.165317478529115  (central difference, top pad only)
  zforce: 17.194671061796875  (single difference, linear move)
  zforce: 17.195067177529683  (central difference, linear move)
  zforce: 17.195068089392858  (central difference, linear move, reconverge)
  zforce: 15.557752604337814  (geometric)
  zforce: 17.195067593386000  (variational integrands)
  

• Some results running in quadratic mode (doing "M 2" before "gogo"):

  zforce: 16.655254887846617  (single difference, top pad only)
  zforce: 16.661739440001117  (central difference, top pad only)
  zforce: 16.729512453868740  (single difference, linear move)
  zforce: 16.729905937632829  (central difference, linear move)
  zforce: 16.729910171283297  (central difference, linear move, reconverge)
  zforce: 14.182445534020912  (geometric)
  zforce: 16.729909603674571  (variational integrands)
  

  After another refinement:

  zforce: 16.683272908342889  (single difference, top pad only)
  zforce: 16.696540795872274  (central difference, top pad only)
  zforce: 16.705052943895755  (single difference, linear move)
  zforce: 16.705446913197381  (central difference, linear move)
  zforce: 16.705446987212248  (central difference, linear move, reconverge)
  zforce: 14.746042425085822  (geometric)
  zforce: 16.705446803680104  (variational integrands)
  

• For me, the combination of simplicity and accuracy makes the central difference, linear move method the best. However, it is a good idea to implement more than one method (say, also the fifth method, reconvergence) at least during testing, and to test against known configurations. Be paranoid! Mistakes in defining the deformations are easy to make. Make big deformations and look at them, to make sure everything is moving the way you want. Don't forget to look at interior surfaces normally invisible.
• Now for a script to automatically generate force data for a range of volumes and heights. The command "runforce" is defined in the file runforce.cmd. (Click here to open the runforce.cmd command file in its own window.) The command is basically a nested loop, with volume the outer loop and height the inner loop. Because the changes in the surface due to changing the height in the inner loop are relatively small, the only iteration necessary is "hessian" three times. But because there are big changes each iteration of the outer loop, as the height is set back down, the command saves the vertex coordinates in a vertex extra attribute "old_coord". To guard against the collapse of the neck if the height gets too large, the script watches for the appearence of negative eigenvalues when "hessian" is run. However, it does not count the negative eigenvalues that appear at low heights, which are due to the bulging of the solder. The output file, "runforce.out" by default, is formatted to be directly input to a spreadsheet, say for conversion to a graph.
• Details of graph making in Excel:
  • File/Open runforce.out
  • Text import wizard: keep defaults all the way.
  • Use mouse to select all numeric data (including label row and column).
  • In Chart Wizard, select XY graph with data connected by lines.
  • Chart Wizard p. 2: series in rows.
  • Chart Wizard p. 3: fill in titles.
  • Chart Wizard p. 4: choose where you want chart.
  • Double click on x axis and use Scale tab to adjust axis limits.
  • Double click on y axis to adjust x axis crossing to -300.