• From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Topic: How to implement surface tension

Hello, I have some practical questions

Just some short background;
I'm need to simulate processes involving surface tension. In it's simplest form it's just a surface traction;
http://www.forkosh.dreamhost.com/mathtex.cgi?{\bf%20t}_s=2H\gamma_s{\bf%20n}
where H is the surface curvature (surface divergence of the normal). γ is a material parameter.

In almost all cases, it's best to rewrite the integral with the divergence theorem, and you end with a load dependent only on the normal, and an edge load where the free  surface ends (this makes it possible to use linear elements, which is very nice)
What you end up with in the weak form is
http://www.forkosh.dreamhost.com/mathtex.cgi?\int_\Gamma{\bf%20w}\cdot{\bf%20t}_s{\rm%20d}a=\int_C\gamma_s{\bf%20w}\cdot{\bf%20m}{\rm%20d}s-\int_\Gamma\gamma_s\nabla_s\cdot{\bf%20w}{\rm%20d}a
In 2D (plane strain or axisymm), Γ is the free surface, C would be the endpoints (if existent) of Γ. m is the tangent of Γ. ∇s would be the derivative along the edge.

Now let's put the theory aside and look at the implementation! I'm seeing three possible options here. I'll list the pros and cons as i see them, but i would appreciate some input on what way is preferable.

  1. Hard code a the edge loads in the input file. This method is of course limited not to have large deformations, which is typical when using surface tension simulations. Not really a suitable solution at all.

  2. Introduce a new type of boundary condition to bcType, i.e. IsotropicSurfaceTensionBC. One problem is the edge term from the divergence theorem, which would need to be applied where the free surface ends. Possible a point

  3. Implement a new boundary element. The element will just as I've just finished such an element, but there is some minor issues.
    This edge element should, just as ConstantEdgeLoad, fit both SM and FM problems just fine, but I either let the Domain default to V_u, V_v, P_f (adding zeros to my load vector) or specify it manually to just V_u, V_v using giveDofManIDMask.

The V_u, D_u independence points towards using a load type, while things like the boundary term seems to fit an element more nicely.

Now, this load would also change if you have axisymmetry. If implemented as an element, I guess it would simply be as with any other element, implement another element for each domain type. Not sure how constantboundaryload handles this sort of thing.

In short; I can't quite decide on what approach i should go for, 2 or 3.

Small note;
One could also generalize this even more, the surface tension doesn't have to be isotropic. In that case, I can't see any elegant solution but to use boundary elements, with respective surface material models. However, this is very uncommon it's not worth considering

  • Registered: 28-07-2008
  • Posts: 71

Re: How to implement surface tension

Hello Mikael,

From my experiences on the problem I am working on, I'd lean towards adding a new boundaryLoad. If I understand correctly, it seems like there will have to be some code at the element level to deal with this load. Depending on how the free surface is implemented, I can imagine that the element would also be able to determine if the free surface endpoint occurs on the element. If not, this approach may get a bit messy for the endpoints.

You will also need to include a way to define which "surface" the load acts on in your element input (similar to existing loads).

I have found that the basic OOFEM framework is pretty flexible. There is generally a function or class of constants that provide or suggest a path to deal with most of the issues I have had implementing the problem I am trying to solve.

If you define something like "Element_SurfaceTensionLoadSupport", you could use "testElementExtension(Element_SurfaceTensionLoadSupport)" to see if a particular element can handle this loading.

I have also found the element functions giveBCGeoType and giveBCValType helpful to figuring out how to handle a load.

Regards,
Erik

3 Reply by Mikael Öhman (edited by Mikael Öhman 16-04-2010 16:58:48)

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: How to implement surface tension

Hi Erik.

I really appreciate your input on this, however, I was working some more on this problem, and currently I'm leaning towards implementing a new type of boundary element.
Surface tension can often be viewed as a sort of skin which is trying to contract, so adding a layer of elements on the surface does seems quite logical to me. I could also create some simple surface energy material models (the simplest being just isotropic, one constant).

This would let me add surface tension to any existing elements on OOFEM, which to me sounds like a big win.

Even if this means i have to code specifically for FM or SM module (perhaps it could be done generally) i still think it's the cleanest solution.

Edit: Fixed a typo, it *does seem quite logical*

4 Reply by bp

  • bp
  • Administrator
  • Offline
  • From: Prague, Czech Rep.
  • Registered: 02-03-2007
  • Posts: 1,047

Re: How to implement surface tension

Hi Mikael,

in my opinion, it would be more natural to implement it as a sort of boundary condition. It is already a some time ago, when I started to think abut some changes. One of them is related to similar issues you have. I would like to make boundary conditions responsible for evaluating resulting terms, instead of elements.  Some steps have already been done, perhaps the most important one is the introduction of interpolation classes at element level. Element provides access to its geometry and interpolation, so that the evaluator (could be a boundary condition, for example) could compute its contribution. It can even dynamically set up its own integration rules (based on element interpolation order and boundary condition nature) and evaluate the response. The element interpolation class probably provides all necessary methods one may need (evaluating shape functions, their derivatives, mapping from local to global (and vice versa) element coordinates, etc.).
Of course, this requires some changes, as up to now boundary conditions are not assumed to evaluate characteristic terms, but this could be easily added. In my point of view, this approach is more clean than providing the element support for evaluating boundary conditions, at least implementation for boundary condition is in one place, not distributed around in many elements.
I have tested this approach already when implementing iso-geometric elements (with bspline and nurbs interpolation), where the structural capability was implemented in a similar fashion. 

Borek

bp's Website

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: How to implement surface tension

Hi Borek, I have some comments

While I can agree that isotropic surface tension can be viewed naturally as a boundary condition, if one considers anisotropic surface potential (not completely unreasonable, for example an anisotropic skin around an object during heating/cooling), one would have a general case with material models.
I don't see how this could be implemented other than introducing more boundary conditions. It is also somewhat strange for a boundary condition to have a material model.

There is also a slight problem with the edge integral. The boundary condition would also need to know where the surface ends. This comes naturally for elements, i.e. some type of boundary flag, but is very odd for a boundary condition.

I think it's logical to view this as a skin with it's own finite elements, as well as a boundary condition, so I guess it's down to what implementation is most practical.

Since I've already implemented this as elements (by far the easiest thing to do for a beginner such as me), I'll continue with that one and we can discuss changes later (the code is dead simple, changes are easy).

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: How to implement surface tension

Upon further investigation, a new type of boundary condition does sound very appealing.

It would have been very nice to handle 2D, 2D-axi, 3D, combined with either displacement or velocity, for any kind of element type (tri, quad, tet, hex). For isoparametric mappings the code is very similar.
If someone needs nonisotropic surface tension, they could always implement that normally with elements. I've decided it's probably not common enough to decide this.

The only problem i have is how one would handle the edge of the free surface. I'm not sure how I would even specify it without adding more information to existing elements.

I could split it up the integrals, and having to manually add an free-surface edge as a boundary condition to effected elements, but i see that as a last resort.