# Cork stopper - ‘Hello World!’¶

## General formulation of the problem¶

```
[2]:
```

```
%%tikz -l patterns -s 600,320
\input{fig/formulation.tikz};
```

### Governing equations¶

The problem for displacement \(\mathbf{u}:\Omega \to \mathbb{R}^3\) of isotropic solid at small strains reads

where

is the Cauchy stress tensor (with \(\varepsilon(\mathbf{u})\) being the symmetric gradient of \(\mathbf{u}\), \(\mathbb{I}\) the identity tensor, and \(\lambda\) with \(\mu\) the Lamé parameters), \(\mathbf{b}\) is the volumetric force density, \(\mathbf{g}\) is the normal traction on the Neumann part of the boundary of \(\Omega\), and \(\mathbf{u}_D\) is the prescribed displacement at the Dirichlet part of \(\partial \Omega\). When reducing the problem for an axial dilatation and a radial compression we suppose

and that \(\mathbf{g}\) is zero on the walls and points in the axial direction at both ends. The system then partially separates; the equation for the component \(u_x\) becomes independent. Its solutions than enters the equations for \(u_y\) and \(u_z\) via the traction on the boundary (here the Poisson’s ratio comes into play); however, we are interested in finding solely \(u_x\). When we integrate its equation over the cross-section with the area \(A\), we obtain

Since we want to describe the material properties rather with Young modulus \(E\) and Poisson’s ratio \(\nu\), we substitute

When supposing \(\nu = 0\), we end up with the reduced problem for \(u:=u_x\)

where \(b(x) = Ab_x(x)\) and \(g = Ag_x\). Weak formulation is

where \(V = \{ v \in H^1(0,L); v(0) = 0\}\)1. Note that \(H^1(0,L) \subset AC[0,L]\)2, thus \(\varphi\) is defined pointwise. Sometimes it is convenient to rewrite the previous variational formulation in the form

where \(a(u,\delta u) = \int_0^L A E u'(x) \delta u'(x)\, \text{d} x\) denotes a bilinear form given by the left-hand side of the formulation and \(L(\delta u) = \int_0^L b(x)\delta u(x)\, \text{d} x + g \delta u(L)\) is a linear functional that contains the external loading of the problem.

#### Geometric derivation*¶

Displacement of the bar cross-section located at \(x\) is given by \(\mathbf{u}(x)\). If we cut an element \([x,x+\Delta x]\) of the bar then the left cross section is displaced by \(u(x)=u_x(x)\) in the \(x\) direction, and the right one is displaced by \(u(x+\Delta x)\), therefore the strain, or the relative deformation, of the element is

Thus, for an infinitesimal element, the deformation at \(x\) is given by

Using 1D stress-strain constitutive relation, i.e. Hooke’s law, we get

where \(\sigma\) now represents only the \(xx\) compononet of the stress tensor \(\boldsymbol{\sigma}\).

```
[3]:
```

```
%%tikz -l arrows -s 400,200
\input{fig/equilibrium.tikz};
```

From static equilibrium of an infinitesimal element, see the previous figure, we get the following force balance

which can be recast into

Using the stress relation we get

If we suppose that \(E(x)=\mathrm{const.}\) and \(A(x)=\mathrm{const.}\), we get

which is the same relation we got in the previous derivation.

## Implementation¶

First of all, we must import `fenics`

library to `python`

interface. Moreover, we import `pyplot`

library that is used for graph rendering. And we will also need some features of the `numpy`

package.

```
[4]:
```

```
import fenics as fe
import matplotlib.pyplot as plt
import numpy as np
```

### Problem parameters¶

Let us set the material parameters of the problem. Following is the table that contains Young modulus and Poisson’s ratio for some common materials.

Material | Young modulus \(E\ [\mathrm{GPa}]\) | Poisson’s ratio \(\nu\ [-]\) |
---|---|---|

Rubber | \(0.01 - 0.1\) | \(0.4999\) |

Cork |
\(0.02 - 0.03\) | \(0\) |

Concrete | \(30\) | \(0.1 - 0.2\) |

Steel | \(200\) | \(0.27 - 0.30\) |

We decided to model a cork stopper, since cork material parameters justify the assumption \(\nu=0\) made in the governing equations. Thus, we choose appropriate dimensions of the cork stopper (radius 1 cm, length 5 cm) and some suitable magnitude of loading \(b\) and \(g\).

```
[5]:
```

```
# --------------------
# Parameters
# --------------------
E = 0.025e9 # Young's modulus [Pa]
A = np.pi*1e-4 # Cross-section area of bar [m^-2]
L = 0.05 # Length of bar [m]
n = 20 # Number of elements
b = 0.03 # Load intensity [Pa*m^-3]
g = 0.0005 # External force [Pa*m^-2]
```

### Geometry of the problem¶

Next, we define the geometry of the problem. There is always a possibility to use an external mesh defined in *.xml* file. In our case we use the predefined mesh provided directly by FEniCS. Initialization of the FEniCS object `IntervalMesh(n, x_1, x_2)`

creates an equidistantly discretized one-dimensional mesh with `n`

elements between the end points `x_1`

and `x_2`

.

You can invoke a basic description of any function/object by`help`

command. Try that for the object`IntervalMesh`

.

```
[6]:
```

```
%%capture
# %%capture supresses the output of this notebook cell, which is quite lenghty for the help command
help(fe.IntervalMesh)
```

```
[7]:
```

```
# --------------------
# Geometry
# --------------------
mesh = fe.IntervalMesh(n, 0.0, L)
```

We can visually check the geometry of the mesh using the FEniCS `plot()`

command. The plot object is then displayed by adding the pyplot command `show()`

.

```
[8]:
```

```
fe.plot(mesh)
plt.show() # this forces the plot to show in a pop-up window when in command line/script environment
```

### Function spaces definition¶

`FunctionSpace(mesh,type,deg)`

with three arguments.`mesh`

variable. Second and third are the type `type`

of the base functions and the order `deg`

of these functions respectively. We will present various types of function elements later in the course. In agreement with the definition of the function spaces given in the first section, we define the conforming space of finite elements, that is the *Lagrange/Continous Galerkin*element (hence the acronym

`CG`

) of order one, in other word, piecewise linear functions.Then we initialize the test and trial functions specified by the function space as `TrialFunction`

and `TestFunction`

objects, respectively.

```
[9]:
```

```
# --------------------
# Function spaces
# --------------------
V = fe.FunctionSpace(mesh, "CG", 1)
u_tr = fe.TrialFunction(V)
u_test = fe.TestFunction(V)
```

### Boundary conditions¶

As a next step, we specify the geometric boundary of the domain through which we apply the Dirichlet boundary condition. It is necessary to mark those entities of the mesh that constitute the boundaries \(\Gamma_{\text{D}}\) and \(\Gamma_{\text{N}}\).

```
[10]:
```

```
# ------------------------
# Boundary markings
# ------------------------
def left(x, on_boundary):
return fe.near(x[0], 0.0) and on_boundary
# analogously, we could define the marking of the right boundary
def right(x, on_boundary):
return fe.near(x[0], L) and on_boundary
```

In the previous, `on_boundary`

flag is `True`

for all points that lie on the boundary of the mesh and `near(x,x0)`

checks whether `x`

is close to `x0`

(by default, within the tolerance given by the global constant `DOLFIN_EPS`

), which respectively gives the left and the right end of the beam.

When you compare two real numbers in`python`

, always use certaintolerance. Value of the tolerance depends on the size of the numbers you compare. For unit sized numbers,`DOLFIN_EPS = 3e-16`

is sufficient.

Try it for yourselves. Print the outcomes e.g., of `0.1+0.1+0.1-0.3`

and `100.1+100.1+100.1-300.3`

:

```
[11]:
```

```
print(0.1+0.1+0.1-0.3)
print(100.1+100.1+100.1-300.3)
```

```
5.551115123125783e-17
-5.684341886080802e-14
```

There is an equivalent way of marking the boundary. It employs the `MeshFunction`

object that represents a numerically/Boolean valued function defined on the mesh. More precisely, the first argument in initialization can either be `‘int’`

, `‘size_t’`

, or `‘double’`

, which denotes the integers, unsigned integers and floats, respectively. The Boolean valued function is initialized with the `‘bool’`

argument.

```
[12]:
```

```
boundary = fe.MeshFunction("size_t", mesh, mesh.topology().dim()-1,0)
for v in fe.facets(mesh):
if fe.near(v.midpoint()[0], 0.0):
boundary[v] = 1 # left boundary
elif fe.near(v.midpoint()[0], L):
boundary[v] = 2 # right boundary
```

The previous `for`

loop iterates through facets of the mesh (vertices in our example) and assigns the value one to the left boundary and the value two to the right boundary of the bar.

note: it is preferable to set some default numerical value (e.g., zero) to all the other mesh facets, this can be done either by giving the last optional argument`0`

in the initialization of the`MeshFunction`

object or by calling`boundary.set_all(0)`

beforethe`for`

loop

#### Dirichlet boundary conditions¶

Dirichlet boundary conditions are sometimes called essential, since they determine the definition of the function spaces. In FEniCS, Dirichlet boundary conditions are implemented through the object `DirichletBC`

initialized with three arguments: the function space, the prescribed value at the boundary, and the corresponding subdomain specified by the previously defined `left`

function (or equivalently by the `boundary`

variable and its numerical value at the desired facets).

```
[13]:
```

```
# --------------------
# Boundary conditions
# --------------------
bc = fe.DirichletBC(V, 0.0, left)
# Equivalently, using boundary variable:
#bc = fe.DirichletBC(V, 0.0, boundary, 1)
```

#### Neumann boundary conditions¶

In the formulation of the problem, we also prescribed a normal force \(\mathbf{g}\) at the right boundary of the domain \(\Gamma_{\text{N}}\). Such boundary interactions falls into the category of Neumann boundary condition, which is sometimes denoted as natural. The natural boundary conditions are reflected directly in the weak formulation of the problem.

In the weak formulation above, the Neumann boundary condition appears as a term from the integration by parts formula. To include this term in our implementation of the problem, we first need to define the integrating measures that correspond to the manifolds representing the geometry of the mesh. To this end, we initialize the `Measure`

object with input arguments that specify the integral type and the respective geometry to which the measure corresponds.

```
[14]:
```

```
dx = fe.Measure("dx", mesh)
ds = fe.Measure("ds", subdomain_data = boundary)
```

### Weak formulation¶

We are now ready to give the weak form of the problem. One of the very nice features of FEniCS is that the implementation closely imitates the weak formulation given above.

```
[15]:
```

```
# --------------------
# Weak form
# --------------------
#a = E*A*fe.inner(fe.grad(u_tr),fe.grad(u_test))*dx
# One can equivalently write the bilinear form in this one-dimensional problem as
a = E*A*u_tr.dx(0)*u_test.dx(0)*fe.dx
l = b*u_test*dx + g*u_test*ds(2)
```

We need to supply the subdomain measure `ds`

with the value of the `MeshFunction`

object that corresponds to the Neumann part of the boundary \(\Gamma_{\text{N}}\), i.e., two.

note: it is possible to access individual elements of the gradient`grad`

of a function`u`

via`u.dx(i)`

, where`i`

specifies the derivative with respect to the`i`

-th coordinate and the enumeration starts withzero, do not confuse this with the measure object`dx`

that specifies the integration domain

`dot(u,v)`

vs`inner(u,v)`

: FEniCS makes distinction between the inner product and dot product of two elements`u`

,`v`

,`dot(u,v)`

contracts (sums) over the last index of the first element and the first index of the second element, whereas`inner(u,v)`

sums over all indices of the elements that must be of the same order

- if
`u`

,`v`

are both vectors (rank 1 tensor), then`dot(u,v)=inner(u,v)`

- matrix-vector multiplication \(\sigma \mathbf{n}\) is computed using
`dot(sigma,n)`

Finally, we call the solver to compute the linear algebraic system generated by the finite element method. FEniCS is prompted to do that by calling the `solve()`

function supplied with three parameters: linear variational equation, function that will store the resulting approximation and the problem specific Dirichlet boundary conditions.

Let us first introduce the function `u`

that will contain the finite element solution. In the code environment we define the object `Function`

pertaining to the function space \(V\). Then we can call `solve(a==l, u, bc)`

.

```
[16]:
```

```
# --------------------
# Solver
# --------------------
u = fe.Function(V)
fe.solve(a == l, u, bc)
# Equivalent implementation:
# F = a-l
# problem = fe.LinearVariationalProblem(fe.lhs(F),fe.rhs(F),u,bc)
# solver = fe.LinearVariationalSolver(problem)
```

## Post processing¶

In the final stage of our simulation we visualize the computational results. This task is easily handled by the `matplotlib`

library that we imported at the beginning.

For many trivial problems, ours included, there is an exact solution. We can use this solution to benchmark our numerical results. The following cell introduces the exact solution and plots it next to the approximative result given by the finite element method.

```
[17]:
```

```
# --------------------
# Exact solution
# --------------------
x_ex = np.linspace(0, L, 10)
u_ex = [-0.5*b*x_ex_i**2/E/A + (g+b*L)*x_ex_i/E/A for x_ex_i in x_ex]
# --------------------
# Post-process
# --------------------
fe.plot(u)
plt.plot(x_ex, u_ex, "x")
plt.xlabel("x [m]")
plt.ylabel("u [m]")
plt.legend(["FEM solution","exact solution"])
plt.show()
```

## Complete code¶

```
[ ]:
```

```
import fenics as fe
import matplotlib.pyplot as plt
import numpy as np
# --------------------
# Parameters
# --------------------
E = 0.025e9 # Young's modulus
A = np.pi*1e-4 # Cross-section area of bar
L = 0.05 # Length of bar
n = 20 # Number of elements
b = 0.03 # Load intensity
g = 0.0005 # External force
# --------------------
# Geometry
# --------------------
mesh = fe.IntervalMesh(n, 0.0, L)
fe.plot(mesh)
plt.show()
# --------------------
# Function spaces
# --------------------
V = fe.FunctionSpace(mesh, "CG", 1)
u_tr = fe.TrialFunction(V)
u_test = fe.TestFunction(V)
# --------------------
# Boundary marking
# --------------------
boundary = fe.MeshFunction("size_t", mesh, mesh.topology().dim()-1,0)
for v in fe.facets(mesh):
if fe.near(v.midpoint()[0], 0.0):
boundary[v] = 1 # left boundary
elif fe.near(v.midpoint()[0], L):
boundary[v] = 2 # right boundary
dx = fe.Measure("dx", mesh)
ds = fe.Measure("ds", subdomain_data = boundary)
# --------------------
# Boundary conditions
# --------------------
bc = fe.DirichletBC(V, 0.0, left)
# --------------------
# Weak form
# --------------------
a = E*A*fe.inner(fe.grad(u_tr),fe.grad(u_test))*dx
l = b*u_test*dx + g*u_test*ds(2)
# --------------------
# Solver
# --------------------
u = fe.Function(V)
fe.solve(a == l, u, bc)
# --------------------
# Exact solution
# --------------------
x_ex = np.linspace(0, L, 10)
u_ex = [-0.5*b*x_ex_i**2/E/A + (g+b*L)*x_ex_i/E/A for x_ex_i in x_ex]
# --------------------
# Post-process
# --------------------
fe.plot(u)
plt.plot(x_ex, u_ex, "x")
plt.xlabel("x [m]")
plt.ylabel("u [m]")
plt.legend(["FEM solution","exact solution"])
plt.show()
```