Approximate Solution
Methods
DR. TAIBA KOUSER
Approximate Solution
There are three ways to estimate approximate solution of PDEs
Least square method
Point collocation method
Galerkin’s method
To have approximate solution we have two functions
Residuals (error)
Weighting functions/trail function
Galerkin’s Function
The actual solution is approximated in the form of polynomials using
residual and weight functions.
Residual is error between exact solution and approximate solution
Linear combination of coefficient with the trial functions. In Galerkin’s
method trial/weight functions are similar.
Weighted-Residual function
Finite Element Method/Analysis
(FEM/FEA)
Fem is a numerical method used to obtain an approximate solution for a given boundary
value problem.
Finite element method can be taken as a discrete version of the Galerkin’s method.
It is widely used method in structural analysis.
•Aerospace industry
•Automotive industry
•Electrical (Transformer,Relays, etc.)
•Electronics (Integrated circuits, Appliances, etc.)
•Power industry (Nuclear, Thermal)
•Construction
Discretization Principles
Divide domain into number of discrete subdomains (every subdomain
being represented by a specific number of discrete points).
Each subdomain is called an element.
Derive algebraic equations from the Governing differential equations valid
at these discrete points
Solve the algebraic equations to obtain the dependent variables at the
discrete points.
Pre-processing
Solution
Post-processing
Example:
j j j
dT dw dT
d
(K
dT
) S
wK (K )dx wSdx
dx dx dx i i
dx dx i
l
d dT
d dT (K ) S wdx 0
(K ) S 0
dx dx
dx dx
1 l l
j
d dT dT dw dT
(K ) S wdx 0 wK (K )dx wSdx
i
dx dx dx 0 0
dx dx 0
Trial function
Tj A Bx j
T A Bx Tj Ti
T j xi Ti x j B
x xi , T Ti A xj xi
xj xi
x x j ,T Tj
Ti x j T j xi Tj Ti
T x
Ti A Bxi xj xi xj xi
xj x x xi Ni 1 at x =xi and 0 at x j
T Ti Tj
xj xi xj xi
T N iTi N jT j
Ti
T Ni N j [ N ][T ]
xj x x xi
Tj
Ni , Nj
xj xi xj xi
Ni , N j
W [ N ][W ]
are shape function Weight function is
same as trial function
in Galerkin’s method
Transpose of a scalar is always a scalar,
t t t
W [W ] [ N ]
j t j
t t j t dN dN
[W ] [ N ] ( q ) [W ] K [T ]dx [W ]t [ N ]t Sdx 0
i
i
dx dx i
j t j
j dN dN W is
t t arbitrary
[N ] ( q ) K [T ]dx [ N ] Sdx 0 matrix
i
i
dx dx i
Term 1:
j j i
t j Ni Ni Ni
[N ] ( q ) q q q
i Nj Nj Nj
i
0 qi qi
qj 0 qj
TERM 2:
EXCLUDING THE “[T]”
xj xi L
t
dN i t
j 1
dx dN i dN j j
K dx L 1 1
dN j dx dx K dx
i
i
1 L L
dx
L
dN i 1 dN j 1
dx x j xi dx xj xi
TERM 3:
j
[ N ]t S
i
By solving substitutions, we have
2
xj x SL SL
j
L 2L 2
S dx 2
i
x xi SL SL
L 2L 2
Assembly
Term 1-Term2[T]+Term3=0
Term2[T]=Term1+Term3
1 -1 Ti qi 1/ 2 These could write for
K number of nodes.
SL
L 1 1 Tj qj 1/ 2
[ K ][T ] [F ]
K will always be a square matrix.
Global matrix for 4 nodes
1 -1 0 0 T1 q1 1/ 2
-1 1+1 -1 0 T2 q2 q2 1/ 2 1/ 2
SL
0 -1 1+1 -1 T3 q3 q3 1/ 2 1/ 2
0 0 -1 1 T4 1/ 2
q4
Activated part of stiffness matrix w.r.t
nodes
Summary
Consider a problem (simply geometry, material properties)
Define loads and contraints (e.g, heat fluxes in heat problems)
Meshing
Weak form
Assembly
BCs
Solution
Summary
Element stiffness
Summary
Finite Element Analysis
Example:
By using FEA
Calculate the global stiffness matrix.
Calculate the reaction force
Solution
Nodes connectivity
Element no. Node i Node j
1 1 2
2 2 3
3 3 4
A1=200mm2 L1=100mm
E=200GPa=200*103MPa
A2=200mm2 L2=100mm
A3=100mm2 L3=200mm
𝐴𝐸 1 −1
𝐾 =
𝐿 −1 1
200 ∗ 200 ∗ 103 1 −1
𝐾1 = 200 ∗ 200 ∗ 103 1 −1
100 −1 1 𝐾2 =
100 −1 1
1 −1
=4∗ 105
1 −1 = 4 ∗ 105
−1 1 −1 1
4 −4
𝐾1 = 105 𝐾2 = 105
4 −4
−4 4 −4 4
100 ∗ 200 ∗ 103 1 −1
𝐾3 =
100 −1 1
1 −1
𝐾3 = 105
−1 1
Local stiffeness matrices are obtained. Now to calculate the global
stiffness matrix
Global stiffness matrix
4 -4 0 0
-4 8 -4 0
0 -4 5 -1
0 0 -1 1
Assembly
[ K ][ D] [F ]
4 -4 0 0 d1 f1
5 -4 8 -4 0 d2 f2
10
0 -4 5 -1 d3 f3
0 0 -1 1 d4 f4
Boundary conditions
Fixed nodes always create a reaction R. Whenever the node is fixed the displacement node is
considered to be zero. HenForce in ce the equilibrium equation is as,
Note:
Force is taken negative if
it is provided in forward
direction. It is taken
4 -4 0 0 0 R negative if its direction is
backward. For unknown
5 -4 8 -4 0 d2 20kN direction, the force
10 variable can be taken as
0 -4 5 -1 d3 0 positive and the value
will reveal its direction.
0 0 -1 1 d4 10kN
Using eliminating method
8 -4 0 d2 0.2 N
-4 5 -1 d 3 0
0 -1 1 d4 0.1N
d2 0.75mm, d3 0.1mm, d 4 0.2mm
Calculation of Reaction force
5
10 (4d1 4d 2 0d 3 0d 4) R
5
10 (4 * 0.075) R
R 30000 N
R 30kN
Finite Difference Method
Taylor Series
f’(x1)=l/h
h
Types of Difference schemes
Governing equation derivatives are expressed in terms of suitable
algebraic differences by using Taylor series expansion.
There are 3 difference schemes
Forward Difference : Approximates the derivative by considering forward
points.
Backward Difference : Uses previous points to estimate the derivative.
Central Difference : Provides a more accurate estimate by averaging
forward and backward differences.
h is spacing between two nodal points.
x
h
Forward difference scheme
From equation (1), we have forward difference scheme
Truncations error is of Order of h.
By truncating the series, we have truncation error of order “h. This term is called
leading order term.
Backward difference scheme
From equation (2)
Leading term of truncation
error
Truncations error is of Order of h.
Central difference scheme
Subtracting (1) from (2)
Leading term of truncation error
Truncation error is of order of O(h2)
Higher the truncation error, higher the accuracy will be
For second order derivative
Adding (1) and (2)
Central difference scheme for second order derivative has truncation error of O(h2)
Central difference from FD for second
order derivative
Adding (1) and (2)
f x h f (x h) 2 f ( x)
f x 2
h
Thus is central difference scheme with Truncation error of O(h2).
Note:
FDM is primarily used in
structured grids and simpler
geometries due to its ease of
implementation, though it is
less flexible for complex
geometries
Finite Difference Scheme
Calculation of “h or ∆x” in grid system
This is required matrix form of linear Algebraic equations.
Python Programming
Solution
Delta spacing effects the solution. (Left H=0.001) Right (H=1000)
2D-grid system and Calculation of “h”
FDM for 2D and 3D
2D and 3D FDM node notation
3D
2D
Grid Location Compass
Grid Location Compass notation
notation i.j,k P
i.j P i-1,j,k W
i-1,j W i+1,j,k E
i+1,j E I,j+1,k N
I,j+1 N I,j-1,k S
I,j-1 S i.j,k+1 T
i-1,j,k-1 B
Finite Difference Method
Forward difference scheme for Time
• When we have time in governing equations, such as
Similar to FD approximations, the time derivatives also can be
approximated in time.
Explicit vs Implicit scheme
Feature Explicit Method Implicit Method
Dependency Uses known values at tnt^n Uses unknown values at t^{n+1}
Stability Conditionally stable (small Δt) Unconditionally stable
Computational Cost Cheaper per step More expensive (solves a system)
Complexity Easy to implement Harder to implement
Use Case Simple or non-stiff problems Stiff or long-term simulations
Poisson equation for pressure
• For pressure, we use central difference scheme
Alternate form
Solvers
• Jacobi iteration scheme
Numerical solution of 2D Poisson Equation
(Gauss-Seidel method)
2D Unsteady Navier-Stokes equation
Number of equations and unknowns must be equal.
Here 3 unknown u,v and p and three equations. X-momemtum, Y-
Momentum and continuity equation.
Continuity equation will be converted to Poisson pressure equation.
• First, we have to find the equation for p.
Discrete form
• Take the divergence
• For incompressible flows (divergence free), there is Poisson equation.
• Solve the Poisson equation for the pressure and then use that
pressure to update the velocity field
• Algorithm to resolve the pressure-velocity coupling include
Boundary condition
• Neumann boundary condition is applied, for stationary/moving and
impermeable wall
Helmholtz-Hedge Decomposition of velocity
vector