Stationary heat flow in 3D

FEM for stationary heat flow

Piotr Pluciński

e-mail: Piotr.Plucinski@pk.edu.pl

Jerzy Pamin

e-mail: Jerzy.Pamin@pk.edu.pl

Chair for Computational Engineering

Faculty of Civil Engineering, Cracow University of Technology URL: www.CCE.pk.edu.pl

Computational Methods, 2020 J.Paminc

Lecture contents

1 Stationary heat flow in 3D Model - strong formulation Model - weak formulation FE equations

2 Selection of approximation functions Shape functions for 1D problems Shape functions for 2D problems Shape functions for 3D problems

Stationary heat flow in 3D

q n

q_n = q^Tn

∇ =







∂

∂x

∂

∂y

∂

∂z







cold hot

q ∇T

Fundamental unknown - temperature T

Fourier’s law of heat conduction q = −D ∇T

heat flux density vector q = {q_x q_y q_z} [W/m²] temperature gradient gradT = ∇T = n

∂T

∂x

∂T

∂y

∂T

∂z

o

[K/m]

heat conduction matrix D = {k_ij} [W/(mK)]

Flux density grows with temperature gradient.

Heat flows from higher to lower temperature.

Computational Methods, 2020 J.Paminc

Heat energy balance

Heat generated equal to heat flowing out

Z

V

f dV = Z

S

q_ndS ∀ V

f – heat source density – energy supplied to the body per unit volume and time [W/m³]

V

S Using Green-Gauss-Ostrogradsky theorem about integration by parts

Z

S

q_ndS = Z

S

q^Tn dS = Z

V

divq dV = Z

V

 ∂q_x

∂x + ∂q_y

∂y + ∂q_z

∂z

 dV Z

V

f dV = Z

V

∇^Tq dV ∀ V ⇒ ∇^Tq = f ∀x ∈ V

Heat flow equations

Heat conduction equation (strong formulation)

∇^T(D∇T ) + f = 0 ∀x ∈ V + boundary conditions

q_n = q^Tn = bq on S_q– natural b.c. (Neumann) T = bT on S_T– essential b.c. (Dirichlet)

V

S_T S_q

Conduction matrix for isotropic materials D = kI

∂²T

∂x² + ∂²T

∂y² + ∂²T

∂z² + f

k = 0 – Poisson equation For isotropic materials without heat source

∂²T

∂x² + ∂²T

∂y² + ∂²T

∂z² = 0 – Laplace equation

Computational Methods, 2020 J.Paminc

Stationary heat flow in 3D

Weighted residual method Z

V

w ∇^T(D∇T ) + f
dV = 0 ∀w 6= 0 Z

V

w∇^T(D∇T )dV + Z

V

wf dV = 0

V

S_T S_q

Weak formulation

− Z

V

(∇w)^TD∇T dV + Z

S

w

D∇T

−q

^T

ndS + Z

V

wf dV = 0

− Z

V

(∇w)^TD∇T dV − Z

S

w q^Tn q_n

dS + Z

V

wf dV = 0 Z

V

(∇w)^TD∇T dV = − Z

Sq

w q dS −b Z

ST

w q_n dS + Z

V

wf dV ∀w natural b.c. secondary unknown

Stationary heat flow in 3D

Strong formulation

∇^T(D∇T ) + f = 0 ∀x ∈ V + boundary conditions

q_n = q^Tn = qb on S_q

T = Tb on S_T

V

S_T S_q

Weak formulation

Z

V

(∇w)^TD∇T dV = − Z

Sq

wqbdS − Z

ST

wq_ndS + Z

V

wf dV ∀w 6= 0

+ boundary condition

T = Tb on S_T

Computational Methods, 2020 J.Paminc

Stationary heat flow in 3D

Set of FE equations

T = NΘ – approximated temperature function N – shape function vector ((global approximation) Θ – nodal temperature (dof) vector

∇T = BΘ – approximated temperature gradient function B = ∇N – shape function derivative matrix

w = w^TN^T – approximation of weight function (∇w = Bw) Z

V

(∇w)^TD∇T dV = − Z

S_q

wbq dS − Z

S_T

wq_ndS + Z

V

wf dV ∀w^T

KΘ = f_b+ f K =

Z

V

B^TDBdV, f_b = − Z

Sq

N^Tq dS−b Z

ST

N^Tq_ndS, f = Z

V

N^Tf dV

Selection of approximation functions

Fundamental steps in FE procedure

1 Build a strong formulation of a problem

2 Convert the formulation into a weak format

3 Select approximation of unknown function

4 Select weighting function (usually Galerkin approach)

Computational Methods, 2020 J.Paminc

Selection of approximation functions in 1D

1 Strong formulation d

dx AkdT dx

!

+ f = 0

+ boundary conditions

q_x = qb at x_q ( e.g. x_q = 0) T = bT at x_T ( e.g. x_T = l)

x_q = 0 x_T = l

2 Conversion into weak formulation Z l

0

dw

dx AkdT dx

!

dx = (wA)     x=0

q − (wAqb _x)     x=l

+ Z l

0

wf dx + b.c.

T = bT at x_T ( e.g. xT = l)

Selection of approximation functions in 1D

3 Selection of approximation functions Linear function

T^e(x) = α₁+ α₂x = Φα^e Φ = [1 x], α^e =

 α₁ α₂



T^e(x) = N_i^e(x^e)T_i + N_j^e(x^e)T_j = N^eΘ^e N^e = [N_i^e(x^e) N_j^e(x^e)], Θ^e =

 T_i T_j



dT^e

dx^e = B^eΘ^e, B^e = dN^e dx^e =

"

dN_i^e dx^e

dN_j^e dx^e

#

x T

T_i Tj

x_i x_j l^e

x^e 0^e l^e 1

N_i^e(x^e)

x^e 0^e l^e

1 N_j^e(x^e)

Computational Methods, 2020 J.Paminc

Selection of approximation functions in 1D

3 Selection of approximation functions Quadratic function

T^e(x) = α₁+ α₂x + α₃x² = Φα^e Φ = [1 x x²], α^e =



 α₁ α₂ α₃





T^e(x) = N_i^e(x^e)T_i + N_j^e(x^e)T_j + N_k^e(x^e)T_k

= N^eΘ^e

N^e = [N_i^e(x^e) N_j^e(x^e) N_k^e(x^e)], Θ^e =



 T_i T_j T_k





dT^e

dx^e = B^eΘ^e, B^e = dN^e dx^e =

"

dN_i^e dx^e

dN_j^e dx^e

dN_k^e dx^e

#

x T

Ti

Tj

T_k

x_i x_j x_k l^e

x^e 0^e x^e_j l^e 1 N_i^e(x^e)

x^e 0^e x^e_j l^e

1 N_j^e(x^e)

x^e 0^e x^e_j l^e

1 N_k^e(x^e)

Selection of approximation functions in 2D

1 Strong formulation

∇^T(D∇T ) + f = 0 ∀x ∈ A + boundary conditions

q_n = q^Tn =qb on Γ_q

T = bT on Γ_T Γ_T

Γ_q

2 Conversion into weak formulation (h - configuration thickness)

Z

A

(∇w)^TDh∇T dA = − Z

Γ_q

whqdΓ −b Z

Γ_T

whq_ndΓ + Z

A

whf dA

+ boundary condition

T = bT on Γ_T

Computational Methods, 2020 J.Paminc

Selection of approximation functions in 2D

3 Selection of approximation functions Three-node element

T^e(x, y) = α₁+ α₂x + α₃y = Φα^e

Φ = [1 x y], α^e =



 α₁ α₂ α₃





T^e(x, y) = N_i^e(x^e, y^e)T_i + N_j^e(x^e, y^e)T_j+ + N_k^e(x^e, y^e)T_k = N^eΘ^e

N^e = [N_i^e(x^e, y^e) N_j^e(x^e, y^e) N_k^e(x^e, y^e)], Θ^e =



 T_i T_j T_k





y T (x, y)

x T_i

T_k

T_j i

k

j

y^e x^e

i k

j e

Selection of approximation functions in 2D

3 Selection of approximation functions Three-node element

N^e = [N_i^e(x^e, y^e) N_j^e(x^e, y^e) N_k^e(x^e, y^e)]

e.g. for N_i(x^e, y^e) N_i(x^e_i, y_i^e) = 1 N_i(x^e_j, y_j^e) = 0 N_i(x^e_k, y_k^e) = 0

y^e Ni(x^e, y^e)

x^e 1 i

k

j

y^e Nj(x^e, y^e)

x^e 1

i k

j y^e

Nk(x^e, y^e)

x^e

1 i

k

j

Determination of shape functions N_i(x^e, y^e) = α_1i+ α_2ix^e + α_3iy^e





1 x_i y_i 1 x_j y_j 1 xk yk







 α_1i α_2i α3i



=



 1 0 0



 =⇒

α1i = ^xjy_2P^k−xkyj

4

α_2i = ^y_2P^j−yk

4

α_3i = ^x_2P^k−xj

4

Computational Methods, 2020 J.Paminc

Selection of approximation functions in 2D

Pascal triangle – three-node element 1

x y

x² xy y²

x³ x²y xy² y³ x⁴ x³y x²y² xy³ y⁴

Pascal triangle – six-node element 1

x y

x² xy y² x³ x²y xy² y³ x⁴ x³y x²y² xy³ y⁴

Selection of approximation functions in 2D

Convergence conditions - requirements for FE approximation

completeness – approximation must be able to represent a constant field and a constant field gradient

continuity (conforming FE) – approximation must be continuous along interelement edges

Computational Methods, 2020 J.Paminc

Selection of approximation functions in 2D

3 Selection of approximation functions Four-node element

T^e(x, y) = α^e₁+ α^e₂x + α^e₃y + α^e₄xy = Φα^e

Φ = [1 x y xy], α^e =







 α^e₁ α^e₂ α^e₃ α^e₄









T^e(x, y) = N_i^e(x^e, y^e)T_i + N_j^e(x^e, y^e)T_j+

+ N_k^e(x^e, y^e)T_k + N_l^e(x^e, y^e)T_l = N^eΘ^e N^e = [N_i^e(x^e, y^e) N_j^e(x^e, y^e) N_k^e(x^e, y^e) N_l^e(x^e, y^e)]

Θ^e = {Ti T_j T_k T_l}

y T (x, y)

x T_j

T_k T_i

T_l j

i

k l

y^e x^e

j

i

k e l

Selection of approximation functions in 2D

3 Selection of approximation functions Four-node element (rectangular)

N^e = [N_i^e(x^e, y^e) N_j^e(x^e, y^e) N_k^e(x^e, y^e) N_l^e(x^e, y^e)]

e.g. for Ni(x^e, y^e) Ni(x^e_i, y_i^e) = 1 Ni(x^e_j, y_j^e) = 0 Ni(x^e_k, y^e_k) = 0 Ni(x^e_l, y_l^e) = 0

y^e Ni(x, y) = ^(x−xj_ab^)(y−yl)

x^e 1

i

j

k l b a

y^e Nl(x, y) = −^(x−xk_ab^)(y−yi)

x^e

i 1 j

k l

b a

y^e Nj(x, y) = −^(x−xi_ab^)(y−yk)

x^e 1 j i

k l b a

y^e Nk(x, y) = ^(x−xl_ab^)(y−yj)

x^e 1

i j

k l

b a

∇T^e = B^eΘ^e B^e =

" ∂N^e

∂x^e

∂N^e

∂y^e

#

Computational Methods, 2020 J.Paminc

Selection of approximation functions in 2D

Pascal triangle – four-node element 1

x y

x² xy y²

x³ x²y xy² y³ x⁴ x³y x²y² xy³ y⁴

Pascal triangle – eight-node element 1

x y

x² xy y² x³ x²y xy² y³ x⁴ x³y x²y² xy³ y⁴

Selection of approximation functions in 3D

1 Strong formulation

∇^T(D∇T ) + f = 0 ∀x ∈ V + boundary conditions

q_n = q^Tn =qb on S_q

T = bT on ST

V

S_T S_q

2 Conversion to weak formulation

Z

V

(∇w)^TD∇T dV = − Z

S_q

wqdS −b Z

S_T

wq_ndS + Z

V

wf dV + boundary condition

T = bT on S_T

Computational Methods, 2020 J.Paminc

Selection of approximation functions in 3D

3 Selection of approximation functions Tetrahedral element

T^e(x, y, z) = α^e₁+ α^e₂x + α^e₃y + α^e₄z T^e(x, y, z) = N_i^e(x^e, y^e, z^e)Ti + N_j^e(x^e, y^e, z^e)Tj

+ N_k^e(x^e, y^e, z^e)T_k + N_l^e(x^e, y^e, z^e)T_l

= N^eΘ^e

Hexahedral element

T^e(x, y, z) = α^e₁+ α^e₂x + α^e₃y + α^e₄z+

+ α^e₅xy + α₆^eyz + α^e₇xz + α^e₈xyz

y z

xi

j

k l

y z

x i

j k

l m

n o

p