1D Problems

1D Problems#

Poisson#

The poisson equation is often solved using finite elements. The poisson equation can be used to describe a variety of physical phenomena such as steady-state heat conduction and flow in permeable media (Darcy’s law). The constitutive relationship depends on the problem being adressed.

The Poisson equation, in its most general form, is defined as follows;

()#\[ ∇² u = -\rho \]

where

$∇²u$ represents the Laplacian of a scalar field u, which is the second spatial derivative of $u$
$\rho$ represents a source or sink term that accounts for the distribution of mass or charge within the system. It can be positive (sources) or negative (sinks) and can vary with spatial coordinates.
rod
soil consolidation

Alternatively, the poisson equation in a continuous body Ω can expressed as follows:

()#\[ ∇² q + f = 0 \]

where q is a flux vector and $f$ is a source term.

Note: In the case of heat conductivity, f would be the heat flux vector. In Darcy's law q would be the flow rate.

The constitutive relationship is given by:

()#\[ q = \kappa ∇ u \]

where u is a potential.

Inserting the constitutive equation into the equation 2.1 yields a poisson equation. As mentioned at the beggining of this chapter, the constitutive relation depends on the the problem being solved. For linear heat conduction, the scalar u is the temperature. For an isotropic medium

κ = $κ$ I

where $\kappa$ is the thermal conductivity.

Note: For Darcy's law u is the hydraulic head and κ is known as the hydraulic conductivity.

Intermezzo: Let's consider the continuous body Ω in the following figure

{figure} .././images/Chapter3/3_1_1.png
---
height: 400px

name: 3_1_1
---
Continuous body 

On part of the boundary, the displacements are prescribed (Dirichlet boundary Conditions). In the figure, this part is denoted as $Γ_g$

u = g on $Γ_g$

On the part of the boundary denoted as $Γ_h$, tractions are applied, which could also be zero.

Traction is the force per unit area on a surface

t = h on $Γ_h$

Back in the poisson problem, the boundary value problem requires boundary conditions. The dirichlet conditions impose

u = g on $Γ_g$

which prescribes the potential on the boundary.

Note: For heat conduction this is the temperature on the boundary $Γ_g$. For Darcy flow this is the hydraulic head.

The Neumann boundary condition requires that -q n = h on $Γ_h$

where n is the outward normal to the surface $Γ_h$

Note: In the case of heat conduction, this imposes the heat flux on $Γ_h$.

For Darcy flow it is the inflow/ outflow across a boundary.

The weak form of the Poisson equation is derived as follows: Multiplying equation (2.1) by a weight function w and integrating over the body Ω.

()#\[ \int_{Ω} w( ∇· \textbf{q})dΩ +\int_{Ω} w fdΩ =0, ∀ w ∈ V \]

Integrating the LHS of the above equation by parts yields:

()#\[ \int_{Ω} ∇w· \textbf{q}dΩ - \int_{Ω} w \textbf{q}· \textbf{n} dΩ +\int_{Ω} w fdΩ =0, ∀ w ∈ V \]

Inserting now the constitutive relationship and the boundary condition − q · n = h leads to the problem: find u ∈ S such that:

()#\[ -\int_{Ω} ∇w ·\textbf{κ} ∇qdΩ + \int_{Γ_h} w h dΓ +\int_{Ω} w fdΩ =0, ∀ w ∈ V \]

which is the weak form of the Poisson equation.

Beam problems#

Typically two types of beams are addressed in literature, the Bernoulli - Euler Beam and the Timoshenko beam. The the Bernoulli - Euler Beam is valid for relatively slender beams, whereas the Timoshenko beam is relevant where shear deformations are present. In order to solve beam elements, a similar process is followed to the one used in continuum elements. Again, since the governing equation is formulated, a weak form of the equation must be developed. At this stage only in-plane, linear beams are considered. All forces and moments are in-plane.

Kinematics of a beam: Beam theory is based upon the assumption that planes which are normal to the beam’s axis remain plane. For a transverse displacement v Figure 2.1 (left panel) shows the rotation of the plane due to a rigid body rotation.

Considering now a fibre which is initially perpendicular to the beam axis. Subjecting a segment from a beam to pure shear, the fibre will not rotate, as illustrated in Figure 2.1 (right panel). Relative to the fibre which has not rotated, a plane which remains perpendicular to the axis undergoes a rotation γ.

For a beam segment subjected to both rotation and shear deformation, the rotation θ of a fibre which is initially perpendicular to the beam’s axis. $$ \theta = v_{,x} - \gamma $$ (beam_kinematics)

{figure} .././images/Chapter3/3_1_2.png
---
height: 400px
name: 3_1_2
---
Plane Beam Element Ω= (x1,x2)

Equilibrium of a beam#

The equilibrium of a beam can be developed in two ways: either elaborating on the kinematics of the beam or considering the equilibrium of a beam directly.

The bending moment $m$ in a beam is defined as

()#\[ m = \int_{-h/2}^{h/2} \sigma_{xx} \: y \: dy \]

and the shear stress $q$ defined by

()#\[ q = \int_{-h/2}^{h/2} \sigma_{xy} \: dy \]

The resultant force $Q$ and moment $M$ can be related to $q$ and $m$ through the equations $Q = q n$ and $M = m n$ where $n$ is the outward unit normal vector of the element.

Considering the vertical equilibrium of a beam element, it follows that

()#\[ \int_{d\Omega} q \: n \: d\Gamma + \int_{\Omega} f_y \: d\Omega = 0 \]

Following that $\int_{d\Omega} q \: n \: d\Gamma = q|_{x=L_2} - q|_{x=L_1}$, this gives

()#\[ \int_{\Omega} q_{,x} \: d\Omega + \int_{\Omega} f_y \: d\Omega = 0 \]

Since the equilibrium described by Equation beam_eq2 must hold for an infinitely small segment of a beam, it must hold that

()#\[ q_{,x} + f_y = 0 \]

For rotational equlibrium, it follows that

()#\[ \int_{d\Omega} m \: n \: d\Gamma - \int_{d\Omega} q \: n \: x \: d\Gamma - \int_{\Omega} f_y \: x \: d\Omega = 0 \]

which can be rearranged such that

()#\[ \int_{\Omega} m_{,x} \: d\Omega - \int_{\Omega} q \: d\Omega - \int_{\Omega} q_{,x} \: x \: d\Omega - \int_{\Omega} f_y \: x \: d\Omega = 0 \]

Since satisfaction of the translational equilibrium implies that $\int_{\Omega} (q_{,x} + f_y) \: x \: d\Omega = 0$, the rotational equlibrium implies that

()#\[ m_{,x} - q = 0 \]

…

Euler-Bernoulli beam#

Strong form equation

The Euler-Bernoulli beam does not allow for shear deformation ($\gamma=0$). As a result, the rotation of the beam can be directly related to the displacement $v$:

()#\[ \theta = \frac{\text{d}v}{\text{d}x} \]

which means that

()#\[ \kappa = \frac{\text{d}^2v}{\text{d}x^2} \]

and the bending moment $m$ in the beam is determined by:

()#\[ m=-EI \cdot \kappa=-EI \cdot \frac{\text{d}^2v}{\text{d}x^2} \]

where

$EI$ is the bending stiffness of the beam.

Taking the derivative of all terms in Equation beam_rot_eq3 to $x$ and substituting in beam_tran_eq3 yields

()#\[ \frac{\partial^2 m}{\partial x^2} + f_y = 0 \]

Assuming $EI$ to be constant and substiuting in bending moment yields

()#\[ -EI \frac{\partial^4 v}{\partial x^4} + f_y = 0 \]

which is the strong form equation of equilibrium for an Euler-Bernoulli beam. Being a fourth-order equation, two boundary conditions at both ends of the beam. Dirichlet boundary conditions involved the prescription of the displacement or rotation and Neumann involves either the shear force or moment. With appropriate boundary conditions, the boundary value problem is complete and can be solved.

Weak governing equation

Following procedures from (REF to weak form eq.), the weak form of equilibrium for a beam can be developed. Multiplying Equation EB_EOM_eq1 by a weight function $\overline{v}$, from an appropriately defined space, which is equal to zero where Dirichlet boundary conditions are applied and integration over the beam $\Omega$ yields:

()#\[ \int_{\Omega} \overline{v} \: m_{,xx} \: d\Omega + \int_{\Omega} \overline{v} f_y \: d\Omega = 0 \]

Integrating by parts the term involving the moment $M$ yields:

()#\[ - \int_{\Omega} \overline{v}_{,x} \: m_{,x} \: d\Omega + \int_{\Gamma} \overline{v} \: m_{,x} \: n \: d\Gamma + \int_{\Omega} \overline{v} \: f_y \: d\Omega = 0 \]

Applying integration by parts again, this time to the term $\int_{\Omega} \overline{v}_{,x} \: m_{,x} \: d\Omega$, yields

()#\[ \int_{\Omega} \overline{v}_{,xx} \: m \: d\Omega - \int_{\Gamma} \overline{v}_{,x} \: m \: n \: d\Gamma + \int_{\Gamma} \overline{v} \: m_{,x} \: n \: d\Gamma + \int_{\Omega} \overline{v} \: f_y \: d\Omega = 0 \]

Inserting now the Neumann boundary conditions and the consitutive relation from Equation bending moment, solving the governing weak form equation for a beam involves: find $v \in S$ such that

()#\[ -\int_{\Omega} \overline{v}_{,xx} \: EI \: v_{,xx} \: d\Omega - \int_{\Gamma_M} \overline{v}_{,x} \: T \: d\Gamma + \int_{\Gamma_Q} \overline{v} \: f_y \: d\Omega + \int_{\Omega} \overline{v} \: f_y \: d\Omega = 0 \]

where $S$ and $V$ are approprately defined spaces.

…

Timoshenko beam#

Strong form equation

The Timoshenko beam allows shear deformation and is more suitable for relatively short beams. In Timoshenko beam theory, planes remain plain during deformation but are not needed to remain normal to the longitudinal axis. This assumption allows for shear deformation. The plane is assumed to have a rotation $\theta$ of aplane given by the equation

()#\[ \theta = \frac{dv}{dx} - \gamma \]

where $\gamma$ is the rotation due to shear. In terms of the shear force and properties of the beam, it is related to the shear force by the following constitutive relationship:

()#\[ \gamma = \frac{q}{GA_s} \]

where $G$ is the shear modulus and $A_s$ is the effective shear area. Relative to the Euler-Bernoulli theory, an additional unknown, the shear strain $\gamma$, has been introduced. The governing equations are given by the equilibrium conditions in equations beam_rot_eq3 and beam_tran_eq3 with the following constitutive relationships:

()#\[ m = -EI\kappa \]

()#\[ q = GA_s \gamma = GA_s (\frac{dv}{dx} - \theta) \]

and considering $v$ and $\theta$ to be the fundamental unknowns. Inserting the constitutive relationships into the equilibrium equations beam_rot_eq3 and beam_tran_eq3 leads to:

()#\[ -EI \frac{d^2 \theta}{dx^2} - GA_s (\frac{dv}{dx} - \theta) = 0 \]

()#\[ GA_s (\frac{d^2 v}{dx^2} - \frac{d\theta}{dx}) + f_y = 0 \]

which are two coupled second-order equations, the first for rotation $\theta$, and the second for the displacement $v$, in the space $\Omega$. Together with the previously defined boundary conditions, this problem is complete.

Weak governing equations

Multiplying equations TS_STRONG_eq1 and TS_STRONG_eq2 by appropriately defined weight functions $\overline{\theta}$ and $\overline{v}$,

()#\[ - \int_{\Omega} \overline{\theta} EI \theta_{,xx} \: d\Omega - \int_{\Omega} \overline{\theta} GA_s (v_{,x} - \theta) \: d\Omega = 0 \]

()#\[ \int_{\Omega} \overline{v} GA_s (v_{,xx} - \theta_{,x}) \: d\Omega + \int_{\Omega} \overline{v} f_y \: d\Omega = 0 \]

Applying integrations by parts once gives

()#\[ \int_{\Omega} \overline{\theta}_{,x} EI \theta_{,x} \: d\Omega - \int_{\Omega} \overline{\theta} GA_s (v_{,x} - \theta) \: d\Omega + \int_{\Gamma_M} \overline{\theta} T \: d\Gamma = 0 \]

()#\[ - \int_{\Omega} \overline{v}_{,x} GA_s (v_{,x} - \theta) \: d\Omega + \int_{\Gamma} \overline{v} GA_s (v_{,x} - \theta) n_x \: d\Gamma + \int_{\Omega} \overline{v} f_y \: d\Omega = 0 \]

and then inserting the Neumann boundary conditions yields the weak problem for a Timoshenko beam of: find $v \in S_v$ and $\theta \in S_{\theta}$ such that

()#\[ \int_{\Omega} \overline{\theta}_{,x} EI \theta_{,x} \: d\Omega - \int_{\Omega} \overline{\theta} GA_s (v_{,x} - \theta) \: d\Omega + \int_{\Gamma_M} \overline{\theta} T \: d\Gamma = 0 \]

()#\[ - \int_{\Omega} \overline{v}_{,x} GA_s (v_{,x} - \theta) \: d\Omega + \int_{\Gamma_Q} \overline{v} F_y \: d\Gamma + \int_{\Omega} \overline{v} f_y \: d\Omega = 0 \]

A more precise definition of the necessary functions spaces should be provided to complete the problem.

Finite element formulations for beams#

Euler beams

A Galerkin problem for a Bernoulli-Euler beam involves: find $v^h \in S^h$ such that:

()#\[ \int_{\Omega} \overline{v}_{,xx}^{h} EI v_{,xx}^h \: d\Omega = - \int_{\Gamma_M} \overline{v}_{,x}^h T \: d\Gamma + \int_{\Gamma_Q} \overline{v}^h F_y \: d\Gamma + \int_{\Omega} \overline{v}^h f_y \: d\Omega = 0 \]

where $S^h \subset S$ and $V^h \subset V$ are finite-dimensional spaces. As for continuum elements, we wish to express the displacement field $v$ in terms of shape functions and nodal degrees of freedom. The problem tat arises is that simple $C^0$ finite element shape functions are not suitable. The above equation requires the evluation of the second derivative of the interpolated displacement field. However the second derivative of a $C^0$ continuous function does not exist in a classical sense. The use of $C^0$ interpolations for fourth-oder problems is not mathimatically consistent and can lead to unpredictable results. Crucially, convergence of the solution is not assured for interpolations with insufficient continuity.

The solution is to use $C^1$ shape functions. Such shape functions can be constructed relatively easily in one dimension (the extension to multiple dimensions is however far from trivial). Hermitian polynomials are $C^1$ functions, and involve both displacement and rotational degrees of freedom. A Hermitian beam element with two nodes has four degrees of freedom (two displacement degrees of freedom and two rotation degrees of freedom). This results in a cubic interpolation of the displacement along the element. The displacement field $v^h$ is given by:

()#\[ v^h (x) = \sum_{i}^{2} (N_i (x) v_i + M_i (x) \theta_i) \]

where $v_i$ and $\theta_i$ are degrees of freedom associated with node $i$, and the shape functions are equal to:

()#\[ N_1 = \frac{-(x-x_2)^2 (-h+2(x_1-x))}{h^3} \]

()#\[ N_2 = \frac{(x-x_1)^2 (h+2(x_2 - x))}{h^3} \]

()#\[ M_1 = \frac{(x-x_1)(x-x_2)^2}{h^2} \]

()#\[ M_2 = \frac{(x-x_1)^2 (x-x_2)}{h^2} \]

for an element of length $h$ with ends from $x_1$ to $x_2$ ($x_2 > x_1$). The displacement at a point in the beam is given by:

()#\[\begin{split} v^h = \textit{\textbf{N}} \textit{\textbf{a}}_e = \ \begin{bmatrix} N_1 & M_1 & N_2 & M_2 \end{bmatrix} \begin{Bmatrix} v_1 \\ \theta_1 \\ v_2 \\ \theta_2 \end{Bmatrix} \end{split}\]

It is necessary to compute both the first and second deriovatives of $v$ with respect to $x$. The first derivative is given by:

()#\[\begin{split} v^h_{,x} = \textit{\textbf{N}}_{,x} \textit{\textbf{a}}_e = \ \begin{bmatrix} N_{1,x} & M_{1,x} & N_{2,x} & M_{2,x} \end{bmatrix} \begin{Bmatrix} v_1 \\ \theta_1 \\ v_2 \\ \theta_2 \end{Bmatrix} \end{split}\]

and the second derivative can be calculated through:

()#\[\begin{split} v^h_{,xx} = \textit{\textbf{N}}_{,xx} \textit{\textbf{a}}_e = \ \begin{bmatrix} N_{1,xx} & M_{1,xx} & N_{2,xx} & M_{2,xx} \end{bmatrix} \begin{Bmatrix} v_1 \\ \theta_1 \\ v_2 \\ \theta_2 \end{Bmatrix} \end{split}\]

Now that $v^h$, $v^h_{,x}$ and $v^h_{,xx}$ can be computed given $\textit{\textbf{a}}_e$, they can be inserted into the Galerkin problem,

()#\[ \int_{\Omega} (\textit{\textbf{N}}_{,xx} \textit{\textbf{b}}_e)^T EI \textit{\textbf{N}}_{,xx} \textit{\textbf{a}}_e \: d\Omega = -\int_{\Gamma_M} (\textit{\textbf{N}}_{,x} \textit{\textbf{b}}_e)^T T \: d\Gamma + \int_{\Gamma_Q} (\textit{\textbf{N}} \textit{\textbf{b}}_e)^T F_y \: d\Gamma + \int_{\Omega} (\textit{\textbf{N}} \textit{\textbf{b}}_e)^T f_y \: d\Omega \]

After some rearranging,

()#\[ \int_{\Omega} \textit{\textbf{N}}_{,xx}^T EI \textit{\textbf{N}}_{,xx} \: d\Omega \: \textit{\textbf{a}}_e = - \int_{\Gamma_M} \textit{\textbf{N}}_{,x}^T T \: d\Gamma + \int_{\Gamma_Q} \textit{\textbf{N}}^T F_y \: d\Gamma + \int_{\Omega} \textit{\textbf{N}}^T f_y \: d\Omega \]

The element stiffness matrix is then given by:

()#\[ \textit{\textbf{k}}_e = \int_{\Omega} \textit{\textbf{N}}_{,xx}^T EI \textit{\textbf{N}}_{,xx} \: d\Omega \]

and the RHS vector is given by:

()#\[ \textit{\textbf{f}}_e = \int_{\Gamma_Q} \textit{\textbf{N}}^T F_y \: d\Gamma - \int_{\Gamma_M} \textit{\textbf{N}}_{,x}^T T \: d\Gamma + \int_{\Omega} \textit{\textbf{N}}^T f_y \: d\Omega \]

The operation to form the RHS vector essentially translates the applied loads into equivalent nodal shear forces and moments. Taking derivatives of the Hermitian shape functions in FE_EB_N1 until FE_EB_M2:

1D elements in 2/3D (space-frame structures)#

Here, a one dimensional rod element is considered in two and three dimensions. In the 2D case the nodes of the truss can translate in both x and y dimensions. As a result, a finite element model will require two degrees of freedom at each of the nodes. Typically, the coordinate system is selected considering that the x-axis is aligned with the axis of a truss element. In the 3D case, every node has three degrees of freedom.

See Garth Wells / Chapter 6.1. Rod in space (Truss) ?