Gull Wing Lead Example gullwing-3 for the Surface Evolver

This file replaces the solder-lead and solder-pad interface facets with edge integrals for surface energy, gravitational energy, and volume. This datafile is finally evolvable.

Design steps:

• First, we delete all the interface faces from the datafile, and the edges and vertices no longer needed. Can't forget to remove the deleted faces from the body, either. Alternatively, you may only want to remove faces one by one as they are replaced by integrals.
• In setting up edge integrals, my approach is through direct slicing arguments rather than a formal Green's Theorem or Stokes' Theorem approach. That is, I look at the geometry and say this bit of edge ought to contribute this much to the integral, and write down the proper formula.
• For the pad contact energy, I choose to slice in y, so each bit of edge contributes a rectangle from the y axis of area x*dy. From the orientation of the edges, I see that dy is positive when x is positive, and the orientation of all the edges is in the same direction all the way around, so the same integral can be used on all. Multiplying by the surface tension factor gives the energy integrand -gamlv*cos(thetap*pi/180)*x*dy. Note that since thetap was defined in degrees, it must be converted to radians here. Adding this to constraint 5 as the energy integrand gives

     constraint 5  // pad surface
     formula: z = 0
     energy:
     e1: 0
     e2: -gamlv*cos(thetap*pi/180)*x
     e3: 0
  

• For the contact energy on the toe of the lead, I slice in x. The height of a slice is (z - gap), and dx is positive on this edge (edge 37), so the area integrand is (z - gap)*dx. Adding the toe surface tension factor gives this energy integrand on constraint 10:

     constraint 10  // toe surface of lead, i.e. -y side
     formula: y = leadtoe+leadheel
     energy:
     e1: -gamlv*cos(theta_toe*pi/180)*(z-gap)
     e2: 0
     e3: 0
  

• For the contact energy on the +x side of the lead, I slice in y. The height of a slice is (z - ZLOWER), where ZLOWER is the height of the bottom curve of the lead as a function of y. This is a complicated piecewise defined curve, so I set up some macros for convenience:

     #define ZLOWER1 gap
     #define ZLOWER2 (orad + gap - sqrt(orad^2-(y-leadheel)^2))
     #define ZLOWER3 (gap+(rad+orad)*(1-cos(bend))-rad + \
               sqrt(rad^2 - (y-leadheel-(rad+orad)*sin(bend))^2) )
     #define ZLOWER ( (y < leadheel) ? ZLOWER1 : \
             (y < leadheel+orad*sin(bend) ? ZLOWER2 : ZLOWER3))
  

  Note the parentheses around macros that expand to compound expressions, so there are no operator precedence problems after expansion. Also note the line continuation backslashes in multi-line macros; be sure these are the very last characters on the lines. Using ZLOWER as a function of y means that things will not be well-behaved if the bend angle were 90 degrees and the solder reach the vertical part of the bend, but hopefully that combination of circumstances will not arise.

  The relevant edge (edge 38) has positive dy, so the area integrand is (z - ZLOWER)*dy. Putting in the surface tension gives the energy integrand on constraint 11:

     constraint 11  // +x side of lead
     formula: x = leadshift + width/2
     energy:
     e1: 0
     e2: -gamlv*cos(theta_side*pi/180)*(z-ZLOWER)
     e3: 0
  

  If you consult the datafile, you will see that this is not the final integrand; gravitational energy will be added below.
• The contact energy on the -x side of the lead follows the same reasoning, except that dy is negative, so there is a change in sign:

     constraint 12  // -x side of lead
     formula: x = leadshift - width/2
     energy:
     e1: 0
     e2: gamlv*cos(theta_side*pi/180)*(z-ZLOWER)
     e3: 0
  

  If you are wondering why the sign changed here and not on constraint 5 above, recall that constraint 5 was a single application of Green's theorem, and here we have two separate surfaces. Also, by checking signs along the various edges, you see that the chosen signs give the proper results.
• For the contact area on the bottom of the lead, I choose to slice in x. Thus each slice of area is a thin strip running along the bottom from the toe up to a bit of edge 39, which is the contact line crossing under the heel of the lead. Since dx is negative on this edge, the area integrand is -striplength*dx, which ultimately gives the integrand for constraint 13:

  energy:
     e1: gamlv*cos(theta_lower*pi/180)*( y < leadheel ? y-leadtoe :
          ( y < leadheel+(orad+rad)*sin(bend) ?
          -leadtoe + orad*asin((y-leadheel)/orad) :
          -leadtoe + orad*bend +
                   rad*(bend-asin((leadheel+(orad+rad)*sin(bend)-y)/rad)) ))
     e2: 0
     e3: 0
  

• Now on to the volume corrections. Since the Evolver calculates body volume with the integral of z dx dy over facets, vertical facets can be omitted. This leaves just the bottom of the lead. We could slice in x, but this leads to some pretty messy formulas, especially in the gravitational energy below. It turns out easier to slice in y. For a bit dx dy of the bottom at height z, the bit of volume is z*dx*dy. We don't have edges running along the bottom of the lead (well, display edges, but we're not using those for energy), so we will integrate along the side edges we do have and use ZLOWER instead of z. And if we integrate the dx out from the y axis, we get a factor of x, so the integrand becomes ZLOWER*x*dy. Checking the orientations of the edges, we see that x*dy always comes out positive, so no sign changes are necessary on the different edges. However, a further sign check is necessary for content integrals. The rule is that an edge content integral adds to a body that has a positive facet with the given edge in positive orientation around the facet. Our edges are negatively oriented with respect to their facets, which are positive on the body, so there is an additional negative sign. Thus the content integrand comes out

     content:
     c1: 0
     c2: -ZLOWER*x
     c3: 0
  

  We make this integrand the content integrand of constraints 11, 12, and 13. Can't forget constraint 13, even though its edge initially has dy = 0; later on it will not. On the other hand, constraint 10 edges always have dy = 0, so we don't need to add to that constraint.
• The gravitational energy works much like volume. For a bit dx dy of the bottom at height z, the bit of gravitational energy is G*den*z^2/2*dx*dy. Again, we slice in y and integrate out dx to get G*den*ZLOWER^2/2*x*dy. Checking the orientations of the edges, we see that x*dy always comes out positive, so no sign changes are necessary on the different edges. So we add this integrand to the energy integrands of constraints 11, 12, and 13.
• One way to check the values of the individual constraint energies is to do convert_to_quantities; show_all-quantities, use the 'v' command to show the constraint energy values, and compare with the energies of the facets on the respective constraints in gullwing-2.fe. For example, for the pad contact in gullwing-2.fe,

     Enter command: print sum(facet where on_constraint 5,area*tension)
        -301.923239058512
  

  This should be equal to constraint_5_energy in gullwing-3.fe. If they disagree, find the error in one or the other.
• Actually, the two numbers won't always agree when curved edges are involved. For those, refine gullwing-2.fe several times to get better approximation to curves, both for shape and for numerical integration accuracy.
• G is used in constraint energy formulas instead of mr_grav, so all gravity can be changed with one 'G' command.
• When the energy is principally surface tension, it turns out the scale factor is about 0.2 when the tension is 1, and varies inversely with tension (simple dimensional analysis). Since the tension here is 0.45, we'll raise the upper bound on the scale factor from the default of 1 to 2:

      scale_limit 2
  

Evolution:

• Enter command: l 40
  This breaks up the longest edges.

• Enter command: r
  Refine to give some flexibility to the surface.

• Enter command: g 10
  Just to get things started, do 10 gradient descent steps. Things look good, and the scale factors are nice and large, near 0.2.

• Enter command: U
  Things are going well so far, so we switch to conjugate gradient mode for faster evolution.

• Enter command: g 20
  The contact line moves up to the tops of the sides, but way down at the toe corners. The scale factor gets a little small at the end, and we notice in the display that there's a skinny triangle we will have to deal with now. Two things we could do: vertex averaging to spread the vertices apart, or delete the skinny triangles. We'll do vertex averaging.

• Enter command: u;V;u;V
  This is a good sequence to use whenever the triangulation gets a little ugly. It usually gets a nice internal triangulation, but the edges need a little more help.

• Enter command: refine edge where valence==1 and not fixed
  Introduces more vertices around the edges for a better triangulation.

• Enter command: g 20; u; V; g 20;
  More iteration, with u; V; just on general principles.

• Enter command: r; V; u; V; u; g 20; V; u; V; u; g 40
  Refining, and more evolution. Note the scale factor is often above 1, which is why we set scale_limit to 2 earlier. At this stage, you can also see that the contact line on the pad is starting to round off in the corners. The solder wets all the way into the corners in the smooth limit only if the contact angle is zero.

• Undercutting
  For some sequences of iterations, a problem can appear on the toe corners. The surface wants to suck under a bit. This is a very small effect here, so no action needed, just keep an eye on it. In the smooth limit, the contact lines will curve to be tangent to the vertical corner edge, and may go all the way down to the bottom of the lead. One cure would be to fix the corner point down at the bottom of the lead. A similar suck under happens on the heel.

  
  

• Enough for this file. For final polishing touches, go on to gullwing-4.

Forces.

Here's the basic x force calculation command:

   dx := 0.0001
   calc_xforce := {
      delta_x := dx;
      move_x;
      hi := total_energy - body[1].pressure*(body[1].volume-body[1].target);
      delta_x := -2*dx;
      move_x;
      lo := total_energy - body[1].pressure*(body[1].volume-body[1].target);
      delta_x := dx;
      move_x;
      xforce := -(hi - lo)/2/dx;
   }

   define vertex attribute alpha real

   calc_alpha := {
      minz := min(vertex where on_constraint 13 or on_constraint 23, z);
      if ( minz <= 0.0 ) then
         printf "WARNING: contact line hits pad.\n";
      foreach vertex vv do
      { if ( vv.on_constraint 1 or vv.on_constraint 5 ) then
          { vv.alpha := 0; continue; };
        if ( vv.fixed ) then { vv.alpha := 1; continue; };
        vv.alpha := vv.z >= minz ?  1 : vv.z/minz ;
      }
   }

   move_x := { calc_alpha;
               leadshift := leadshift + delta_x;  // do parameters first
               set vertex x x+alpha*delta_x;
             }