xbe.at

Python / Scientific and math computing

Finite Element Analysis with Python:

5 min read, 15 sections

Contents
  1. Introduction to Finite Element Analysis (FEA)
  2. Setting Up the FEA Environment
  3. Defining the Problem Domain
  4. Element Stiffness Matrix
  5. Global Stiffness Matrix Assembly
  6. Applying Boundary Conditions
  7. Solving the System
  8. Visualizing the Results
  9. Real-life Example 1 - Beam Deflection
  10. Element Types in FEA
  11. Nonlinear FEA
  12. Real-life Example 2 - Stress Analysis
  13. Mesh Refinement
  14. Error Estimation and Convergence
  15. Additional Resources

Introduction to Finite Element Analysis (FEA)

Finite Element Analysis is a numerical method for solving complex engineering problems. It involves dividing a large problem into smaller, simpler parts called finite elements. Let's visualize a simple mesh generation:

import numpy as np
import matplotlib.pyplot as plt

def create_mesh(x_range, y_range, num_elements):
    x = np.linspace(x_range[0], x_range[1], num_elements + 1)
    y = np.linspace(y_range[0], y_range[1], num_elements + 1)
    return np.meshgrid(x, y)

x, y = create_mesh((0, 1), (0, 1), 10)
plt.figure(figsize=(8, 6))
plt.plot(x, y, 'b')
plt.plot(x.T, y.T, 'b')
plt.title('Simple 2D Mesh for FEA')
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.show()

Setting Up the FEA Environment

Before diving into FEA, we need to set up our Python environment with necessary libraries:

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt

# Example: Check versions
print(f"NumPy version: {np.__version__}")
print(f"SciPy version: {linalg.__version__}")
print(f"Matplotlib version: {plt.__version__}")

Defining the Problem Domain

Let's consider a simple 1D heat conduction problem. We'll define the domain and boundary conditions:

# Problem parameters
L = 1.0  # Length of the domain
k = 1.0  # Thermal conductivity
T_left = 100.0  # Temperature at left boundary
T_right = 0.0  # Temperature at right boundary

# Discretization
n_elements = 10
element_length = L / n_elements

# Node coordinates
x = np.linspace(0, L, n_elements + 1)

Element Stiffness Matrix

For our 1D heat conduction problem, let's create the element stiffness matrix:

def element_stiffness_matrix(k, l):
    return k / l * np.array([[1, -1], [-1, 1]])

K_e = element_stiffness_matrix(k, element_length)
print("Element Stiffness Matrix:")
print(K_e)

Global Stiffness Matrix Assembly

Now, we'll assemble the global stiffness matrix:

def assemble_global_stiffness(n_elements, K_e):
    K = np.zeros((n_elements + 1, n_elements + 1))
    for i in range(n_elements):
        K[i:i+2, i:i+2] += K_e
    return K

K = assemble_global_stiffness(n_elements, K_e)
print("Global Stiffness Matrix:")
print(K)

Applying Boundary Conditions

Let's apply the Dirichlet boundary conditions to our problem:

def apply_boundary_conditions(K, T_left, T_right):
    K_bc = K.()
    K_bc[0, :] = 0
    K_bc[0, 0] = 1
    K_bc[-1, :] = 0
    K_bc[-1, -1] = 1
    
    F = np.zeros(K.shape[0])
    F[0] = T_left
    F[-1] = T_right
    
    return K_bc, F

K_bc, F = apply_boundary_conditions(K, T_left, T_right)

Solving the System

Now we can solve the system of equations:

def solve_system(K, F):
    return linalg.solve(K, F)

T = solve_system(K_bc, F)
print("Temperature distribution:")
print(T)

Visualizing the Results

Let's plot the temperature distribution:

plt.figure(figsize=(10, 6))
plt.plot(x, T, 'ro-')
plt.title('Temperature Distribution')
plt.xlabel('Position (m)')
plt.ylabel('Temperature (°C)')
plt.grid(True)
plt.show()

Real-life Example 1 - Beam Deflection

Let's consider a simple beam deflection problem:

def beam_deflection(L, E, I, w):
    x = np.linspace(0, L, 100)
    y = w * x**2 * (L - x)**2 / (24 * E * I)
    return x, y

L = 5  # Length of beam (m)
E = 200e9  # Young's modulus (Pa)
I = 3e-5  # Moment of inertia (m^4)
w = 10000  # Distributed load (N/m)

x, y = beam_deflection(L, E, I, w)

plt.figure(figsize=(10, 6))
plt.plot(x, y)
plt.title('Beam Deflection')
plt.xlabel('Position along beam (m)')
plt.ylabel('Deflection (m)')
plt.grid(True)
plt.show()

Element Types in FEA

Different element types are used for various problems. Let's visualize some common 2D elements:

def plot_element(vertices, element_type):
    plt.figure(figsize=(4, 4))
    plt.plot(vertices[:, 0], vertices[:, 1], 'ro-')
    plt.title(f'{element_type} Element')
    plt.axis('equal')
    plt.grid(True)
    plt.show()

# Triangle element
triangle = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3)/2]])
plot_element(triangle, 'Triangle')

# Quadrilateral element
quad = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
plot_element(quad, 'Quadrilateral')

Nonlinear FEA

Nonlinear FEA is used for problems with material nonlinearity, geometric nonlinearity, or both. Here's a simple example of a nonlinear spring:

def nonlinear_spring(k, x):
    return k * x + 0.1 * k * x**3

x = np.linspace(-2, 2, 100)
F = nonlinear_spring(1000, x)

plt.figure(figsize=(8, 6))
plt.plot(x, F)
plt.title('Nonlinear Spring Force-Displacement Curve')
plt.xlabel('Displacement')
plt.ylabel('Force')
plt.grid(True)
plt.show()

Real-life Example 2 - Stress Analysis

Let's simulate a simple stress analysis for a plate with a hole:

import numpy as np
import matplotlib.pyplot as plt

def stress_concentration(r, a, sigma):
    theta = np.linspace(0, 2*np.pi, 100)
    x = r * np.cos(theta)
    y = r * np.sin(theta)
    
    sigma_xx = sigma * (1 - (a**2 / r**2) * (3/2 * np.cos(2*theta) + np.cos(4*theta)) + 
                        (3*a**4 / (2*r**4)) * np.cos(4*theta))
    
    return x, y, sigma_xx

r = 1.0  # Radius of the plate
a = 0.2  # Radius of the hole
sigma = 100  # Applied stress

x, y, sigma_xx = stress_concentration(r, a, sigma)

plt.figure(figsize=(10, 8))
sc = plt.scatter(x, y, c=sigma_xx, cmap='jet')
plt.colorbar(sc, label='Stress')
plt.title('Stress Distribution Around a Circular Hole')
plt.xlabel('X')
plt.ylabel('Y')
plt.axis('equal')
plt.show()

Mesh Refinement

Mesh refinement is crucial for accurate FEA results. Let's visualize different mesh densities:

def create_refined_mesh(x_range, y_range, num_elements):
    x = np.linspace(x_range[0], x_range[1], num_elements + 1)
    y = np.linspace(y_range[0], y_range[1], num_elements + 1)
    return np.meshgrid(x, y)

fig, axs = plt.subplots(1, 3, figsize=(15, 5))
mesh_sizes = [5, 10, 20]

for i, size in enumerate(mesh_sizes):
    x, y = create_refined_mesh((0, 1), (0, 1), size)
    axs[i].plot(x, y, 'b')
    axs[i].plot(x.T, y.T, 'b')
    axs[i].set_title(f'Mesh with {size}x{size} elements')
    axs[i].set_aspect('equal')

plt.tight_layout()
plt.show()

Error Estimation and Convergence

Error estimation is important in FEA. Let's simulate convergence for our 1D heat conduction problem:

def analytical_solution(x, L, T_left, T_right):
    return T_left + (T_right - T_left) * x / L

def compute_error(numerical, analytical):
    return np.max(np.abs(numerical - analytical))

errors = []
element_counts = [5, 10, 20, 40, 80]

for n in element_counts:
    x = np.linspace(0, L, n + 1)
    K = assemble_global_stiffness(n, element_stiffness_matrix(k, L/n))
    K_bc, F = apply_boundary_conditions(K, T_left, T_right)
    T_num = solve_system(K_bc, F)
    T_ana = analytical_solution(x, L, T_left, T_right)
    errors.append(compute_error(T_num, T_ana))

plt.figure(figsize=(10, 6))
plt.loglog(element_counts, errors, 'bo-')
plt.title('Convergence of FEA Solution')
plt.xlabel('Number of Elements')
plt.ylabel('Maximum Error')
plt.grid(True)
plt.show()

Additional Resources

For further reading on Finite Element Analysis using Python, consider the following resources:

  1. "A Tutorial on the Implementation of Finite Element Methods using Python" by A. Sridhar et al. (arXiv:2005.07742)
  2. "Finite Element Analysis with Python" by R. Jaiman (arXiv:1910.04772)
  3. "Numerical Methods in Engineering with Python 3" by J. Kiusalaas (Cambridge University Press)

These resources provide in-depth explanations and more advanced examples of FEA implementation in Python.

More in Scientific and math computing

  1. Harmonic Function Theory with Python11 min

    Harmonic functions are solutions to Laplace's equation and play a crucial role in many areas of mathematics and physics. They are smooth functions with...

  2. Exponential Equations and Functions Using Python5 min

    Exponential equations involve variables in the exponent. They are crucial in modeling growth and decay processes in various fields like finance, biology, and...

  3. Homological Algebra in Python9 min

    Homological Algebra is a branch of abstract algebra that studies algebraic structures and their properties using techniques from homology theory. It provides...

  4. Exploring Lie Groups and Lie Algebras with Python7 min

    Lie groups are continuous symmetry groups that play a crucial role in physics and mathematics. They are named after Sophus Lie, who studied their properties...

  5. Introduction to Boolean Algebra with Python4 min

    Boolean algebra is a mathematical system that deals with binary values, True and False, or 1 and 0. It provides a set of rules and operations used in...

  6. Exploring Elliptic Curve Arithmetic with Python6 min

    Elliptic curves are algebraic structures with applications in cryptography and number theory. In this presentation, we'll explore their arithmetic using...