Finite Element
Analysis
MODULE 5: Practical Applications and Software Integration
Guidelines for modeling, meshing, and convergence, Introduction to ANSYS
Workbench and MATLAB FEA scripting, Applications to aerospace
structures: wing spars, fuselage panels, Case studies and FEA validation,
Introduction to nonlinear and contact problems (overview).
Convergence Study in Finite Element Analysis
Convergence refers to the tendency of the finite element solution to approach the exact
(analytical) solution as the discretization is improved.
1. h-type Convergence/refinement
➢ Achieved by increasing the number of elements (reducing element size h)
➢ Polynomial order of shape function remains constant
➢ Commonly used in practical engineering problems
➢ Increases computational cost due to higher DOFs
➢ Example: Refining mesh near a hole in a plate
2. p-type Convergence/refinement
➢ Achieved by increasing the order of shape functions (linear → quadratic → cubic)
➢ Mesh size remains unchanged
➢ Effective for smooth stress/strain fields
➢ Requires higher-order elements
3. hp-type Convergence/refinement
➢ Combination of mesh refinement (h) and higher-order interpolation (p)
➢ Most efficient and accurate approach
➢ Used in advanced commercial solvers
Error Classification in Finite Element Method (FEM)
In the Finite Element Method, the obtained solution is always an approximation of the
exact (analytical) solution. The difference between the exact solution and the FEM
solution is called error.
1. Discretization Error (Approximation Error)
Discretization error arises because a continuous domain is replaced by a finite number
of elements, and the displacement field is approximated using shape functions.
Cause:
➢ Finite number of elements
➢ Approximate interpolation within elements
➢ Low-order shape functions
Characteristics:
➢ Dominant error in FEM,
➢ Reduces with mesh refinement,
➢ Depends on element size (h) and polynomial order (p)
Reduction Methods
➢ h-refinement (increase number of elements)
➢ p-refinement (use higher-order elements)
➢ hp-refinement (combined approach)
Example: Coarse mesh near a stress concentration region produces inaccurate stress values.
2. Interpolation Error
Interpolation error occurs because the actual displacement variation is
approximated by polynomial shape functions within an element.
Cause:
➢ Use of linear elements for curved or highly varying displacement fields
➢ Inadequate shape function order
Characteristics:
➢ Stress and strain errors are higher than displacement errors
➢ Strongly dependent on element type
Reduction Methods:
➢ Use quadratic or higher-order elements
➢ Refine mesh in regions with high gradients
Example: Using a linear bar element for a problem where displacement varies
nonlinearly.
3. Modeling Error
Modeling error arises when the physical problem is idealized or simplified
during FEM formulation.
Cause
➢ Incorrect material properties
➢ Simplified geometry
➢ Improper boundary conditions
➢ Ignoring temperature, contact, or nonlinear effects
Characteristics
➢ Cannot be reduced by mesh refinement
➢ Depends on engineering judgment
Reduction Methods
➢ Accurate material modeling
➢ Realistic boundary and loading conditions
➢ Validation with experiments or analytical solutions
Example: Assuming linear elastic behavior for a material that exhibits plasticity.
4. Numerical Error
Numerical error results from finite precision arithmetic used by computers
during calculations.
Cause
➢ Round-off errors
➢ Truncation errors
➢ Ill-conditioned stiffness matrices
Characteristics
➢ Becomes significant for very large systems
➢ Accumulates during matrix operations
Reduction Methods
➢ Use double precision
➢ Proper scaling of equations
➢ Avoid extremely distorted meshes
5. Ill-Conditioning Error
Ill-conditioning occurs when the stiffness matrix has a very large condition
number, causing small numerical errors to produce large solution errors.
Cause
➢ Large variation in material properties
➢ Poor element quality
➢ Excessive mesh refinement in limited regions
Effects
➢ Slow solver convergence
➢ Inaccurate displacements and stresses
Reduction Methods
➢ Improve mesh quality
➢ Use appropriate element types
➢ Apply proper constraints
6. Post-Processing Error
Post-processing error occurs during interpretation and visualization of results.
Cause
➢ Stress extrapolation from Gauss points to nodes
➢ Averaging of stresses
➢ Misinterpretation of contour plots
Characteristics
➢ Stresses are less accurate than displacements
➢ Peak stresses may be exaggerated
Reduction Methods
➢ Evaluate stresses at integration points
➢ Use mesh convergence checks
➢ Avoid relying on single-node peak values
Element Quality and Result Improvement
In FEM, solution accuracy strongly depends on mesh quality. Even with correct
material properties, boundary conditions, and solvers, poor-quality elements can lead
to inaccurate results, slow convergence, or complete solution failure. Therefore,
assessing element quality and applying result-improvement techniques is a critical part
of FEM analysis.
Element Quality Parameters
1. Aspect ratio: Aspect ratio is the ratio of the longest dimension of an
element to its shortest dimension. Elements with a high aspect ratio are
elongated and can lead to inaccurate stress and strain calculations,
especially in bending-dominated problems. Lower aspect ratios (close to
unity) are preferred, as they provide better numerical accuracy and faster
convergence.
2. Skewness: Skewness measures the deviation of an element’s internal angles
from the ideal values (60° for triangles and 90° for quadrilaterals). High
skewness indicates distorted elements, which reduce interpolation accuracy
and may cause numerical instability. Low skewness improves strain
distribution and solution reliability.
3. Jacobian: The Jacobian determinant represents the transformation from
natural coordinates to physical coordinates in an element. A positive Jacobian
ensures a valid, properly oriented element, while a zero or negative Jacobian
indicates an inverted or collapsed element that can cause solver failure.
Maintaining a positive and well-scaled Jacobian is critical for solution stability.
4. Element Distortion: Element distortion refers to excessive deviation of an
element’s shape from its ideal reference geometry due to skewness, stretching,
or collapsing edges. Distorted elements produce inaccurate strain-displacement
relations and may fail patch tests. Controlling distortion is essential for
achieving accurate and convergent FEM results.
Patch Test in Finite Element Method (FEM)
▪ The patch test is a basic verification test used in the Finite Element Method to
check the correctness and convergence capability of an element formulation. It
ensures that a finite element can represent constant strain and stress states without
producing spurious results.
▪ Passing the patch test is considered a minimum requirement for any finite element
to be reliable.
Purpose of the Patch Test
The patch test is performed to verify whether:
➢ The element formulation is mathematically correct
➢ The element can reproduce rigid body motion
➢ The element can represent a constant strain field
➢ The element will converge as the mesh is refined
Concept of a Patch
➢ A patch is a small group of elements (usually 2–6 elements) connected together,
forming a simple mesh.
➢ This patch is subjected to boundary conditions that theoretically produce:
▪ Zero strain (rigid body motion)
▪ Constant strain (linear displacement field)
Types of Patch Tests
➢ Rigid Body Motion Patch Test
▪ The entire patch is given a uniform displacement or rotation
▪ No deformation occurs
Expected Result: Zero strain, Zero stress
Significance:
▪ Confirms absence of spurious stiffness
▪ Verifies translational and rotational invariance
➢ Constant Strain Patch Test
▪ Linear displacement boundary conditions applied
▪ Produces uniform strain and stress
Expected Result: Constant strain across all elements, Constant stress across the patch
Significance:
▪ Confirms correct strain-displacement relationship
▪ Essential for convergence
Solver Types in Finite Element Method (FEM)
▪ In the Finite Element Method, a solver is a numerical technique used to compute
unknown nodal displacements, velocities, and accelerations by solving the
governing equilibrium equations.
▪ Based on the time integration and solution strategy, FEM solvers are broadly
classified into: (1) Implicit solvers, (2) Explicit solvers
➢ Governing Dynamic Equilibrium Equation
The general FEM equilibrium equation is:
𝑀 𝑢ሷ + 𝐶 𝑢ሶ + 𝐾 𝑢 = 𝐹 𝑡
where: 𝑀 =Mass matrix, 𝐶 =Damping matrix, 𝐾 =Stiffness matrix, 𝑢 =
Displacement vector, 𝐹 𝑡 =Time-dependent load
The way this equation is solved defines whether the solver is implicit or explicit.
Implicit Solver
An implicit solver computes the unknown response at the current time step by solving a
system of simultaneous equations. The stiffness matrix must be inverted or factorized.
➢ Mathematical Form
For static analysis:
𝐾 𝑢 = 𝐹
For nonlinear problems:
𝐾𝑇 𝑢 Δ𝑢 = Δ𝐹
Where 𝐾𝑇 =Tangent stiffness matrix
Iterative solution is required
➢ Key Characteristics of Implicit Solvers
▪ Requires global stiffness matrix
▪ Uses Newton–Raphson iteration
▪ Suitable for static and quasi-static problems
▪ Large time steps possible
▪ Computationally expensive per step
➢ Stability of Implicit Solvers
▪ Unconditionally stable
▪ Stability does not depend on time step size
▪ Accuracy still depends on time step
Advantages of Implicit Solvers
▪ Stable for large time steps
▪ Accurate for static equilibrium
▪ Well-suited for nonlinear material and contact problems
▪ Fewer total time steps required
Limitations of Implicit Solvers
▪ High memory requirement
▪ Slow for very large models
▪ Convergence difficulties in severe nonlinearity
▪ Requires matrix inversion
Applications of Implicit Solvers
▪ Static structural analysis
▪ Low-speed dynamic problems
▪ Thermal analysis
▪ Creep and plasticity
▪ Contact problems with slow motion
Implicit Solver in ANSYS
▪ Used in Static Structural, Transient Structural
▪ Newton–Raphson based
▪ Automatic time stepping
▪ Convergence monitoring
Explicit Solver
An explicit solver computes the response at the next time step directly from known
values at the current step without solving global equations.
➢ Mathematical Form
Acceleration is computed explicitly:
𝑢ሷ = ൣ𝑀]−1 𝐹 − 𝐶 𝑢ሶ − 𝐾𝑢
Displacement and velocity updated using:
1
𝑢𝑡+Δ𝑡 = 𝑢𝑡 + 𝑢ሶ 𝑡 Δ𝑡 + 𝑢ሷ 𝑡 Δ𝑡 2
2
➢ Key Characteristics of Explicit Solvers
▪ No global stiffness matrix inversion
▪ Uses diagonal (lumped) mass matrix
▪ Very small time steps required
▪ Condition-based stability
➢ Stability of Explicit Solvers
▪ Conditionally stable
▪ Time step must satisfy: Δ𝑡 ≤ Δ𝑡𝑐𝑟𝑖𝑡𝑖𝑐𝑎𝑙
▪ Critical time step depends on:
▪ Smallest element size
▪ Material wave speed
Advantages of Explicit Solvers
▪ Extremely efficient for large models
▪ Handles severe nonlinearity and contact well
▪ No convergence issues
▪ Ideal for short-duration events
Limitations of Explicit Solvers
▪ Requires very small time steps
▪ Computationally expensive for long-duration problems
▪ Less accurate for static problems
Applications of Explicit Solvers
▪ Impact and crash analysis
▪ Drop tests
▪ Metal forming
▪ Explosions
▪ High-speed contact problems
Explicit Solver in ANSYS
▪ Used in ANSYS Explicit Dynamics
▪ Automatic time step control
▪ Efficient for contact-dominated problems
Comparison: Implicit vs Explicit Solvers
Feature Implicit Solver Explicit Solver
Stability Unconditionally stable Conditionally stable
Time step Large Very small
Matrix inversion Required Not required
Convergence Iterative No iterations
Contact handling Moderate Excellent
Best suited for Static / slow dynamics Impact / crash
Solver Selection Guidelines
Feature Implicit Solver Explicit Solver
Stability Unconditionally stable Conditionally stable
Time step Large Very small
Matrix inversion Required Not required
Convergence Iterative No iterations
Contact handling Moderate Excellent
Best suited for Static / slow dynamics Impact / crash
Introduction to Nonlinear Problems in Finite Element Method (FEM)
➢ Linear vs Nonlinear FEM Problems
In linear FEM, the relationship between forces and displacements is linear, and the principle of
superposition is valid. In contrast, nonlinear FEM problems exhibit a nonlinear relationship
between applied loads and structural response, and superposition is no longer applicable.
A problem is said to be nonlinear if any one of the following is nonlinear:
▪ Material behavior
▪ Geometry (deformation)
▪ Boundary conditions
➢ Characteristics of Nonlinear FEM Problems
▪ Load–displacement relationship is nonlinear
▪ Stiffness matrix depends on displacement
▪ Solution must be obtained incrementally
▪ Requires iterative solvers
▪ Computationally expensive
▪ Path-dependent response
➢ Governing Equation for Nonlinear FEM
For linear problems: 𝐾 𝑢 = 𝐹
For nonlinear problems: 𝐾 𝑢 𝑢 = 𝐹 𝑢
Where:𝐾 𝑢 =displacement-dependent stiffness matrix, 𝑢= unknown displacement vector, 𝐹 𝑢 =
load vector may also depend on deformation
Types of Nonlinearity in FEM
1. Material Nonlinearity
Material nonlinearity occurs when stress is not proportional to strain.
Causes: Plastic deformation, Creep, Viscoelasticity, Hyperelasticity
Examples:
▪ Metal yielding beyond elastic limit
▪ Rubber components
▪ Polymer materials
FEM Treatment:
▪ Incremental stress–strain relations
▪ Use of constitutive models (von Mises, Mohr–Coulomb, etc.)
2. Geometric Nonlinearity
Geometric nonlinearity arises when deformations are large, causing changes in
geometry during loading.
Causes: Large displacement, Large rotation, Large strain
Effects: Change in stiffness, Change in load direction
Examples: Buckling of columns, Cable structures, Thin shells under large deflection
FEM Treatment:
•Updated or total Lagrangian formulation
•Nonlinear strain–displacement relations
3. Boundary Condition Nonlinearity
Boundary nonlinearity occurs when boundary conditions change with
deformation.
Causes:
▪ Contact and separation
▪ Friction
▪ Impact problems
Examples:
▪ Gear tooth contact
▪ Bolted joints
▪ Sheet metal forming
FEM Treatment:
▪ Contact algorithms
▪ Penalty and Lagrange multiplier methods
Comparison: Linear vs Nonlinear FEM
Feature Linear FEM Nonlinear FEM
Stiffness matrix Constant Variable
Superposition Valid Invalid
Load application Single step Incremental
Solver Direct Iterative
Computation cost Low High
Applications of Nonlinear FEM
➢ Plastic deformation of metals
➢ Buckling and post-buckling analysis
➢ Contact analysis
➢ Rubber and elastomer components
➢ Crash and impact analysis
➢ Aerospace structural collapse analysis
Contact Problems in Finite Element Method
A contact problem in FEM arises when two or more bodies interact by touching, separating, or
sliding against each other during deformation. Unlike fixed boundary conditions, contact
conditions change with time and deformation, making the problem highly nonlinear.
Contact problems are a major source of boundary nonlinearity in FEM.
Contact problems are nonlinear because:
▪ Contact status changes continuously (open closed)
▪ Contact forces depend on deformation
▪ Friction introduces nonlinear tangential forces
▪ Stiffness matrix changes during solution
As a result:
▪ Superposition principle is invalid
▪ Incremental-iterative solution is required
Physical Characteristics of Contact
➢ Contact Conditions
At the interface between two bodies: No penetration (normal constraint), Possible separation,
Possible sliding with friction
➢ Contact Constraints
▪ Normal direction: prevents penetration
▪ Tangential direction: governs frictional behavior
Types of Contact in FEM
➢ Based on Geometry: Node-to-Node contact, Node-to-Surface contact, Surface-to-Surface
contact (most accurate)
➢ Based on Motion: Frictionless contact, Frictional contact, Rough contact (no sliding)
➢ Based on Interaction: Bonded contact (no separation or sliding), No-separation contact,
Standard contact (separation allowed)
Contact Kinematics
➢ Normal Gap
𝑔𝑛 = distance between contacting surfaces
▪ 𝑔𝑛 > 0 :Open (no contact)
▪ 𝑔𝑛 = 0 :Contact
▪ 𝑔𝑛 < 0 :Penetration (not allowed)
➢ Tangential Slip
▪ Occurs when friction is present
▪ Governed by Coulomb’s friction law
➢ Friction in Contact Problems
Coulomb Friction Law
𝐹𝑡 ≤ 𝜇𝐹𝑛
where:𝐹𝑡 =tangential (friction) force, 𝐹𝑛 =normal contact force, 𝜇= coefficient of friction
➢ Friction States
▪ Stick condition: No relative sliding
▪ Slip condition: Relative motion occurs
Problem 1: A compound bar of total length 2 m is fixed at the left end and
subjected to an axial tensile load of 60 kN at the right end.
The bar consists of two segments with different material and cross-sectional
properties. Using the Finite Element Method with two-noded bar elements,
write a MATLAB script to determine: (i) nodal displacements, (ii) element
strains, and (iii) element stresses.
Code in MATLAB script (2-noded elements):
clc; clear;
%% INPUT DATA:
E1 = 200e9; E2 = 100e9; % Young's modulus of element 1 and 2 (Pa)
A1 = 400e-4; A2 = 200e-4; % Area of element 1 and 2 (m^2)
L1 = 1; L2 = 1; % Length of element 1 and 2 (m)
P = 60e3; % Point load at node 3 (N)
%% ELEMENT STIFFNESS MATRICES:
k1 = (E1*A1/L1) * [ 1 -1; -1 1];
k2 = (E2*A2/L2) * [ 1 -1; -1 1];
%% GLOBAL STIFFNESS MATRIX:
K = zeros(3,3);
% Assembly
K(1:2,1:2) = K(1:2,1:2) + k1;
K(2:3,2:3) = K(2:3,2:3) + k2;
%% GLOBAL FORCE VECTOR
F = zeros(3,1);
F(3,1) = P; % Point load at node 3
%% BOUNDARY CONDITION
% Node 1 is fixed → u1 = 0
K_reduced = K(2:3,2:3);
F_reduced = F(2:3);
%% SOLVE FOR DISPLACEMENTS
U_reduced = K_reduced \ F_reduced;
U = zeros(3,1);
U(2:3) = U_reduced;
%% ELEMENT STRAINS
strain1 = (U(2) - U(1)) / L1;
strain2 = (U(3) - U(2)) / L2;
%% ELEMENT STRESS
stress1 = E1 * (U(2) - U(1)) / L1;
stress2 = E2 * (U(3) - U(2)) / L2;
%% DISPLAY RESULTS
disp('Nodal Displacements (m):’); disp(U);
disp('Element Strains:’); disp(strain1); disp(strain2);
disp('Element Stresses (Pa):’); disp(stress1); disp(stress2);
Problem: AA compound bar is rigidly fixed to a wall at one end and subjected to
axial load P = 5 kN as illustrated in Figure. The bar consists of two segments with
varying cross-sectional areas (A), lengths (L), and Young’s moduli (E). Using the
Finite Element Method with two-noded bar elements, write a MATLAB
script to determine: (i) nodal displacements, (ii) element strains, and (iii)
element stresses.
Code in MATLAB script (2-noded elements):
clc; clear;
%% INPUT DATA:
E1 = 70e9; E2 = 200e9; % Young's modulus of element 1 and 2 (Pa)
A1 = 900e-4; A2 = 2000e-4; % Area of element 1 and 2 (m^2)
L1 = 500e-3; L2 = 250e-3; % Length of element 1 and 2 (m)
P = -5e3; % Compressive Point load at node 2 (N)
%% ELEMENT STIFFNESS MATRICES:
k1 = (E1*A1/L1) * [ 1 -1; -1 1];
k2 = (E2*A2/L2) * [ 1 -1; -1 1];
%% GLOBAL STIFFNESS MATRIX:
K = zeros(3,3);
% Assembly
K(1:2,1:2) = K(1:2,1:2) + k1;
K(2:3,2:3) = K(2:3,2:3) + k2;
%% GLOBAL FORCE VECTOR
F = zeros(3,1);
F(2,1) = P; % Point load at node 2
%% BOUNDARY CONDITION
% Node 1 is fixed → u1 = 0
K_reduced = K(2:3,2:3);
F_reduced = F(2:3);
%% SOLVE FOR DISPLACEMENTS
U_reduced = K_reduced \ F_reduced;
U = zeros(3,1);
U(2:3) = U_reduced;
%% ELEMENT STRAINS
strain1 = (U(2) - U(1)) / L1;
strain2 = (U(3) - U(2)) / L2;
%% ELEMENT STRESS
stress1 = E1 * (U(2) - U(1)) / L1;
stress2 = E2 * (U(3) - U(2)) / L2;
%% DISPLAY RESULTS
disp('Nodal Displacements (m):’); disp(U);
disp('Element Strains:’); disp(strain1); disp(strain2);
disp('Element Stresses (Pa):’); disp(stress1); disp(stress2);
Problem 2: A compound bar of total length 2 m is fixed at the left end and
subjected to an axial tensile load of 60 kN at the right end.
The bar consists of two segments with different material and cross-sectional
properties. Using the Finite Element Method with three-noded bar elements,
write a MATLAB script to determine: (i) nodal displacements, (ii) element
strains, and (iii) element stresses.
Code in MATLAB script (3-noded elements):
clc; clear;
%% INPUT DATA:
E1 = 200e9; E2 = 100e9; % Young's modulus of element 1 and 2 (Pa)
A1 = 400e-4; A2 = 200e-4; % Area of element 1 and 2 (m^2)
L1 = 1; L2 = 1; % Length of element 1 and 2 (m)
P = 60e3; % Point load at node 3 (N)
%% ELEMENT STIFFNESS MATRICES:
k1 = ((E1*A1)/(3*L1)) * [ 7 -8 1; -8 16 -8; 1 -8 7];
k2 = ((E2*A2/(3*L2)) * [ 7 -8 1; -8 16 -8; 1 -8 7];
%% GLOBAL STIFFNESS MATRIX:
K = zeros(5,5);
% Assembly
K(1:3,1:3) = K(1:3,1:3) + k1; % Element 1 → Nodes (1,2,3)
K(3:5,3:5) = K(3:5,3:5) + k2; % Element 2 → Nodes (3,4,5)
%% GLOBAL FORCE VECTOR
F = zeros(5,1);
F(5,1) = P; % Point load at node 5
%% BOUNDARY CONDITION
% Node 1 is fixed → u1 = 0
K_reduced = K(2:5,2:5);
F_reduced = F(2:5);
%% SOLVE FOR DISPLACEMENTS
U_reduced = K_reduced \ F_reduced;
U = zeros(5,1);
U(2:5) = U_reduced;
%% ELEMENT STRAINS
strain1 = (U(3) - U(1)) / L1;
strain2 = (U(5) - U(3)) / L2;
%% ELEMENT STRESS
stress1 = E1 * (U(3) - U(1)) / L1;
stress2 = E2 * (U(5) - U(3)) / L2;
%% DISPLAY RESULTS
disp('Nodal Displacements (m):’); disp(U);
disp('Element Strains:’); disp(strain1); disp(strain2);
disp('Element Stresses (Pa):’); disp(stress1); disp(stress2);
Problem 3: A straight prismatic bar of total length (L) = 2 m is subjected to an
axial uniformly distributed tensile load of intensity (q) = 60 kN/m along its
length. The bar has a constant cross-sectional area (A) of 400 mm² and is made of
a linear elastic material with Young’s modulus (E) 200 GPa. The left end of the
bar is fixed, while the right end is free. Using the Finite Element Method, the
bar is discretized into two equal two-noded bar elements. Write a MATLAB
script to determine: (i) nodal displacements, (ii) element strains, and (iii) element
stresses.
Code in MATLAB script (2-noded elements):
clc; clear;
%% INPUT DATA
E = 200e9; % Young's modulus (Pa)
A = 400e-6; % Cross-sectional area (m^2)
L = 2; % Total length of bar (m)
q = 60e3; % UDL (N/m)
nElem = 2; % Number of elements
nNode = nElem + 1; % Number of nodes
Le = L / nElem; % Length of each element
%% INITIALIZE GLOBAL MATRICES
K = zeros(nNode,nNode);
F = zeros(nNode,1);
%% ELEMENT STIFFNESS MATRIX
ke = (E*A/Le) * [ 1 -1; -1 1];
%% ELEMENT FORCE VECTOR (UDL)
fe = q*Le/2 * [1; 1]; % Equivalent nodal forces for 2-noded bar under UDL
%% ASSEMBLY PROCESS
for e = 1:nElem
node1 = e; node2 = e + 1;
K(node1:node2, node1:node2) = ...
K(node1:node2, node1:node2) + ke;
F(node1:node2) = F(node1:node2) + fe;
end
%% APPLY BOUNDARY CONDITION
K_reduced = K(2:end,2:end); % Node 1 is fixed (u1 = 0)
F_reduced = F(2:end);
%% SOLVE FOR NODAL DISPLACEMENTS
U_reduced = K_reduced \ F_reduced;
U = zeros(nNode,1);
U(2:end) = U_reduced;
%% STRAIN AND STRESS CALCULATION
strain = zeros(nElem,1);
stress = zeros(nElem,1);
for e = 1:nElem
strain(e) = (U(e+1) - U(e)) / Le;
stress(e) = E * strain(e);
end
%% DISPLAY RESULTS
disp('Nodal Displacements (m):’); disp(U);
disp('Element Strains:’); disp(strain);
disp('Element Stresses (Pa):’); disp(stress);
Problem 3: A straight prismatic bar of total length (L) = 2 m is subjected to an
axial uniformly distributed tensile load of intensity (q) = 60 kN/m along its
length. The bar has a constant cross-sectional area (A) of 400 mm² and is made of
a linear elastic material with Young’s modulus (E) 200 GPa. The left end of the
bar is fixed, while the right end is free. Using the Finite Element Method, the
bar is discretized into two equal three-noded bar elements. Write a MATLAB
script to determine: (i) nodal displacements, (ii) element strains, and (iii) element
stresses.
Code in MATLAB script (3-noded elements):
clc; clear;
%% INPUT DATA
E = 200e9; % Young's modulus (Pa)
A = 400e-4; % Cross-sectional area (m^2)
L = 2; % Length of each element (m)
q = 60e3; % UDL (N/m)
%% NUMBER OF NODES
nNodes = 5; % Two 3-noded elements → 5 nodes total
%% INITIALIZE GLOBAL MATRICES
K = zeros(nNodes);
F = zeros(nNodes,1);
%% ELEMENT STIFFNESS MATRIX (Quadratic Bar)
ke = (E*A/(3*L)) * [ 7 -8 1; -8 16 -8; 1 -8 7];
%% ELEMENT FORCE VECTOR (UDL)
fe = (q*L/6) * [1; 4; 1];
%% ASSEMBLY
% Element 1 → Nodes (1,2,3)
K(1:3,1:3) = K(1:3,1:3) + ke;
F(1:3) = F(1:3) + fe;
K(3:5,3:5) = K(3:5,3:5) + ke; % Element 2 → Nodes (3,4,5)
F(3:5) = F(3:5) + fe;
%% BOUNDARY CONDITION
% Node 1 fixed → u1 = 0
K_red = K(2:end,2:end);
F_red = F(2:end);
%% SOLVE FOR DISPLACEMENTS
U_red = K_red \ F_red;
U = zeros(nNodes,1);
U(2:end) = U_red;
%% STRAIN CALCULATION
strain1 = (U(3) - U(1)) / (2*L);
strain2 = (U(5) - U(3)) / (2*L);
%% STRESS CALCULATION
stress1 = E * strain1;
stress2 = E * strain2;
%% DISPLAY RESULTS
disp('Nodal Displacements (m):');
disp(U);
disp('Element Strains:’);
disp(strain1); disp(strain2);
disp('Element Stresses (Pa):’);
disp(stress1); disp(stress2);