An Elastically Supported Rod

An Elastically Supported Rod#

Part 1: A modification to the PDE: continuous elastic support#

In the book, the finite element derivation and implementation of rod extension (the 1D Poisson equation) is presented. In this workshop, you are asked to do the same for a slightly different problem.

For this exercise we consider an elastically supported rod. An example of this would be a foundation pile in soil. The equilibrium of an elastically supported rod can be described with the following differential equation:

\[ -EA \frac{\partial^2 u}{\partial x^2} + ku = q \]

where \(EA\) is the axial stiffness of the rod, \(u(x)\) is the unknown displacement, \(k\) is the stiffness of the elastic support and \(q\) is a distributed load. We consider a case with Dirichlet boundary condition at \(x=0\) and a Neumann boundary condition at \(x=L\).

\[\begin{split} u = 0, \quad \text{at} \quad x = 0 \\ EA\frac{\partial u}{{\partial x}} = F, \quad \text{at} \quad x = L \end{split}\]

This differential equation is the inhomogeneous Helmholtz equation, which also has applications in dynamics and electromagnetics.

The finite element discretized version of this PDE can be obtained following the same steps as shown for the unsupported rod in the book. The additional term with respect to the case without elastic support is the second term on the left hand side: \(ku\). Using Neumann boundary condition (i.e. an applied load) at \(x=L\), the following expression is found for the discretized form:

\[\left[\int \mathbf{B}^T EA \mathbf{B} + \mathbf{N}^T k \mathbf{N} \,dx\right]\mathbf{u} = \int \mathbf{N}^T q \,d x + \mathbf{N}^T F \Bigg|_{x=L} \]

Task 1.1: Derive the discrete form

Derive the discrete form of the PDE given above. You can follow the same steps as in the book for the term with \(EA\) and the right hand side, but now carrying along the additional term \(ku\) from the PDE. Note that there are no derivatives in this term \(ku\) which means that integration by parts does not need to be applied on it. Show that this term leads to the \(\int\mathbf{N}^Tk\mathbf{N}\,dx\) term in the \(\mathbf{K}\)-matrix.

Part 2: Modification to the FE implementation#

The only change with respect to the procedure as implemented in the book is the formulation of the \(\mathbf{K}\)-matrix, which now consists of two terms:

\[ \mathbf{K} = \int \mathbf{B}^TEA\mathbf{B} + \mathbf{N}^Tk\mathbf{N}\,dx \]

To calculate the integral exactly we must use two integration points.

\[ \mathbf{K_e} = \sum_{i=1}^{n_\mathrm{ip}} \left(\mathbf{B}^T(x_i)EA\mathbf{B}(x_i) + \mathbf{N}^T(x_i) k\mathbf{N}(x_i) \right) w_i\]

Task 2.1: Code implementation

The only change needed with respect to the implementation of the book is in the calculation of the element stiffness matrix. Copy the code from the book and add the term related to the distributed support in the right position.

Remarks:

  • The function evaluate_N is already present in the code in the book

  • The get_element_matrix function already included a loop over two integration points

  • You need to define \(k\) somewhere. To allow for varying \(k\) as required below, it is convenient to make \(k\) an additional argument of the simulate function and pass it on to lower level functions from there (cf. how \(EA\) is passed on)

import numpy as np
import matplotlib.pyplot as plt

### YOUR_CODE_HERE ###

Task 2.2: Analyze the influence of \(k\) on the results

Use the following parameters: \(L=3\text{ m}\), \(EA=1000\text{ N}\), \(F=10\text{ N}\), \(q=0 \text{ N/m}\) (all values are the same as in the book, except for \(q\)). Additionally, use \(k=1000\) N/m2 and a discretization of 10 elements.

Check the influence of the distributed support on the solution:

  • First use \(q=0\) N/m and \(k=1000\) N/m2

  • Then set \(k\) to zero and compare the results

  • Does the influence of the distributed support on the solution make sense?

# use the simulate function to get displacements
### YOUR CODE HERE ###

# plot u versus x for both cases
### YOUR CODE HERE ###

Task 2.3: Investigate the influence of discretization on the quality of the solution

  • How many elements do you need to get a good solution with \(k=1000\) N/m?

  • How about when the stiffness of the distributed support is increased to \(k=10^6\) N/m2 ?

Simulate and plot results for different numbers of elements for each of the above questions.

### YOUR_CODE_HERE ###

Part 3: Change boundary condition#

Now suppose the end of the rod at \(x=L\) is loaded not by a force but by a prescribed displacement. What do you need to change to the simulate function to achieve this? Hint: look at the box Eliminiation of Dirichlet boundary conditions on the Finite element implementation page in the book for a description of how to deal with inhomogeneous Dirichlet boundary conditions in a more general case.

Task 3.1: Implement changes required to prescribe \(u(L)\)

  • Code the simulate_new function below such that a nonzero Dirichlet boundary condition is applied at \(x=L\) instead of a force

  • Perform the simulation with the same parameters as before and \(k=1000\) and \(u_\mathrm{right}=0.01\) m.

  • Compare the obtained field \(u(x)\) to that from task 2.2.

def simulate_new(n_element,
                 length=3,
                 EA=1e3,
                 k=1000,
                 u_right=0.01,
                 q_load=0
    ):
### YOUR CODE HERE ###

    
x,u = simulate_new(10, q_load=0, u_right=0.01)
plt.plot(x,u)
plt.xlabel('x [m]')
plt.ylabel('u [m]')

By Frans van der Meer, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.

Contents