Modelling steam methane reforming in a fixed bed reactor using - - PowerPoint PPT Presentation

Modelling steam methane reforming in a fixed bed reactor

using orthogonal collocation Kasper Linnestad

Department of Chemical Engineering Norwegian University of Science and Technology

December 3, 2014

Kasper Linnestad Reactor modelling December 3, 2014 1 / 36

Outline

1

Theory

2

Governing equations

3

Implementation

4

Results

5

Conclusion

Kasper Linnestad Reactor modelling December 3, 2014 2 / 36

Outline

1

Theory Weighted residuals method Collocation points Orthogonal collocation method

2

Governing equations

3

Implementation

4

Results

5

Conclusion

Kasper Linnestad Reactor modelling December 3, 2014 3 / 36

Theory

Weighted residuals method

General problem

L

• f(z, r)
• = g(z, r)

B

• fb(zb, rb)
• = gb(z, r)

Kasper Linnestad Reactor modelling December 3, 2014 4 / 36

Theory

Weighted residuals method

General problem

L

• f(z, r)
• = g(z, r)

B

• fb(zb, rb)
• = gb(z, r)

Approximation

f(z, r) ≈

Pz

• jz=0

Pr

• jr=0

ajz,jrljz(z)ljr(r)

Kasper Linnestad Reactor modelling December 3, 2014 4 / 36

Theory

Weighted residuals method

General problem

L

• f(z, r)
• = g(z, r)

B

• fb(zb, rb)
• = gb(z, r)

Approximation

f(z, r) ≈

Pz

• jz=0

Pr

• jr=0

ajz,jrljz(z)ljr(r)

Residual

R(z, r) = L

• f(z, r)
• − g(z, r)

Kasper Linnestad Reactor modelling December 3, 2014 4 / 36

Theory

Weighted residuals method

General problem

L

• f(z, r)
• = g(z, r)

B

• fb(zb, rb)
• = gb(z, r)

Approximation

f(z, r) ≈

Pz

• jz=0

Pr

• jr=0

ajz,jrljz(z)ljr(r)

Residual

R(z, r) = L

• f(z, r)
• − g(z, r)

Minimize residuals by

• R(z, r)wi(z, r) dz dr = 0,

∀i

Kasper Linnestad Reactor modelling December 3, 2014 4 / 36

Theory

Weighted residuals method

General problem

L

• f(z, r)
• = g(z, r)

B

• fb(zb, rb)
• = gb(z, r)

Approximation

f(z, r) ≈

Pz

• jz=0

Pr

• jr=0

ajz,jrljz(z)ljr(r)

Residual

R(z, r) = L

• f(z, r)
• − g(z, r)

Minimize residuals by

• R(z, r)wi(z, r) dz dr = 0,

∀i

Orthogonal collocation uses

ajz,jr = f(zjz, rjr)

Kasper Linnestad Reactor modelling December 3, 2014 4 / 36

Theory

Weighted residuals method

General problem

L

• f(z, r)
• = g(z, r)

B

• fb(zb, rb)
• = gb(z, r)

Approximation

f(z, r) ≈

Pz

• jz=0

Pr

• jr=0

ajz,jrljz(z)ljr(r)

Residual

R(z, r) = L

• f(z, r)
• − g(z, r)

Minimize residuals by

• R(z, r)wi(z, r) dz dr = 0,

∀i

Orthogonal collocation uses

ajz,jr = f(zjz, rjr) ln(x) =

P

• i=0

i=n

x − xi xn − xi

Kasper Linnestad Reactor modelling December 3, 2014 4 / 36

Theory

Weighted residuals method

General problem

L

• f(z, r)
• = g(z, r)

B

• fb(zb, rb)
• = gb(z, r)

Approximation

f(z, r) ≈

Pz

• jz=0

Pr

• jr=0

ajz,jrljz(z)ljr(r)

Residual

R(z, r) = L

• f(z, r)
• − g(z, r)

Minimize residuals by

• R(z, r)wi(z, r) dz dr = 0,

∀i

Orthogonal collocation uses

ajz,jr = f(zjz, rjr) ln(x) =

P

• i=0

i=n

x − xi xn − xi wi(z, r) = δ(z − ziz)δ(r − rir)

Kasper Linnestad Reactor modelling December 3, 2014 4 / 36

Theory

Collocation points

Legendre polynomials

1

−1

Ln(x)Lm(x) dx = 0, m = n −1 −0.5 0.5 1 −2 2 x Lm(x)

Kasper Linnestad Reactor modelling December 3, 2014 5 / 36

Theory

Collocation points

Legendre polynomials

1

−1

Ln(x)Lm(x) dx = 0, m = n

Collocation points

Lm(xi) = 0, ∀i ∈ {0, . . . , m} −1 −0.5 0.5 1 −2 2 x Lm(x)

Kasper Linnestad Reactor modelling December 3, 2014 5 / 36

Theory

Orthogonal collocation method

Minimize residuals by

• R(z, r)wi(z, r) dz dr = 0,

∀i

Orthogonal collocation uses

ajz,jr = f(zjz, rjr) ln(x) =

P

• i=0

i=n

x − xi xn − xi wi(z, r) = δ(z − ziz)δ(r − rir)

Kasper Linnestad Reactor modelling December 3, 2014 6 / 36

Theory

Orthogonal collocation method

Minimize residuals by

• R(z, r)wi(z, r) dz dr = 0,

∀i

Orthogonal collocation uses

ajz,jr = f(zjz, rjr) ln(x) =

P

• i=0

i=n

x − xi xn − xi wi(z, r) = δ(z − ziz)δ(r − rir)

Inserted

R(ziz, rir) = 0, ∀ iz ∈ {0, . . . , Pz} ir ∈ {0, . . . , Pr}

Kasper Linnestad Reactor modelling December 3, 2014 6 / 36

Theory

Orthogonal collocation method

Minimize residuals by

• R(z, r)wi(z, r) dz dr = 0,

∀i

Orthogonal collocation uses

ajz,jr = f(zjz, rjr) ln(x) =

P

• i=0

i=n

x − xi xn − xi wi(z, r) = δ(z − ziz)δ(r − rir)

Inserted

R(ziz, rir) = 0, ∀ iz ∈ {0, . . . , Pz} ir ∈ {0, . . . , Pr}

System of equations

L

• f(z, r)
• = g(z, r)

B

• fb(zb, rb)
• = gb(zb, rb)

Kasper Linnestad Reactor modelling December 3, 2014 6 / 36

Theory

Orthogonal collocation method

Minimize residuals by

• R(z, r)wi(z, r) dz dr = 0,

∀i

Orthogonal collocation uses

ajz,jr = f(zjz, rjr) ln(x) =

P

• i=0

i=n

x − xi xn − xi wi(z, r) = δ(z − ziz)δ(r − rir)

Inserted

R(ziz, rir) = 0, ∀ iz ∈ {0, . . . , Pz} ir ∈ {0, . . . , Pr}

System of equations

L

• f(z, r)
• = g(z, r)

B

• fb(zb, rb)
• = gb(zb, rb)

Linearisation

Af = b

Kasper Linnestad Reactor modelling December 3, 2014 6 / 36

Theory

Linearisation

Example function

f(ψ) = ψ2 + ψ exp (−ψ)

Kasper Linnestad Reactor modelling December 3, 2014 7 / 36

Theory

Linearisation

Example function

f(ψ) = ψ2 + ψ exp (−ψ)

Fixed-point (Picard-iteration)

f(ψ) ≈

• ψ⋆ + exp
• −ψ⋆

ψ solves f(ψ) = 3 in 114 iterations

Kasper Linnestad Reactor modelling December 3, 2014 7 / 36

Theory

Linearisation

Example function

f(ψ) = ψ2 + ψ exp (−ψ)

Fixed-point (Picard-iteration)

f(ψ) ≈

• ψ⋆ + exp
• −ψ⋆

ψ solves f(ψ) = 3 in 114 iterations

Taylor (Newton-Raphson-iteration)

f(ψ) ≈ f(ψ⋆) + ∂f ∂ψ

• ψ⋆
• ψ − ψ⋆

solves f(ψ) = 3 in 8 iterations

Kasper Linnestad Reactor modelling December 3, 2014 7 / 36

Outline

1

Theory

2

Governing equations Continuity equation Energy equation Species mass balance Ergun’s equation Initial and boundary conditions

3

Implementation

4

Results

5

Conclusion

Kasper Linnestad Reactor modelling December 3, 2014 8 / 36

Governing equations

Fixed bed reactor

z r

Kasper Linnestad Reactor modelling December 3, 2014 9 / 36

Governing equations

Steam methane reforming

Assumptions

Pseudo-homogeneous Efficiency factor = 10−3

Reactions

CH4 + H2O = CO + 3 H2 CO + H2O = CO2 + H2 CH4 + 2 H2O = CO2 + 4 H2

Kasper Linnestad Reactor modelling December 3, 2014 10 / 36

Governing equations

Continuity equation

General form

∂ρ ∂t + ∇ · (ρu) = 0

Kasper Linnestad Reactor modelling December 3, 2014 11 / 36

Governing equations

Continuity equation

General form

∂ρ ∂t + ∇ · (ρu) = 0

Simplified

∂ρuz ∂z = 0

Kasper Linnestad Reactor modelling December 3, 2014 11 / 36

Governing equations

Continuity equation

General form

∂ρ ∂t + ∇ · (ρu) = 0

Simplified

∂ρuz ∂z = 0

Linearised

ρ⋆ ∂uz ∂z + uz ∂ρ⋆ ∂z

• ⇒Auz

=

• ⇒buz

Kasper Linnestad Reactor modelling December 3, 2014 11 / 36

Governing equations

Energy equation

General form

ρcp ∂T ∂t + u · ∇T

• = −∇ · q − ∆rxH

Kasper Linnestad Reactor modelling December 3, 2014 12 / 36

Governing equations

Energy equation

General form

ρcp ∂T ∂t + u · ∇T

• = −∇ · q − ∆rxH

Simplified

ρcpuz ∂T ∂z = λeff r ∂ ∂r

• r ∂T

∂r

• − ∆rxH

Kasper Linnestad Reactor modelling December 3, 2014 12 / 36

Governing equations

Energy equation

General form

ρcp ∂T ∂t + u · ∇T

• = −∇ · q − ∆rxH

Simplified

ρcpuz ∂T ∂z = λeff r ∂ ∂r

• r ∂T

∂r

• − ∆rxH

Linearised

ρ⋆c⋆

pu⋆ z

∂T ∂z

• ⇒AT

= λeff,⋆

• 1

r ∂T ∂r + ∂2T ∂r2

• ⇒AT

− ∆rxH⋆

⇒bT

Kasper Linnestad Reactor modelling December 3, 2014 12 / 36

Governing equations

Species mass balance

General form

∂ρωi ∂t + ∇ · (ρuωi) = −∇ · ji + Ri

Kasper Linnestad Reactor modelling December 3, 2014 13 / 36

Governing equations

Species mass balance

General form

∂ρωi ∂t + ∇ · (ρuωi) = −∇ · ji + Ri

Simplified

∂ρuzωi ∂z = Deff r ∂ ∂r

• rρ ∂ωi

∂r

• − Ri

Kasper Linnestad Reactor modelling December 3, 2014 13 / 36

Governing equations

Species mass balance

General form

∂ρωi ∂t + ∇ · (ρuωi) = −∇ · ji + Ri

Simplified

∂ρuzωi ∂z = Deff r ∂ ∂r

• rρ ∂ωi

∂r

• − Ri

Linearised

ρ⋆u⋆

z

∂ωi ∂z + ωiu⋆

z

∂ρ⋆ ∂z + ρ⋆ωi ∂u⋆

z

∂z

• ⇒Aωi

= R⋆

i

• ⇒bωi

+ Deff,⋆

• ρ⋆

r ∂ωi ∂r + ∂ρ⋆ ∂r ∂ωi ∂r + ρ⋆ ∂2ωi ∂r2

• ⇒Aωi

Kasper Linnestad Reactor modelling December 3, 2014 13 / 36

Governing equations

Ergun’s equation

General form

dp dz = −f ρu2

z

dp

Kasper Linnestad Reactor modelling December 3, 2014 14 / 36

Governing equations

Ergun’s equation

General form

dp dz = −f ρu2

z

dp

Linearised

dp dz

• ⇒Ap

= −f⋆ ρ⋆(u⋆

z)2

dp

• ⇒bp

Kasper Linnestad Reactor modelling December 3, 2014 14 / 36

Governing equations

Initial and boundary conditions

Velocity

uz(z = 0, r)

• ⇒Auz

= uz,inlet

⇒buz

Kasper Linnestad Reactor modelling December 3, 2014 15 / 36

Governing equations

Initial and boundary conditions

Velocity

uz(z = 0, r)

• ⇒Auz

= uz,inlet

⇒buz

Pressure

p(z = 0, r)

• ⇒Ap

= pinlet

• ⇒bp

Kasper Linnestad Reactor modelling December 3, 2014 15 / 36

Governing equations

Initial and boundary conditions

Velocity

uz(z = 0, r)

• ⇒Auz

= uz,inlet

⇒buz

Temperature

T(z = 0, r)

• ⇒AT

= Tinlet

⇒bT

∂T ∂r

• r=0

⇒AT

=

• ⇒bT

λeff,⋆ ∂T ∂r

• r=R
• ⇒AT

= −U ⋆(T ⋆(r = R) − Ta)

• ⇒bT

Pressure

p(z = 0, r)

• ⇒Ap

= pinlet

• ⇒bp

Kasper Linnestad Reactor modelling December 3, 2014 15 / 36

Governing equations

Initial and boundary conditions

Velocity

uz(z = 0, r)

• ⇒Auz

= uz,inlet

⇒buz

Temperature

T(z = 0, r)

• ⇒AT

= Tinlet

⇒bT

∂T ∂r

• r=0

⇒AT

=

• ⇒bT

λeff,⋆ ∂T ∂r

• r=R
• ⇒AT

= −U ⋆(T ⋆(r = R) − Ta)

• ⇒bT

Pressure

p(z = 0, r)

• ⇒Ap

= pinlet

• ⇒bp

Mass fractions

ωi(z = 0, r)

• ⇒Aωi

= ωi,inlet

⇒bωi

∂ωi ∂r

• r=0
• ⇒Aωi

=

• ⇒bωi

∂ωi ∂r

• r=R
• ⇒Aωi

=

• ⇒bωi

Kasper Linnestad Reactor modelling December 3, 2014 15 / 36

Outline

1

Theory

2

Governing equations

3

Implementation Julia Under-relaxation Segregated approach Combined approach Coupled approach

4

Results

5

Conclusion

Kasper Linnestad Reactor modelling December 3, 2014 16 / 36

Implementation

Julia

Open-source Just-in-time compilation Constant global variables Pass by reference (!) Easy to call programs written in other languages (Fortran, C, Python)

Kasper Linnestad Reactor modelling December 3, 2014 17 / 36

Implementation

Under-relaxation

Trial-and-error ωi and T T = γT T + (1 − γT )T ⋆ γT = γω = 5 · 10−2

Kasper Linnestad Reactor modelling December 3, 2014 18 / 36

Implementation

Segregated approach

T ωi ||R||2 < tol uz p ||R||2 < tol Update Initial guess ||R||2 < tol

Kasper Linnestad Reactor modelling December 3, 2014 19 / 36

Implementation

Combined approach

Coupled mass fractions and temperature

        Aω1 · · · Aω2 · · · · · · ... I I · · · I · · · AT                 ω1 ω2 . . . ωCO T         =         bω1 bω2 . . . 1 bT        

Kasper Linnestad Reactor modelling December 3, 2014 20 / 36

Implementation

Combined approach

Coupled mass fractions and temperature

        Aω1 · · · Aω2 · · · · · · ... I I · · · I · · · AT                 ω1 ω2 . . . ωCO T         =         bω1 bω2 . . . 1 bT        

Coupled velocity, density and pressure

    Auz I −diag

• M⋆

RT⋆

• Ap

       uz ρ p    =    buz bp   

Kasper Linnestad Reactor modelling December 3, 2014 20 / 36

Implementation

Coupled approach

Fully coupled system of equations

               Aω1 · · · Aω2 · · · · · · ... I I · · · I · · · · · · AT · · · · · · Auz · · · I −diag

• M⋆

RT⋆

• · · ·

Ap                              ω1 ω2 . . . ωCO T uz ρ p               =               bω1 bω2 . . . 1 bT buz bp              

Kasper Linnestad Reactor modelling December 3, 2014 21 / 36

Outline

1

Theory

2

Governing equations

3

Implementation

4

Results Profiles Run-time Hindsights

5

Conclusion

Kasper Linnestad Reactor modelling December 3, 2014 22 / 36

Results

Velocity

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 z [m] uz

• m s−1

Axial velocity in the reactor r = 0,00 m r = 1,01 · 10−3 m r = 5,19 · 10−3 m r = 1,21 · 10−2 m r = 2,08 · 10−2 m r = 3,02 · 10−2 m r = 3,89 · 10−2 m r = 4,58 · 10−2 m r = 5,00 · 10−2 m r = 5,10 · 10−2 m Kasper Linnestad Reactor modelling December 3, 2014 23 / 36

Results

Velocity

0.00 0.01 0.02 0.03 0.04 0.05 1 2 3 4 5 6 7 2 2.5 3 3.5 r [m] z [m] uz

• m s−1

Velocity in the reactor Kasper Linnestad Reactor modelling December 3, 2014 23 / 36

Results

Pressure

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 2.76 · 106 2.78 · 106 2.80 · 106 2.82 · 106 2.84 · 106 2.86 · 106 2.88 · 106 2.90 · 106 z [m] p [Pa] Pressure in the reactor Kasper Linnestad Reactor modelling December 3, 2014 24 / 36

Results

Temperature

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 750 800 850 900 950 1,000 1,050 1,100 z [m] T [K] Temperature in the reactor r = 0,00 m r = 1,01 · 10−3 m r = 5,19 · 10−3 m r = 1,21 · 10−2 m r = 2,08 · 10−2 m r = 3,02 · 10−2 m r = 3,89 · 10−2 m r = 4,58 · 10−2 m r = 5,00 · 10−2 m r = 5,10 · 10−2 m Kasper Linnestad Reactor modelling December 3, 2014 25 / 36

Results

Temperature

0.00 0.01 0.02 0.03 0.04 0.05 1 2 3 4 5 6 7 800 900 1,000 r [m] z [m] T [K] Temperature in the reactor Kasper Linnestad Reactor modelling December 3, 2014 25 / 36

Results

Hydrogen

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 0.000 0.050 0.100 0.150 0.200 0.250 0.300 0.350 0.400 0.450 z [m] xH2 Mole fraction of H2 in the reactor r = 0,00 m r = 1,01 · 10−3 m r = 5,19 · 10−3 m r = 1,21 · 10−2 m r = 2,08 · 10−2 m r = 3,02 · 10−2 m r = 3,89 · 10−2 m r = 4,58 · 10−2 m r = 5,00 · 10−2 m r = 5,10 · 10−2 m

Kasper Linnestad Reactor modelling December 3, 2014 26 / 36

Results

Methane

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 0.040 0.060 0.080 0.100 0.120 0.140 0.160 0.180 0.200 0.220 z [m] xCH4 Mole fraction of CH4 in the reactor r = 0,00 m r = 1,01 · 10−3 m r = 5,19 · 10−3 m r = 1,21 · 10−2 m r = 2,08 · 10−2 m r = 3,02 · 10−2 m r = 3,89 · 10−2 m r = 4,58 · 10−2 m r = 5,00 · 10−2 m r = 5,10 · 10−2 m

Kasper Linnestad Reactor modelling December 3, 2014 27 / 36

Results

Carbon dioxide

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 0.000 0.010 0.020 0.030 0.040 0.050 0.060 z [m] xCO2 Mole fraction of CO2 in the reactor r = 0,00 m r = 1,01 · 10−3 m r = 5,19 · 10−3 m r = 1,21 · 10−2 m r = 2,08 · 10−2 m r = 3,02 · 10−2 m r = 3,89 · 10−2 m r = 4,58 · 10−2 m r = 5,00 · 10−2 m r = 5,10 · 10−2 m

Kasper Linnestad Reactor modelling December 3, 2014 28 / 36

Results

Steam

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 0.350 0.400 0.450 0.500 0.550 0.600 0.650 0.700 0.750 z [m] xH2O Mole fraction of H2O in the reactor r = 0,00 m r = 1,01 · 10−3 m r = 5,19 · 10−3 m r = 1,21 · 10−2 m r = 2,08 · 10−2 m r = 3,02 · 10−2 m r = 3,89 · 10−2 m r = 4,58 · 10−2 m r = 5,00 · 10−2 m r = 5,10 · 10−2 m

Kasper Linnestad Reactor modelling December 3, 2014 29 / 36

Results

Carbon monoxide

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 0.000 0.010 0.020 0.030 0.040 0.050 0.060 0.070 z [m] xCO Mole fraction of CO in the reactor r = 0,00 m r = 1,01 · 10−3 m r = 5,19 · 10−3 m r = 1,21 · 10−2 m r = 2,08 · 10−2 m r = 3,02 · 10−2 m r = 3,89 · 10−2 m r = 4,58 · 10−2 m r = 5,00 · 10−2 m r = 5,10 · 10−2 m

Kasper Linnestad Reactor modelling December 3, 2014 30 / 36

Results

Nitrogen

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 0.031 0.032 0.033 0.034 0.035 0.036 0.037 0.038 0.039 0.040 0.041 z [m] xN2 Mole fraction of N2 in the reactor r = 0,00 m r = 1,01 · 10−3 m r = 5,19 · 10−3 m r = 1,21 · 10−2 m r = 2,08 · 10−2 m r = 3,02 · 10−2 m r = 3,89 · 10−2 m r = 4,58 · 10−2 m r = 5,00 · 10−2 m r = 5,10 · 10−2 m

Kasper Linnestad Reactor modelling December 3, 2014 31 / 36

Results

Run-time

Solution method Run-time [s] Finite-difference & ode15s (Matlab) 1.8 Collocation – segregated 129.9 Collocation – combined 51.1 Collocation – coupled 52.0 Collocation – coupled & SparseMatrixCSC 14.8

Kasper Linnestad Reactor modelling December 3, 2014 32 / 36

Results

Hindsight

Non-dimensional equations Sparse matrices Higher collocation point density near inlet Linearise with Newton Iterative matrix solvers (CG, GMRES)

Kasper Linnestad Reactor modelling December 3, 2014 33 / 36

Outline

1

Theory

2

Governing equations

3

Implementation

4

Results

5

Conclusion

Kasper Linnestad Reactor modelling December 3, 2014 34 / 36

Conclusion

FDM with ode15s

Fast Simple Automatic step-size control No support for boundary conditions in z

Orthogonal collocation

Slow Not as simple (but not that hard) Supports boundary conditions in z Easy differentiation (Kronecker product) Trial and error for relaxation factors Global index notation Collocation points must be given beforehand

Kasper Linnestad Reactor modelling December 3, 2014 35 / 36

Questions?

Kasper Linnestad Reactor modelling December 3, 2014 36 / 36