KIT Workshop

Dr. Pritam Chakraborty

Asst. Prof., AE, IIT Kanpur

16th January, 2020

Outreach Audi., IIT Kanpur
 Outline

▪ Nonlinearity
Definition, Examples.

▪ Newton Raphson
▪ Nonlinear 1D Bar Problem
ODEs, weak Form, solution strategy, MATLAB code

▪ Incremental Method
▪ Postprocessing
▪ Other Direct Solvers
▪ Problems to Solve
 References

▪ Nonlinear Finite Elements for Continua and Structures

By Ted Belytschko, Wing Kam Liu, Brian Moran, Khalil Elkhodary

▪ Finite element procedures

By Klaus-Jürgen Bathe

▪ Non‐Linear Finite Element Analysis of Solids and Structures, Vol 1

By M. A. Crisfield

▪ Nonlinear Finite Element Methods

By Peter Wriggers
 Nonlinearity

A physical system (a body or assembly, etc.) is defined as nonlinear

if the output(s) or response(s), y, of the system due to input(s), x, do
not possess a linear relationship, or,

𝑦 ≠ 𝑚𝑥 + 𝑐

In structural mechanics, the outputs and inputs are forces, stresses,

displacements, strains, etc.

Nonlinearity implies that force do not follow a linear relationship with

displacement, or vice-versa.
 Types of Nonlinearity

▪ Geometric
Involves large displacement and rotation, and/or deformation of structural
elements

Examples: cables, membranes, frames, tyres, metal forming

Buckling
Hyperelastic seal Cable
 Types of Nonlinearity

▪ Physical or Material
Nonlinear response function between stress and strain
Examples: metal plasticity, viscoelasticity, soil, damage

Stress-strain Damage of Porous

Material
 Types of Nonlinearity

▪ Boundary Conditions
Deformation dependent boundary condition, contact

Contact Rolling

Geometry, material and

boundary nonlinearity
 Geometric Nonlinearity

Example: Large displacement of a rigid beam [Wriggers, 2008]

Initial Configuration From Moment Balance 𝐹𝑙 cos ∅ = 𝑐 ∅

Hence, spring rotation, f, output, is a nonlinear

function of applied force, F, input.

Equilibrated Configuration
 Geometric Nonlinearity

Example: Large displacements of elastic springs [Wriggers, 2008]

Initial Configuration Force Balance

𝐹
Equilibrated Configuration 2𝑐𝑙
𝑤
𝐹 = 2𝑇𝑠𝑖𝑛 𝜃 = 2𝑇
𝑙+𝑓

For a linear spring 𝑇 = 𝑐𝑓

Kinematics Hence, 𝑤 1 𝐹
1− =
𝑙 𝑤 2 2𝑐𝑙
Deflection 1+
𝑙

𝑤 2 Or, deflection of spring mid-point, w, output, is a

Elongation 𝑓=𝑙 1+ −1
𝑙 nonlinear function of applied force, F, input.
 Geometric Nonlinearity

Example: Snap through [Wriggers, 2008]

For a linear spring 𝑇 = 𝑐𝑓
Initial Configuration Equilibrated Configuration
𝑤−ℎ 𝐿0 𝐹
Truss members
Hence, 𝑙
1− =
2𝑐𝑙
2
𝑤−ℎ
𝑙 1+
𝑙

𝐿0
= 1.25
𝑙
Kinematics Force Balance
𝐹
2𝑐𝑙

ℎ−𝑤
𝐹 = 2𝑇𝑠𝑖𝑛 𝛼 − 𝜑 = 2𝑇
𝐿
Original length - L0 ; Current Length - L
Deflection of A - w
2 2 Loss of equilibrium at D and restored at E
ℎ − 𝑤 2 + 𝑙 2 = 𝐿2 ℎ−𝑤 ℎ Snap through –equilibrium at D or E for a F
𝑓 =𝑙 1+ − 1+
ℎ 2 + 𝑙 2 = 𝐿0 2 𝑙 𝑙 and different w.
 Material Nonlinearity

Example: Response of elasto-plastic bar assembly [Wriggers, 2008]

Elasto-Perfectly Plastic When bars 1 and 2 are elastic
𝐹 2𝐹 𝐹
𝜎1 = 2𝐸𝜀 𝜎2 = 𝐸𝜀 𝜀= 𝜎1 = 𝜎2 =
3𝐸𝐴 3𝐴 3𝐴
3
ሖ
Σ = 𝐸𝜀 = 𝐸𝜀
2

Onset of plasticity in bar 2

𝜎𝑌 3
𝜎1 = 2𝐸𝜀 ∗ 𝜎2 = 𝜎𝑌 𝜀∗ = Σ∗ = 𝜎𝑌
𝐸 2
Bars 1 and 2 experience same strain, e, as
that of assembly Bar 1 elastic and bar 2 plastic
Force Equilibrium: F = F1 + F2 = s1 A + s2 A 𝐹 = 𝜎𝑌 𝐴 + 2𝐸𝜀𝐴 𝜀Ƹ = 𝜀 − 𝜀 ∗ 𝐹 = 3𝜎𝑌 𝐴 + 2𝐸𝐴𝜀Ƹ
Σ෠ = Σ − Σ ∗ = 𝐸 𝜀Ƹ
Stress of the assembly: S = F/2A
 Material Nonlinearity

Onset of plasticity in bar 1

𝜎1 = 3𝜎𝑌 𝜎2 = 𝜎𝑌

3𝜎𝑌 3 ∗
𝜀 ∗∗ = = 𝜀
2𝐸 2
𝐹 = 𝜎𝑌 𝐴 + 3𝜎𝑌 𝐴
Σ ∗∗ = 2𝜎𝑌

Bar 1 and 2 plastic

Σ = 2𝜎𝑌
 Newton Raphson

In the examples, nonlinear algebraic equations describe the relationship

between forces (stresses) and displacements (deflections, strains)

1st example: 𝐹𝑙 = ∅
𝑐 cos∅
𝐹𝑙
Evaluating for a ∅ is trivial. However the reverse requires solving of roots.
𝑐

Newton-Raphson is a widely used method for finding out roots.

Iteratively solves by linearizing a nonlinear equation.

 Algorithm

Problem: f(u) is some function of u. Find u0 such that f(u0) = 0.

Solution: Let us guess a u = u1. No that lucky, so, f(u1) ≠ 0
But at u0 = u1 + Δu1, which we don’t know, f(u0) = 0 or f(u1 + Δu1) = 0
Taylor series expansion of f(u) about u1:
𝑑𝑓 1 𝑑2 𝑓
𝑓 𝑢1 + ∆𝑢1 = 𝑓 𝑢1 + ቚ ∆𝑢1 + ቚ ∆𝑢12 + ⋯ = 0
𝑑𝑢 𝑢1 2 𝑑𝑢2 𝑢1

Truncate the series after linear term and equate to zero (Linearization)
𝑑𝑓 𝑓 𝑢1
Step 1: 𝑓𝐿 𝑢1 + ∆෦
𝑢1 = 𝑓 𝑢1 + 𝑑𝑢ቚ
𝑢1
∆෦
𝑢1 = 0 Update: ∆෦
𝑢1 = −
𝑑𝑓 𝑢2 = 𝑢1 + ∆෦
𝑢1
ฬ
𝑑𝑢 𝑢1

However, u2 need not be u0 or f(u2) ≠ 0

Repeat step 1 about u2 to obtain an updated value
 Flow chart

Initial Guess Is
Y

Update
i = Iteration number
 Example

𝐹𝑙 𝑢 𝑢0 𝑓 𝑢 = 2 cos 𝑢 − 𝑢
= Find u0 such that = 2 Or
𝑐 cos(𝑢) cos(𝑢0 ) And find u0 such that f(u0) = 0

Initial guess of u = 0 and df/du = -2sin(u) - 1

𝑓 𝑢1

𝑓𝐿 𝑢2 = 0 𝑓 𝑢3

𝑓𝐿 𝑢4 = 0
𝑓𝐿 𝑢5 = 0
𝑓𝐿 𝑢3 = 0 𝑓 𝑢4
𝑓 𝑢2

u 0 2 0.99514 1.030107 1.029867

f(u) 2 -2.83229 0.093631 -0.00065 -2.98E-08
MATLAB Code
 Boundary Value Problems

Mechanical behaviour of continuous bodies or their assembly can be suitably

represented by boundary value or initio-boundary value problems
Example: Elastic Cantilever Beam Boundary Conditions
𝜕𝜎𝑖𝑗 • At x = y = z = 0
Governing =0 u=v=w=0
z 𝜕𝑥𝑗
• At x = 400 y = -60 z = 40
Constitutive 𝜎ሶ𝑖𝑗 = 𝜎𝑖𝑗 (𝐷𝑘𝑙 )
Fz = -2500
𝜕𝑣𝑖
y Kinematic 𝐿𝑖𝑗 = 𝐷𝑖𝑗 + 𝑊𝑖𝑗 = • On all the faces (exclude
𝜕𝑥𝑗 the fixed face)
x
tx = ty = tz = 0
Problem can be nonlinear
• geometrically if cantilever beam undergoes large bending
• physically (material) if the stress or stress rate is a nonlinear function of displacement
or velocity gradients
Hence, solve nonlinear PDE(s) to model deformation behaviour of continuous bodies
 Nonlinear Differential Equations

Let L is a differential operator and operates on some scalar or vector or tensor valued
𝜕
function, u. (e.g. 𝐿 = )
𝜕𝑥
The operator is said to be linear if for w = u + v, Lw = Lu + Lv

If the above doesn’t hold then the ordinary or partial differential equation is nonlinear
Example: Linear Homogeneous ODE
𝑑2𝑢 𝑑2𝑣 𝑑 2𝑤
+ 𝑢 = 0 for some u(x), and + 𝑣 = 0 for some v(x). Then if w(x) = u(x) + v(x), then +𝑤 =0
𝑑𝑥 2 𝑑𝑥 2 𝑑𝑥 2

Example: Nonlinear Homogeneous ODE

𝑑2𝑢 2 𝑑2𝑣
2
+ 𝑢 = 0 for some u(x), and + 𝑣 2 = 0 for some v(x). Then if w(x) = u(x) + v(x), then
𝑑𝑥 𝑑𝑥 2
𝑑2𝑤 2
𝑑2𝑢 2
𝑑2𝑣 2 + 2𝑢𝑣 ≠ 0
+ 𝑤 = + 𝑢 + + 𝑣
𝑑𝑥 2 𝑑𝑥 2 𝑑𝑥 2
 FEM for Nonlinear Differential Equations
▪ A numerical technique for solving ordinary (ODEs) and partial differential equations (PDEs)

▪ Based on decomposition of domain into finite regions or elements

▪ Interpolative approximation of the primary variable(s) in every element

▪ Weak form of ODE or PDE in every element gives algebraic equations

▪ Linear: Stiffness matrix and force vector (RHS vector) derivable from weak form
▪ Nonlinear: Linearization to obtain the tangent stiffness (Jacobian) matrix and RHS vector (Residual)

▪ Assembly of tangent stiffness matrix and RHS vector for all elements

▪ Apply boundary conditions

▪ Linear: Solve to obtain unknown nodal values

▪ Nonlinear: Solve to obtain update on initial guess of nodal values
Calculate RHS vector and check for convergence (Newton Raphson)
 Example: Nonlinear bar in 1D

A0(X) A tapered bar fixed on left face and displaced

A(x) by u0 on right face.

u0 The bar is incompressible elastic and

𝑒1Ƹ
undergoes large deformation - nonlinear

L
L0
Assumption: Stress and strain components,
other than sxx and exx , are zero
x(X,t)

X: Position vector of a
particle, A, at t = 0.
Deformed or Current x: Position vector of A at t
Initial or Reference
Configuration (t)
Configuration (t = 0)
 Kinematics
Motion: Maps any particle from the initial position to its current position or x = x( X, t)
For a continuous body (without discontinuities such as cracks), the mapping is unique
Thus, inverse exists, i.e. X = x-1( x, t)
Displacement of particle A: u = x – X
In terms of initial position, u = 𝑢ො 𝑋, 𝑡 = 𝑥 𝑋, 𝑡 − 𝑋
current position, u = 𝑢෤ 𝑥, 𝑡 = 𝑥 − 𝑥 −1 (𝑥, 𝑡)

Deformation Gradient (F) of A and its small neighbourhood

t=0 Taylor series till linear term
t
𝑑𝑥
∆𝑥 = 𝑥 𝑋 + ∆𝑋, 𝑡 − 𝑥 𝑋, 𝑡 ∆𝑥 ≈ 𝑥 𝑋, 𝑡 + (𝑋, 𝑡)∆𝑋 − 𝑥 𝑋, 𝑡
𝑑𝑋
A A’ A A’
(X) (X+DX) (x) (x+Dx) 𝑑𝑥 Deformation 𝑑𝑥 𝑑𝑢ො
∆𝑥 = 𝑋, 𝑡 ∆𝑋 = 𝐹∆𝑋 𝐹= =1+
𝑑𝑋 Gradient 𝑑𝑋 𝑑𝑋
Initial Current
 Kinematics

Adopting the Green-Lagrange strain definition to

account for large stretching of material filament

This strain is defined w.r.t. to the infinitesimal length of filament in initial position or DX.
and strain is defined in the initial configuration of the body (t = 0).

In terms of u:

Under small strain assumption, x ~ X, giving

Thus the Green-Lagrange strain introduces a nonlinear term to account for large
deformation
 Force Balance

Deformed body
under force P P P P(x)
equilibrium Cutting Plane
x
In current configuration: P(x) = P = constant Or P(x(X,t)) = P’(X,t) = constant
Condition of force equilibrium
in initial configuration

In terms of stress
sN – Nominal stress component sxx; s – True stress component sxx

From incompressibility: A0 DX = A Dx or A0 = AF
 Constitutive Equation

Define a linear nominal stress vs Green Lagrange strain relation

Both the stress and strain are defined in the initial configuration
 Nonlinear ODE

Governing Equation

Strain - Displacement

Stress - Strain

Boundary Conditions: u(x = 0) = 0 and u(x = L0) = u0

 Lagrangian and Eulerian Descriptions

Lagrangian Description
Balance Laws, Kinematics and Constitutive equations are described in the initial
configuration (X)

Eulerian Description
Balance Laws, Kinematics and Constitutive equations are described in the current
configuration (x)

Equilibrium

The constitutive and kinematics are described in the rate form

Hypoelastic
 Lagrangian FEM

In structural/solid mechanics, the interest is the deforming body and its state of stress, etc.
Hence, the solving methodology must track the motion of the domain or the initial mesh
(discretized domain)

Two widely used approaches:

Total Lagrangian :

▪The Lagrangian description of mass and momentum balance, kinematics and constitutive
equations are solved
▪Interpolation and weak form are defined over the initial configuration
▪Mesh remains undistorted
 Lagrangian FEM

Updated Lagrangian :

▪The Eulerian description of mass and momentum balance, kinematics and constitutive
equations are solved
▪The deformation gradient is retained to connect the configurations
▪Interpolation and weak form are defined over the current configuration
▪Mesh distorts with deformation

▪Used in FEM softwares

 Demonstration of solving the 1-D problem using the
Total Lagrangian Method and related derivations
▪Derive the residual vector and tangent stiffness matrix for 1 element
▪Assembly for n elements
▪Application of BCs
▪MATLAB CODE

▪Use of Gauss Quadrature

▪General representation of residual vector and tangent stiffness matrix
▪MATLAB CODE
Assembly – Pseudo Code
Assembly – Pseudo Code
 Incremental Method

In structural problems, the loads and/or displacements are applied over a finite interval
of time.
Applied
When solving, the interest is only at
PF dF certain times or state
P d In previous examples, stress, strain,
OR etc. were obtained at (PF, dF)

t t For nonlinear problems, achieving

the final state from initial state can
PF, dF lead to convergence difficulties
P
Response One reason – Due to an initial guess
far from the actual solution
d
 Incremental Method

Divide the applied loads and/or displacements into finite number of intervals

PF Time points - t1, t2, ...

P Time Increment – Dtn = tn+1- tn , n – discrete points

For static problems with rate independent

constitutive models, time points and increments
are markers
t1 t2 t3 t4 t5 tf
t The applied load and/or displacement increments
determine convergence of N-R iterations
 Example: Incremental Method

Example: Using the tapered bar problem subjected to load at the right end

The load is assumed to vary linearly with time followed by

PF
equal divisions of the total
P
The final state is reached in the following manner:
P2
P1 Knowing the initial state, first solve for applied load of P1
using NR (initial guess t = 0)
t1 t2 t3 t4 t5 tf
Knowing the state at t1, solve for applied load of P2 using
t
NR (initial guess t = t1)

So on and so forth till PF

 Demonstration of solving the 1-D problem using the
Incremental method
▪MATLAB code
Flow chart
 Flow chart

IN

OUT
 Postprocessing
▪ Nodal values of only the primary variables are accurate without any further
manipulation

▪ Stress, strain, etc. are accurate only at the Gauss points

▪ Nodal stresses, etc. are extrapolation or obtained through other specialized

techniques (superconvergent patch recovery [Zienkiewicz, 1992])

▪ Nodal forces (reactions) can be obtained as follows for the 2 element 2 Dirichlet
BC problem - once the NR iterations have converged.
 Some Other Direct Solvers

▪ NR method has quadratic convergence

▪ However, requires
▪The exact tangent stiffness matrix to be calculated in each iteration
▪Assembly of the local matrices
▪Inversion (Factorization) of the tangent stiffness matrix in every iteration –
n3 operation

▪ While still the most effective and widely used, other modified versions of
Newton’s update can help
▪Reducing number of factorization
▪Calculation of exact tangent stiffness
 Modified Newton Raphson

▪The tangent stiffness matrix is derived only in the first iteration of the iterative
update of a time point
▪ The choice of time points depends on the nonlinearity (iterations to converge)
 Quasi Newton Method

▪The iterations in a time point starts with an approximate stiffness matrix

▪ Updated iteratively by the enforcing the secant condition

Secant
condition

▪The Broyden-Fletcher-Goldfarb-Shanno (BFGS) method is the most popular for

the secant update of the matrix – rank two update
BFGS Update
 Problem: 1d-linear bar with nonlinear body force

Boundary Conditions: u(x = 0) = 0 and u(x = L0) = u0

▪ Derive the residual vector and tangent stiffness matrix for an element
using linear Lagrange shape functions.
▪ For the bar with 2 elements, show the equation(s) to solve.
 Problem 2: 1d nonlinear bar with C2 elements

Governing Equation

Strain - Displacement

Stress - Strain
Boundary Conditions: u(x = 0) = 0 and u(x = L0) = u0

▪ Use quadratic shape functions to derive the residual vector and

tangent stiffness matrix.
▪ For a single element obtain the nonlinear equation(s) to solve.
▪ Find the suitable n-point rule if quadratic shape functions are used in
the master element.