FEM stands for Finite Element Method. It is a numerical technique used to solve partial differential equations (PDEs) and boundary value problems. FEM divides the problem domain into a finite number of sub-domains, known as elements, and approximates the behavior of the PDEs within each element using simple functions.
These simple functions are called basis functions, and they allow for the representation of the behavior of the PDE within each element using a small number of parameters. By assembling the element-wise solutions, a global approximation of the solution to the PDE can be obtained.
FEM has a wide range of applications in many areas of engineering and physics, including structural mechanics, fluid mechanics, heat transfer, electromagnetics, and more. It is widely used in the design and analysis of structures, such as buildings, bridges, aircraft, and ships. FEM has become an essential tool for engineers and scientists, as it provides a powerful and flexible method for solving complex problems in many different fields.
The FEM involves the following steps:
- Discretization: The problem domain is discretized into a finite number of sub-domains, called elements, which are typically defined by a set of nodes.
- Element formulation: Within each element, the behavior of the PDE is approximated using simple functions, known as basis functions. The basis functions are typically polynomials or piecewise polynomials, and they are used to interpolate the solution within the element.
- Assembly: The element-wise solutions are assembled into a global approximation of the solution to the PDE.
- Solution: The global system of equations is solved to obtain the numerical solution to the PDE.
The FEM can be expressed mathematically as follows:
Discretization:
- The domain is divided into a finite number of elements:
- Ω = ⋃Ωᵢ (where Ω is the domain and Ωᵢ are the elements)
- The boundary of the domain is defined as ∂Ω.
Element formulation:
- Within each element, the solution is approximated as a linear combination of basis functions:
- u(x) ≈ ∑ᵢ uᵢ φᵢ(x)
- where uᵢ are the nodal values, φᵢ(x) are the basis functions, and x is the spatial location within the element.
- The basis functions are typically chosen to be polynomials or piecewise polynomials that satisfy certain continuity and differentiability conditions.
- The behavior of the PDE within each element is approximated using a weak form of the PDE that involves integrating over the element and multiplying by a test function:
- ∫ᵩ (∇u⋅∇v + fv) dx = ∫ᵩ fv dx
- where v is a test function that satisfies certain continuity and differentiability conditions, and f is the source term.
Assembly:
- The element-wise solutions are assembled into a global system of equations:
- Ku = f
- where K is the global stiffness matrix, u is the vector of nodal values, and f is the vector of source terms.
Solution:
- The global system of equations is solved to obtain the numerical solution to the PDE.
- This involves solving a system of linear equations, which can be done using a variety of numerical methods, such as Gaussian elimination, LU decomposition, or iterative methods like conjugate gradient or GMRES.
Example of solving a 1D heat conduction problem using the FEM in Python.
Problem statement: Consider a rod of length L = 1 meter, with thermal conductivity k = 1 W/m-K and heat generation rate q = 0 W/m3. The rod is initially at a uniform temperature of T = 0°C. At time t = 0, the left end of the rod is suddenly raised to a temperature of T = 100°C, while the right end is kept at T = 0°C. Find the temperature distribution within the rod at time t = 10 seconds.
Mathematical formulation: The heat conduction equation in 1D is given by:
∂T/∂t = k∂²T/∂x² + q
Using the FEM, we can write the weak form of the equation as:
∫ᵩ ρCp∂T/∂tφ + k∇T⋅∇φ dx = ∫ᵩ qφ dx
where ρ is the density, Cp is the specific heat capacity, φ is the test function, and the boundary conditions are:
T(0) = 100°C
T(L) = 0°C
We can discretize the domain into N elements, with each element i having a length of h = L/N. Within each element, we can approximate the temperature distribution as:
T(x) ≈ ∑ᵢ Tᵢφᵢ(x)
where Tᵢ are the nodal values, and φᵢ(x) are the basis functions.
We can then substitute this approximation into the weak form and integrate over each element to obtain the element-wise equations:
[Mᵢ]{dTᵢ}/{dt} + [Kᵢ]{Tᵢ} = {fᵢ}
where [Mᵢ] is the element mass matrix, [Kᵢ] is the element stiffness matrix, and {fᵢ} is the element force vector.
The global system of equations is then obtained by assembling the element-wise equations into a global mass matrix [M], stiffness matrix [K], and force vector {f}:
[M]{dT}/{dt} + [K]{T} = {f}
We can then solve this system of equations numerically using the forward Euler method, which is a first-order explicit method for time integration.These simple functions are called basis functions, and they allow for the representation of the behavior of the PDE within each element using a small number of parameters. By assembling the element-wise solutions, a global approximation of the solution to the PDE can be obtained.
import numpy as np
import matplotlib.pyplot as plt
# Define problem parameters
L = 1.0 # Length of rod (m)
k = 1.0 # Thermal conductivity (W/m-K)
q = 0.0 # Heat generation rate (W/m^3)
rho = 1.0 # Density (kg/m^3)
Cp = 1.0 # Specific heat capacity (J/kg-K)
T_left = 100.0 # Temperature at left end (°C)
T_right = 0.0 # Temperature at right end (°C)
t_end = 10.0 # End time (s)
N = 10 # Number of elements
dt = 0.01 # Time step size (s)
# Define basis functions (linear)
phi = lambda x: np.array([1-x, x])
# Define element mass matrix
def element_mass(h):
return (rho*Cp*h/6.0)*np.array([[2, 1], [1, 2]])
# Define element stiffness matrix
def element_stiffness(k, h):
return (k/h)*np.array([[2, -1], [-1, 2]])
# Define element force vector
def element_force(q, h):
return (q*h/2.0)*np.array([1, 1])
# Define global mass matrix
M = np.zeros((N+1, N+1))
for i in range(N):
h = L/N
M[i:i+2, i:i+2] += element_mass(h)
# Define global stiffness matrix
K = np.zeros((N+1, N+1))
for i in range(N):
h = L/N
K[i:i+2, i:i+2] += element_stiffness(k, h)
# Define global force vector
f = np.zeros(N+1)
for i in range(N):
h = L/N
f[i:i+2] += element_force(q, h)
# Apply boundary conditions
K[0, :] = 0
K[0, 0] = 1
f[0] = T_left
K[N, :] = 0
K[N, N] = 1
f[N] = T_right
# Solve system of equations numerically
T = np.zeros(N+1)
t = 0.0
while t <t_end:
T += dt*np.linalg.solve(M, f - K.dot(T))
t += dt
# Plot results
x = np.linspace(0, L, N+1)
plt.plot(x, T)
plt.xlabel('Position (m)')
plt.ylabel('Temperature (°C)')
plt.show()
The strong form of a partial differential equation (PDE) typically involves second-order derivatives, and is expressed as follows:
L(u) = f
where L is a linear differential operator, u is the unknown function, and f is a given function.
For example, the strong form of the heat conduction equation in 1D is:
∂T/∂t = k∂²T/∂x² + q
where T is the temperature, k is the thermal conductivity, q is the heat generation rate, t is the time, and x is the position.
The weak form of the same PDE is obtained by multiplying both sides of the equation by a test function v, integrating over the domain Ω, and using integration by parts:
∫Ω vL(u) dΩ = ∫Ω vf dΩ - ∫∂Ω v(αu + β∂u/∂n) ds
where α and β are constants, and ∂Ω is the boundary of the domain.
The weak form typically involves lower-order derivatives than the strong form, and only requires that the solution satisfies the equation when multiplied by the test function and integrated over the domain.
For example, the weak form of the heat conduction equation in 1D is:
∫Ω (ρCp∂T/∂tφ + k∂T/∂x∂φ) dx = ∫Ω qφ dx
where ρ is the density, Cp is the specific heat capacity, φ is the test function, and the integral is over the domain of interest.
The weak form allows for a wider range of solutions that may not necessarily be differentiable, and is often used in numerical methods such as the Finite Element Method.
In summary, the strong form expresses the differential equation in its original form, while the weak form is obtained by multiplying the differential equation by a test function and integrating over the domain. The strong form typically involves higher-order derivatives than the weak form, and provides more information about the behavior of the solution. The weak form, on the other hand, allows for a wider range of solutions and is often used for numerical solutions.