Advanced Computational Mechanics - Appunti

Politecnico di Milano

Facoltà di Ingegneria

Corso di laurea in Ingegneria Civile

Questions

Advanced computational mechanics

Student: Lorenzo Sostegni

Matricola: 996088

Academic Year 2022-2023

Contents

1 Potential problems 1

1.1 Question 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1.1 Taut membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1.2 Anti-plane shear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Question 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2.1 Heat conduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2.2 Torsion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.3 Question 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.4 Question 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2 Nonlinear problems 9

2.1 Question 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.1.1 Newton-Raphson method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.1.2 Modified Newton-Raphson method . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2 Question 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3 Question 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.4 Question 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.5 Question 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3 Elastoplasticity 17

3.1 Question 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2 Question 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.2.1 Return mapping algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.3 Question 12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.3.1 von Mises associative model with linear isotropic hardening . . . . . . . . . . . . . 20

3.4 Question 13 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

4 Integral equations 23

4.1 Question 14 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

4.2 Question 15 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

4.2.1 Derivation of the fundamental solution . . . . . . . . . . . . . . . . . . . . . . . . . 26

4.3 Question 16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4.3.1 BIE for potential problem in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

4.4 Question 17 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

Advanced computational mechanics I

1

Potential problems

1.1 Question 1

“Potential problems” governed by Laplace/Poisson equation: taut membrane and anti-

plane shear: derive and discuss the differential equations and the boundary conditions for

these two problems.

1.1.1 Taut membrane

Let us imagine to have a taut membrane subject to a pressure p(x, y) in direction z, as showed in the

following figure: Figure 1.1: Taut membrane problem

Where S is the tension per unit length and w(x, y) is the transversal displacement. On the boundary Γ

we impose the condition w = 0. The problem is studied under the following initial assumptions:

• ≃ ≃

small displacements and slopes =⇒ α sin α tan α

• tension S remains constant in the deformed configuration

Figure 1.2: Slope of the membrane

Advanced computational mechanics 1

1. Potential problems 1.1. Question 1

The slope of the membrane is equal to:

∂w ∂w ∂ ∂w

∂w α = α + dα = + dx = w + w dx

α = y ,x ,xx

x ∂x ∂y ∂x ∂x ∂x

Knowing the value of the slopes in the two planes we can write vertical equilibrium:

− −

S dy w + S dy (w + w dx) S dx w + S dx (w + w dy) + p dx dy = 0

,x ,x ,xx ,y ,y ,yy

p

2

∇

=⇒ S (w + w ) + p = 0 =⇒ w + = 0 in Ω

,xx ,yy S

The dynamic equation of taut membrane problem reads: 2

∂ w(x, y, t)

2

∇

S w(x, y, t) + p(x, y, t) = ρ 2

∂t

Where ρ is the mass per unit area. The condition w = 0 in Γ is the so-called kinematic condition. We

can have different kinematic boundary conditions like that one represented in the following figure.

Figure 1.3: Different boundary condition

In this case the membrane edge is free to slide vertically on boundary Γ . Writing the equilibrium of the

V

slider we get: V

∂w ∂w

− = in Γ

S ds + V ds = 0 =⇒ V

∂n ∂n S

This is the so-called Neumann boundary condition.

1.1.2 Anti-plane shear

Let us consider a homogeneous isotropic elastic cylinder with generators parallel to the z-axis and constant

cross section Ω. The cylinder is subjected to:

• prescribed surface tractions on its lateral boundary with nonzero component only in z-direction

• prescribed body forces, only in the z-direction

Figure 1.4: Anti-plane shear problem

Advanced computational mechanics 2

1. Potential problems 1.2. Question 2

The displacement field is equal to u = v = 0 and w = w(x, y). The only nonzero strain are:

γ (x, y) = w (x, y) γ (x, y) = w (x, y)

zx ,x zy ,y

The stresses are:

τ (x, y) = G γ (x, y) = G w (x, y) τ (x, y) = G γ (x, y) = G w (x, y)

zx zx ,x zy zy ,y

Writing the indefinite equilibrium equations we get:

σ + b = 0 =⇒ τ (x, y) + τ (x, y) + b (x, y) = 0 =⇒ G (w + w ) + b = 0

ij,i j zx,x zy,y z ,xx ,yy z

The boundary equilibrium equations are: ∂w

σ n = p =⇒ τ n + τ n = p =⇒ G = p

ij j i zx x zy y z z

∂n

Thus, the governing equations for the anti-plane shear problem are:

∂w p

z

2

∇

G w(x, y) + b (x, y) = 0 in Ω = in Γ

z p

∂n G

For a portion of the elastic cylinder of length t, the elastic energy is:

Z Z Z

t Gt

1 h i

2 2

σ ε dV = (τ γ + τ γ ) dΩ = (w ) + (w ) dΩ

= ij ij zx zx zy zy ,x ,y

2 2 2

V Ω Ω

The potential energy of the loads is: Z Z

−t −

= b w dΩ t p w dΓ

z z

Ω Γ

p

Therefore, the total potential energy is:

Z Z Z

Gt 2

|grad − −

π [w(x, y)] = w| dΩ t b w dΩ t p w dΓ

z z

2 Ω Ω Γ

p

1.2 Question 2

“Potential problems” governed by Laplace/Poisson equation: heat conduction and torsion

with stress function approach. Derive and discuss the differential equations and the bound-

ary conditions for these two problems.

1.2.1 Heat conduction

Let us consider a solid cylinder with generators parallel to the z-axis and constant cross section Ω. The

cylinder, homogeneous and thermally isotropic, is subjected to thermal external actions which do not

vary with z.

The objective of thermal analysis is to determine the temperature and heat flux within the solid body.

Here the focus is on steady-state problem (time-independent). In one-dimensional heat flow (along

x-axis) the heat flux is: Q

q =

x A ∆t

Where Q is the heat energy transferred through surface A in the direction of x during time interval ∆t.

The following relation holds for the fluxes: q = q n + q n

n x x y y

The Fourier’s law reads: ∂T

∂T −k −k ∇T

−k q = q =

q =

x y

∂x ∂y

Advanced computational mechanics 3

1. Potential problems 1.2. Question 2

Where k is the thermal conductivity. We write the energy balance for an infinitesimal element and for a

time interval dt:

∂ q ∂ q

x y

− + dx dy dt + g dx dy dt = c dm dT = c ρ dx dy dT

∂x ∂y

∂ q ∂T

∂ q y

x + + cρ = g

=⇒ ∂x ∂y ∂t

Substituting Fourier’s law we get: 2

∇

k T + g = c ρ Ṫ

For steady-state heat conduction Ṫ = 0 and the governing equation is:

2

∇

k T (x, y) + g(x, y) = 0

The boundary conditions can be:

• given temperature (Dirichlet boundary condition) T = T on Γ T

q

• ∂T −

given heat flux (Neumann boundary condition) = on Γ

n q

∂n k

• ∂T

−k −

given linear combination of temperature and flux (Robin boundary condition) = h (T T )

f

∂n

on Γ R

1.2.2 Torsion

Let us have a homogeneous, isotropic, elastic beam of constant cross section loaded only by two equal

and opposite twisting moments acting on the beam end cross sections. No restraints against the warping

of the originally plane sections are present at either beam ends.

Figure 1.5: Torsion problem

The displacement formulation is: −θ

u = y z v = θxz w = θ Ψ(x, y)

The strain are:

∂ Ψ ∂ Ψ

−y

γ = θ + γ = θ x +

yz

xz ∂x ∂y

The stress are:

∂ Ψ ∂ Ψ

−y

τ = G θ + τ = Gθ x +

xz yz

∂x ∂y

The indefinite equilibrium equations read:

2

∇

τ + τ = 0 =⇒ Ψ(x, y) = 0

xz,x yz,y ∂ Ψ

−

τ n + τ n = t = 0 =⇒ Ψ n + Ψ n = y n x n =⇒ = f (s) on Γ

xz x yz y z ,x x ,y y x y ∂n

For equilibrium: Z Z 2 2

− −

M = (τ x τ y) dΩ = G θ x + y + Ψ x Ψ y dΩ = G θ J

z zy zx ,y ,x t

Ω Ω

Advanced computational mechanics 4

1. Potential problems 1.3. Question 3

The problem can be solved by using Prandtl stress function. The initial assumptions are:

∂ F (x, y) ∂ F (x, y)

−

σ = σ = σ = 0 τ = 0 τ = τ =

x y z xy xz yz

∂y ∂x

To determine the value of the stress function we compare shear strain with those found with the displace-

ment approach: F F

,y ,x

−

γ = θ (−y + Ψ ) = γ = θ (x + Ψ ) =

xz ,x yz ,y

G G

And deriving by x and y: −F

G θ (−1 + Ψ ) = F G θ (1 + Ψ ) =

,xy ,yy ,xx ,xx

Subtracting the second equation from the first:

2

∇ −2G

F (x, y) = θ in Ω

The boundary equation gives: −

τ n + τ n = 0 =⇒ F y F (−x ) = F = 0 on Γ

xz x yz y ,y ,s ,x ,s ,s

This means that on the boundary F is constant and it can be chosen equal to zero.

Z Z Z Z

h i

− − − −

M = (τ x τ y) dΩ = (F x F y) dΩ = (F x) + (F y) dΩ + 2 F dΩ

z yz xz ,x ,y ,x ,y

Ω Ω Ω Ω

Using the divergence theorem:

Z Z Z

−

M = [F (x n + y n )] dΓ + 2 F dΩ = 2 F dΩ on Γ

z x y

Ω Ω Ω

Analogy with the heat conduction problem:

2 2

∇ F (x, y) + 2G θ = 0 k∇ T (x, y) + g(x, y) = 0

To solve the problem of torsion is sufficient to assume g = 2G θ and k = 1.

1.3 Question 3

Describe the generation of the discretized Finite Element equations for one of the (scalar)

potential problem dealt with in the course, using as starting point the stationarity of total

potential energy or the weak formulation of the problem (virtual displacement principle,

for elastostatic problems). Describe also the solution procedure.

Let us write the total potential energy for the taut membrane problem. The TPE depends on three

contributions:

• elastic energy associated to pretension S, which does not depend on w(x, y)

• contribution of elastic energy due to w (and stretching of the membrane at constant S), which can

be expressed as the product of S with the increase of element area.

1 1

q

q 2 2 2 2

≃ ≃

1 + (w ) dx 1 + (w ) dx ds = 1 + (w ) dy 1 + (w ) dy

ds = ,x y ,y

,x ,y

x 2 2

The energy contribution is:

1 1 S h i

2 2 2 2

· −

= S 1 + (w ) 1 + (w ) dx dy dx dy = (w ) + (w ) dx dy

,x

,x ,y ,y

2 2 2

• potential energy of the loads Z Z

− −

= p w dΩ V w dΓ

Ω Γ

V

Advanced computational mechanics 5

1. Potential problems 1.4. Question 4

Therefore, the total potential energy is:

Z Z Z

S 2

|grad − −

w| dΩ V w dΓ

π [w (x, y)] = p w dΩ

2 Ω Ω Γ

V

Discretizing the domain Ω in finite elements Ω we can write the total potential energy as follows:

n T

∂

∂w ∂w ∂w

2 T

n n n n,T n

|grad

w(x) = N (x) w = [N (x)] w w = w B(x) w

= B(x) =⇒ w| = B (x)

∂x ∂x ∂x ∂x

( )

Z Z Z

S

X T T T

n,T n n,T n,T

− −

dΩ w N V dΓ

π = B B w N p dΩ w =

w

2 Ω Ω Γ

n n

n Vn

1

X n n n

n,T n n,T

−

= w w + P

k w P p V

2

n

Assembling the elements of the domain we get:

1 1

n T T T T

n

− −

w = C W =⇒ π = W K W W P + P = W K W W P

p V

2 2

Imposing stationarity of total potential energy we find:

∂π −

= 0 =⇒ K W P = 0

∂W

1.4 Question 4

Explain the basis of the Lagrange multipliers method and illustrate a couple of applications

of the technique in computational mechanics (with comments on the mechanical interpre-

tation of the Lagrange multipliers). n n

⊂ → ⊂ →

Method of Lagrange Multipliers: suppose that f : U and g : U are given

R R R R

1 ∈

C U and g(x ) = c, and let S be the level set for g with value c (recall

real-valued functions. Let x 0 0

n

∈ ∇g(x ̸

that this is the set of points x f satisfying g(x) = c). Assume ) = 0. If f S, which denotes ”f

R 0

restricted to S”, has a local maximum or minimum on S at x , then there is a real number λ such that:

0

∇f ∇g(x

(x ) = λ )

0 0

±

We define Lagrangian the quantity L (x, f ) = f (x) λ g(x).

The method of Lagrange multipliers is a strategy for finding local maxima and minima of a function

subjected to equality constraints.

Example

For a function f (x, y) in two variables maximize f (x, y) subject to g(x, y) = c. Let us denote as (x , y )

0 0

the solution to this optimization problem. Intuively, at (x , y ) the constraint curve tangentially touches

0 0

∇f ∥ ∇g.

the contour line of f , therefore

Example on the sphere surface of radius R such that its distance from another point in the space

Find the point x

x is maximum/minimum.

0 h i

2 2 2 2 2 2

− − − −

min (x x ) + (y y ) + (z z ) subject to constraint x + y + z 1 = 0

0 0 0

The Lagrange function is: L (x, y, z, λ) = f (x, y, z) + λ g (x, y, z)

Advanced computational mechanics 6

1. Potential problems 1.4. Question 4

We impose stationarity: ∂L = 0 =⇒ x (1 + λ) = x 0

∂x

∂L = 0 =⇒ y (1 + λ) = y

0

∂y

∂L = 0 =⇒ z (1 + λ) = z

0

∂z

∂L 2 2 2 −

= 0 =⇒ x + y + z 1 = 0

∂λ

Solving the system we get:

2 2

20 y z

x 2 2

0 0 − |x | − ± |x |

+ + 1 = 0 =⇒ (1 + λ) = 0 =⇒ (1 + λ) =

0 0

2 2 2

(1 + λ) (1 + λ) (1 + λ)

And therefore: x y z

0 0 0

± ± ±

x = y = z =

| | |

|x |x |x

0 0 0

Example 2

A rectangular box without a lid is to be made from 12 m of cardboard. Find the dimensions x, y, z

of the box so as to maximize its volume. The function is the volume f = x y z and the constraint is the

−

surface g = x y + 2x z + 2y z 12 = 0. −

L = x y z + λ [x y + 2x z + 2y z 12]

∂L = y z + λ (y + 2z) = 0

∂x

∂L = x z + λ (x + 2z) = 0

∂y

∂L = x y + λ (2x + 2y) = 0

∂z

∂L −

= x y + 2x z + 2y z 12 = 0

∂λ

Solution: 3 2

−

x + 2z y 2z y z

yz − −

− =⇒ x z y z = 0 =⇒ x = y =⇒ = 0 =⇒ x = y = 2z λ =

λ = y + 2z y + 2z y + 2z 2

1

2 − −

=⇒ 12z 12 = 0 =⇒ z = 1 x = y = 2 λ = 2

Let us analyze the following problem both with finite elements and Lagrange multipliers.

Example Figure 1.6: Elastic structure with a rigid link

     

−k

2k 0 u 0

1

−k −k

2k u Q

=

2

     

−k

0 2k u Q

3

Adding the rigid link we have:      

−k

2k 0 1 u 0

1

−k −k

2k 0 u Q

2

