Let’s talk about boundary conditions

A simulation using finite element analysis involves solving a large set of equations. To get a meaningful solution we need to ensure that we are solving both the correct equations and correctly solving the equations. And an essential aspect of a correct set of equations is a correct set of boundary conditions, so let’s make sure we are setting the right boundary conditions!

This post may seem very basic at first, but keep on reading because we will dig into some subtleties of boundary conditions. Those of you who already know who Dirichlet and Neumann were might skip down to section 4 directly.

On the other hand, if you just finished your CAD model in Onshape and now want to run your first simulation, then I have two things to say:

1. Good for you!
2. This post is exactly what you want to read first.

1 Type of simulations

You want to “perform a simulation,” but what kind of simulation? Do you want to simulate…

• thermal conduction?
• solid mechanics?
• thermo-mechanics?
• fluid mechanics?
• conjugate heat transfer?
• electromagnetism?
• particle or radiation transport?
• chemical diffusion?
• quantum mechanics (i.e. the Schrödinger equation)?
• …

As you can imagine the list goes on. There are a lot of things one can simulate, but first things first: the type of simulation you need to perform depends on…

a. the quantity you’re interested in computing, such as
• stress
• temperature
• electrical field
• fission power
• etc.

b. what the quantity you're interested in depends on, such as does
• stress depend on temperature, loading rate or other quantities
• the electric field depend on temperature and vice versa
• the fission power depend on temperature
• etc.

Onscale Solve hosts several solvers specialized for different simulation types, all requiring boundary conditions. For thermomechanical problems it uses Onscale’s next-generation finite-element solver, Reflex.

2 Mechanical resistance

Let us assume for generality that we designed part of a super ultra-high-tech device that has to withstand the weight of a person standing on it (a step to reach the top of a tall tank?). So the natural question is how much weight can it support? Well, it is indeed a valid question, but not fully formulated. Why? Because we have to say how the part is

a. supported, and
b. loaded.

Figure 1. 1 \text{m} \times 0.5 \text{m} \times 0.1 \text{m} step designed in Onshape

2.1 Supports

When people stand on our step, their weight is transferred to the ground, and how the weight (or force) of the person is transferred affects the stresses that our device will experience. Let us assume the person weighs 1,000 Newtons. One option for supporting the step is to lay the step on the floor. This way, when somebody steps on it, the weight is directly transferred to the floor through the base. The compressive stress on the step is:

\begin{align*} \sigma = \frac{1000}{1 \cdot 0.5} = 0.002 \ \text{MPa} \end{align*}

The drawback is that our user only gained 10 cm in height.

Option two, if the layout allows it, is to fix the two lateral ends so the weight would be transferred to the ground by the two lateral surfaces. The maximum stress now rises to approximately 0.1 MPa. But we still need to be able to fix (to weld, perhaps) the lateral surfaces to two fixed locations in the tank. Chances are that we can only fix one of the sides and leave the other one free, configuring a cantilever. The maximum stress is now 0.6 MPa, but we have many more options on where to put the step. And of course there are far more variants about how to support the device, but you get the general idea.

Figure 2. Schematics of the three possible ways of supporting the step

Figure 3. Step leaning on the floor. Due to a non-zero Poisson’s ratio, the step expands horizontally as it is compressed vertically. The stress is uniformly equal to the force divided by the area.

Figure 4. Step fixed through both ends. Maximum displacement is at the center of the step and stresses distribute symmetrically with maximum values at the fixed ends.

Figure 5. Cantilevered step. Maximum displacement is at the free end, maximum stresses occur at the fixed end.

The discussion above should make crystal clear the following two points:

1. It is mandatory to define at least one fixation in static mechanical problems (i.e. where we solve for the stationary state of an object) unless you are solving something floating in space such as a satellite, but that’s out of the scope of this post.
2. The results depend significantly on what the fixations are.

2.2 Loads

The other requirement to fully define the problem is to specify the loads applied to the part. In the previous section, how were we able to give figures for stress without saying how the loads were distributed? Well, we assumed the weight was uniformly distributed across the top face. Indeed, that is what happens when we click on a surface in Solve and say “this face has a force.”

But this is not true, is it? People will be standing on their two feet, which will not span a rectangle of 1~\text{m} \times 0.5~\text{m}. Even more, at some point the person will be standing on only one foot when climbing over the step. So we might want to define a shoe-shaped surface and apply a load only to that face instead of doing it over the whole rectangle. Then, who says that the weight distributes uniformly on the shoe-shaped surface? And so on. You get the point.

3 Partial conclusions

Point number one above is the main explanation of why the apparently reasonable and valid question “how much weight can my part support?” is ill defined. It has to come with at least a description of the fixations. And, even more, a correctly-posed question ought to also include details about where the loads are prescribed and how they are distributed over the surfaces.

To finish the first part of the article, allow me to say that fixations and restraints are also called “essential” or “Dirichlet” boundary conditions. And loads are called “natural” or “Neumann” boundary conditions. Hence, this explains the title of the post and the links to those mathematicians entries in Wikipedia in the introduction.

4 Some advanced topics

Now that we have introduced boundary conditions, let us briefly dig into a couple of advanced features included in OnScale Solve.

4.1 Released degrees of freedom

Thoughtful readers might have wondered how the results from fig. 3 were obtained. The usual (and easy) way would be to fully restrain the base. However, if the Poisson’s ratio \nu is not zero—and it never is—a fully restrained base leads to a solution which is rarely found in practice (fig. 6). Our previous solution from fig. 3 shows an homogeneous stress distribution and allows the step to freely expand horizontally—which may not be fully physical since some friction is expected between the expanding face and the soil where it is standing, but it is still better than not allowing any expansion at all.

Figure 6. Step leaning on the floor with a fully fixed base.

Now, how did we manage to get the analytical stress distribution \sigma = F/A? The answer is in the essential boundary condition of the base. It has two distinctive features:

a. It has a cylindrical coordinate system attached.
b. It restrains only two out of the three cylindrical degrees of freedom: it restrains the angular coordinate \theta and the axial coordinate z but not the radial coordinate r.

These settings have to be defined using code in the SimAPI Python console (click on the <> icon in the right-bottom toolbar in OnScale Solve). First, set the cylindrical coordinate system as:

coord_cyl = on.CoordinateSystem(
   vector_1=[1, 0, 0],
   vector_2=[0, 1, 0],
   center=[0.5, 0.25, 0],
   system_type="cylindrical",
   alias= "cyl"
)

Then, set the boundary condition at the bottom face to release the first degree of freedom (called x for now even though it is the r coordinate) and fix the other two. Make sure you attach the cylindrical coordinate system to the BC:

   restraint = on.loads.Restraint(x= False, y=True, z=True, coord_system=coord_cyl, alias='Base' )
   restraint >> geometry.parts[0].faces[5]

Voilà! This trick can be used to create uniaxial compressive states even with a non-zero \nu. The location of the coordinate system center does not affect the stress distribution, but can affect the displacement. For example, if we choose one of the vertex locations as the coordinate system center, that vertex will remain fixed while the other three vertices will expand radially outward.

4.2 Space-dependent loads

What if we want to distribute the weight acting as the vertical load in a non-uniform way? We can use algebraic expressions to set the natural boundary condition at the top face. Let us assume that we want to “simulate” two fuzzy feet and that we come up with the following expression for the pressure in the upper face:

\begin{align} p(x,y) = 2000~\text{Pa} \cdot 2\pi \cdot \sin( 2 \pi x)^2 \cdot \sin(2 \pi y) \end{align}

Figure 7. Pressure distribution from an algebraic expression for two fuzzy feet.

This expression gives the pressure distribution depicted in fig. 7. Let me explain eq. 1 above:

• 2000~\text{Pa} is the “uniform” pressure from the previous cases (recall it was a force of 1,000 Newtons uniformly distributed over an area half a squared meter).
• 2\pi is a normalization factor needed to account for the fact that the distribution is not uniform any more (it is the inverse of the integral over the rectangle of the two factors that follow, feel free to do the math).
• \sin(2\pi x)^2 gives the two bumps in the x coordinate over a length of one meter (i.e. two fuzzy feet from left to right)
• \sin(2\pi y) gives the single bump in the y coordinate over a length of half a meter (i.e. front-back symmetry)

We can fire up the SimAPI console again and instead of a natural BC setting a total force, write:

   pressure = on.loads.Pressure(“2000 *2*pi * sin(2*pi*x)^2 *sin(2*pi*y)”,  alias= ‘Pressure 1')
   pressure >> geometry.parts[0].faces[4]

Fig. 8 shows how the results of the first case change in this case. The expansion was uniform both in the vertical and horizontal directions. Now the part is both compressed and expanded at different locations in both directions. Watch out!

(a) Top view

(b) Bottom view

Figure 8. Step leaning on the floor with a non-uniform pressure given by eq. 1

5 Final conclusions and homework

Besides the partial conclusions already discussed in sec. 3, let us add a couple of important points here:

• Care has to be taken about how both the restraints and the loads of a mechanical case (or the boundary conditions of a general partial differential equation) are defined.
• Sometimes fully fixing a face is not physically correct and there are more accurate alternatives available.
• Onscale Solve allows to define boundary conditions using algebraic expressions.

5.1 Homework

a. Solve the other two steps with the funny expression above.
b. All of the cases discussed so far involve symmetries in the geometry and in the loads. Exploit them to use less nodes to get the same accuracy or, conversely, get more accuracy with the same number of nodes.
c. Explore using other expressions in your BCs.