General Balance Equation

n S

V

Consider an arbitrary fluid volume V, bounded by a surface S within a flow stream of G G velocity v . The volume is located at a position r and may be moving or stationary such G G that any point on S has velocity vs (rs ) ≠ 0 . Define ψ as some general property per unit volume of fluid, e.g. Mass → ψ = ρ (scalar) Energy E = Me → ψ = ρ e (scalar)

G G Momentum Mv → ψ = ρ v (vector) A balance on ψ gives Time rate of change of ψ within V= net flow rate of ψ into (out of) V across S plus generation rate of ψ within V. G G d ψ dV = − ψ ⋅ ndS + ψ g dV dt

∫

∫

V

∫

S

V

G Note: ψ = Net flux of ψ into V across S and contains a convective component and a surface flux component G

G

G G

ψ = ψ (v − vs ) +ψ s G G d G G G ψ dv = − ψ (v − vs ) ⋅ nds − ψ ⋅ ndS + ψ g dv dt

∫ V

∫ S

∫ S

∫ V

1

Leibnitz Rule ∂ψ

∫

G G

∫ ∂t dV + ∫ψ v ⋅ ndS

d ψ dV = dt

s

V

V

S

which may be substituted into the balance equation to give ∂ψ

G G

G G

G G

G

∫ ∂t dV + ∫ψ v ⋅ ndS = −∫ψ (v − v ) ⋅ nds − ∫ψ ⋅ ndS + ∫ψ dv s

V

∫ V

s

S

g

S

S

V

G G ∂ψ G G dV + ψ v ⋅ nds = − ψ ⋅ ndS + ψ g dv ∂t

∫ S

∫ S

∫ V

Gauss’s Theorem

G For any vector quantity B

∫ S

G G G G B ⋅ nds = ∇ ⋅ BdV

∫ V

such that the surface integrals can be eliminated giving

∫ V

G G G G ∂ψ dV + ∇ ⋅ψ vdV = − ∇ ⋅ψ s dV + ψ g dV ∂t

∫ V

∫ V

∫ V

or

∫ V

⎡ ∂ψ G G G G ⎤ ⎢⎣ ∂t + ∇ ⋅ψ v + ∇ ⋅ψ s −ψ g ⎥⎦dV = 0

For this to be true for any arbitrary function, and any arbitrary volume requires the integrand to be zero giving the final form of the general balance equation.

G G ∂ψ G G + ∇ ⋅ψ v = −∇ ⋅ψ s + ψ g ∂t The general balance equation provides the basis for the derivation of the fluid conservation equations.

2

Derivation of the Fluid Conservation Equations The conservation principles governing fluid flow are the same as those for any dynamic system, i.e. mass, energy and momentum. Mass Conservation (Continuity Equation)

ψ = ρ (Scalar) G ψ s = 0 (No surface sources of mass) ψ g = 0 ∂ρ G G + ∇ ⋅ ρv = 0 ∂t Momentum Conservation

G

ψ = ρ v (mass flux) G ψ g = Volumetric sources of Force → Weight = ρ g

GG ψ s = T Surface Stress Tensor containing both normal and shear stresses acting on S G

G G GG ∂ρ v G GG G + ∇ ⋅ ρ vv = −∇ ⋅ T + ρ g ∂t GG vv = diadic product

⎛ v1v1 v1v2 GG ⎜ vv = ⎜ v2 v1 v2 v2 ⎜v v v v ⎝ 31 3 2

v1v3 ⎞ iˆ ⎟ v2 v3 ⎟ ˆj v3v3 ⎟⎠ kˆ

GG GG GG T = PI − σ

ˆ GG ⎛ 1 0 0 ⎞ i GG ⎜ ⎟ I = ⎜ 0 1 0 ⎟ ˆj is the unit Tensor. The dot product of I and a vector returns the ⎜0 0 1⎟ ˆ ⎝ ⎠k vector.

G G GG GG ∂ρ v G GG G + ∇ ⋅ ρ vv = −∇ ⋅ ( P ⋅ I − σ ) + ρ g ∂t For any coordinate direction i ∈ ( x, y, z )

3

G G ∂ρ vi G G ∂P + ∇ ⋅ ρ vvi = − + ρ gi + ∇ ⋅ σ i ∂t ∂xi

or ∂ρ vi + ∂t

∑

∂ ( ρ v j vi )

j

∂x j

=−

G G ∂P + ρ gi + ∇ ⋅ σ i ∂xi

Energy Conservation

1 G G 2 G G ψ g = q′′′ + ρ v ⋅ g G G GG G ψ s = q′′ + T ⋅ v (Neglect Radiative Heat Transfer within the medium itself) G G q′′ = − k ∇T

ψ = ρu + ρ v ⋅ v

which yields the Total Energy Equation G G GG G ∂ ⎡ ⎛ 1 G G ⎞⎤ G ⎡ ⎛ 1 G G ⎞ G⎤ G G G ρ u v v ρ u v v v + ⋅ + ∇ ⋅ + ⋅ = −∇ ⋅ (q′′ + T ⋅ v ) + q′′′(r , t ) + ρ g ⋅ v ⎜ ⎟⎥ ⎜ ⎟ ⎥ ⎢ ⎢ 2 2 ∂t ⎣ ⎝ ⎠⎦ ⎠ ⎦ ⎣ ⎝

or G G G G G GG G ∂ ⎡ ⎛ 1 G G ⎞⎤ G ⎡ ⎛ 1 G G ⎞ G⎤ G G G ρ ⎜ u + v ⋅ v ⎟ ⎥ + ∇ ⋅ ⎢ ρ ⎜ u + v ⋅ v ⎟ v ⎥ = −∇ ⋅ q′′ − ∇ ⋅ Pv + ∇ ⋅ σ ⋅ v + q′′′(r , t ) + ρ g ⋅ v ⎢ 2 2 ∂t ⎣ ⎝ ⎠⎦ ⎠ ⎦ ⎣ ⎝

Unknowns G

GG

ρ , v , P, u , k , T , σ

Assume k = k(T) GG GG G σ = σ (T , v )

G Fundamental Variables: ρ , v , P, u, T ⇒ 7 Equations: Mass Momentum Energy

1 3 1 5

4

Two additional equations are then required for closure. Equations of State

ρ = ρ (u, T ) or ρ = ρ (T , P) u = u (T , P)

Additional Forms of the Energy Equation Kinetic Energy Equation Consider again the momentum equation for coordinate direction i G G ∂ρ vi G G ∂P + ∇ ⋅ ρ v ⋅ vi = − + ρ gi + ∇ ⋅ σ i ∂t ∂xi

Multiply by the scalar vi vi

G G G G ∂ρ vi ∂P + vi ∇ ⋅ ρ v ⋅ vi = −vi + vi ρ gi + vi∇ ⋅ σ i ∂t ∂xi

Note: vi

∂ρ vi 1 ∂ v 2 ∂ρ = ρ vi2 + i 2 ∂t 2 ∂t ∂t

Proof 1 ∂ ρ ∂vi2 vi2 ∂ρ 2 ρ vi = + 2 ∂t 2 ∂t 2 ∂t ∂v v 2 ∂ρ = ρ vi i + i ∂t 2 ∂t Also vi

∂ ∂ρ ⎤ ⎡ ∂v ρ vi = vi ⎢ ρ i + vi ⎥ ∂t ∂t ⎦ ⎣ ∂t ∂v ∂ρ = ρ vi i + vi2 ∂t ∂t

⇒ ρ vi

∂vi ∂ ∂ρ = vi ( ρ vi ) − vi2 ∂t ∂t ∂t

Substituting above gives 5

∂v v 2 ∂ρ 1 ∂ ρ vi2 = ρ vi i + i ∂t 2 ∂t 2 ∂t vi2 ∂ρ ∂ 2 ∂ρ = vi ( ρ vi ) − vi + ∂t ∂t 2 ∂t v 2 ∂ρ ∂ = vi ( ρ vi ) − i ∂t 2 ∂t ∂ρ vi 1 ∂ v 2 ∂ρ ⇒ vi = ρ vi2 + i 2 ∂t 2 ∂t ∂t

Similarly G G v2 G G 1G G vi ∇ ⋅ ρ vvi = ∇ ⋅ ρ vvi2 + i ∇ ⋅ ρ v 2 2

The momentum equation multiplied by vi is then G G vi2 ∂ρ vi2 G G 1 G G 2 1 ∂ ∂P 2 ( ρ vi ) + + ∇ ⋅ ρ v + ∇ ⋅ ρ vvi + vi − vi ∇ ⋅ σ i − vi ρ gi = 0 2 ∂t 2 ∂t 2 2 ∂x  

Ni 2 vi ⎛ d ρ G G ⎞ +∇⋅ ρ v ⎟ = 0 ⎜ 2 ⎝ dt ⎠

∂Pvi ∂v −P i ∂xi ∂xi

Examine the viscous term G G vi ∇ ⋅ σ i = vi

∂σ ij

∑ ∂x = ∑ v

i

j

vi

∂σ ij ∂x j

=

∂viσ ij ∂x j

j

− σ ij

j

∂σ ij ∂x j

∂vi ∂x j

Giving for the viscous term G G vi ∇ ⋅ σ i =

∑ j

vi

∂σ ij ∂x j

=

∑ j

∂ σ ij vi − ∂x j

∑ j

σ ij

∂vi G G = ∇ ⋅ σ i vi − ∂x j

∑σ j

ij

∂vi ∂x j

The momentum equation multiplied by vi may then be written as ∂Pvi ∂v G G 1 ∂ 1G G ρ vi2 + ∇ ⋅ ρ vvi2 + − P i − ∇ ⋅ σ i vi + 2 ∂t 2 ∂xi ∂xi

∑σ j

ij

∂vi − vi ρ gi = 0 ∂x j

for i ∈ (1, 2,3)

6

Sum over all values of i G G G GG G 1 ∂ G G 1 G GG G G G ρ v ⋅ v + ∇ ⋅ ρ vv ⋅ v + ∇ ⋅ Pv − P∇ ⋅ v − ∇ ⋅ σ ⋅ v + 2 ∂t 2

∑∑σ i

ij

j

∂vi G G − ρv ⋅ g = 0 ∂x j

which is the final form of the Kinetic Energy Equation. NOTE: The Kinetic Energy Equation was obtained by linear operations on the Momentum Equation and therefore is not linearly independent Internal Energy Equation Consider again the Total Energy Equation ∂ ⎡ ⎛ 1 G G ⎞⎤ G ⎡ ⎛ 1 G G ⎞ G ⎤ G G G G G GG G G G G ρ ⎜ u + v ⋅ v ⎟ ⎥ + ∇ ⋅ ⎢ ρ ⎜ u + v ⋅ v ⎟ v ⎥ + ∇ ⋅ q′′ + ∇ ⋅ Pv − ∇ ⋅ σ ⋅ v − q′′′(r , t ) − ρ g ⋅ v = 0 ⎢ 2 2 ∂t ⎣ ⎝ ⎠⎦ ⎠ ⎦ ⎣ ⎝ and the Kinetic Energy Equation

G G G GG G 1 ∂ G G 1 G GG G G G ρ v ⋅ v + ∇ ⋅ ρ vv ⋅ v + ∇ ⋅ Pv − P∇ ⋅ v − ∇ ⋅ σ ⋅ v + 2 ∂t 2

∑∑σ i

ij

j

∂vi G G − ρv ⋅ g = 0 ∂x j

Subtract the Kinetic Energy Equation from the Total Energy Equation giving the Internal Energy Equation G G G ∂ G G G ( ρ u ) + ∇ ⋅ ρ uv + ∇ ⋅ q′′ + P∇ ⋅ v − ∂t

∑∑σ i

j

ij

∂vi G − q′′′(r , t ) = 0 ∂x j

The term

∑∑ i

j

σ ij

∂vi ∂x j

represents the component of the viscous work term that results in internal energy (thermal energy). In comparison to the other terms in the internal energy equation, this term is generally insignificant and as such can normally be neglected. The Internal Energy Equation is also not linearly independent from the Total Energy Equation. This implies we are free to use either the Total Energy Equation or Internal Energy Equation as convenient. The Internal Energy Equation has computational advantages and will be the Energy Equation of choice for our remaining discussions

7

The fluids equations constitute a system of seven nonlinear, coupled, partial differential equations. In there general form, analytic solutions are not possible and numerical solutions can be extremely difficult. When modeling large engineering systems, such as reactor coolant systems, the overall system behavior is governed by the macroscopic fluid behavior. Flow systems tend to be primarily one dimensional with one dimensional segments connected at manifolds. The one-dimensional flow equations can be obtained by averaging the general fluids equations over the cross sectional flow area.

8

Area Averaged Single Phase Flow Equations

We start with the arbitrary flow geometry illustrated below. nz c a x

nc

z

n

θ

z – Distance along flow path (not elevation) a – Cross sectional area C – Interface between the fluid and the wall θ - Angle of the flow path relative to the horizontal G n - Unit outward normal vector to the surface at C G nz - Unit vector in the flow direction G nc - Outward normal to the plane G G (If a = constant then n = nc ) To derive the area averaged equations we make use of limiting forms of Leibnitz Rule and Gauss’s Theorem Leibnitz Rule ∂ ∂t

∫

G F (r , t )dS =

a

∫ a

∂F G G dC dS + Fvs ⋅ n G G n ⋅ nc ∂t

∫ c

Gauss’s Theorem

∫ A

G G ∂ ∇ ⋅ BdS = ∂z

∫ A

G G G G dC B ⋅ nz dS + B ⋅ n G G n ⋅ nc

∫ c

G G Special Case: B = nz

∫ a

G G ∂ G G ∂a G G dC G G dC nz ⋅nz dS + nz ⋅ n G G = ∇ ⋅ nz dS = + nz ⋅ n G G n ⋅ nc ∂z n ⋅ nc ∂z

∫ a

∫ C

∫ c

9

For the gradient operator projected along the flow direction G ∂ G G ∇ = nz (projects nz along flowpath) ∂z

∫

∂ ∂a G G dC (1)dS = 0 = + nz ⋅ n G G n ⋅ nc ∂z ∂z

−

∂a = ∂z

∫

a

c

∫ c

G G nz ⋅ n G G dC n ⋅ nc

Integrate the General Balance Equation over the flow area a. G G ∂ψ G G + ∇ ⋅ψ v = −∇ ⋅ψ s + ψ g ∂t 1)

∂ψ dS + ∂t

∫ a

2)

∫

G G ∇ ⋅ψ vdS =

3)

a

∫

4) G −∇ ⋅ψ s dS + ψ g dS

∫

a

a

Apply Leibnitz Rule and Gauss’s Theorem 1)

∂ G G dC ψ dS − ψ vs ⋅ n G G ∂t n ⋅ nc

∫

∫

c

2)

+

c

∂ G G G G dC ψ v ⋅ nz dS + ψ v ⋅ n G G ∂z n ⋅ nc

∫

∫

c

=

c

G G ndC ∂ G 3) − ψ s ⋅ nz ds − ψ s ⋅ G G + ∂z ns ⋅ nc

∫

∫

a

c

∫

4) ψ g dS a

G G G ndC ∂ ∂ ∂ G G G G G G ndC ψ ds + ψ v ⋅ nz dS + ψ (v − vs ) ⋅ G G = − ψ s ⋅ nz dS − ψ s ⋅ G G + ψ g dS ∂t ∂z n ⋅ nc ∂z n ⋅ nc

∫ a

∫ a

∫ c

∫ c

∫ c

∫ c

The integral around the interface (contour) can be broken into integrals over fixed and moving surfaces

10

∫= ∫ +∫ c

cN

cN → Fixed

cM → Moving

cM

G G G ndC ψ (v − vs ) ⋅ G G = 0 on the wall n ⋅ nc cN G G v ⋅ n = 0 (Fixed impermeable surface)

∫

∫

cM

G G G ndC G G G ψ (v − vs ) ⋅ G G (If moving surface is impermeable) (v − vs ) ⋅ n = 0 n ⋅ nc

so that the integrated balance equation is G G ndC ∂ ∂ ∂ G G G G ψ ds + ψ v ⋅ nz dS = − ψ s ⋅ nz dS − ψ s ⋅ G G + ψ g dS ∂t ∂z ∂z n ⋅ nc

∫

∫

a

∫

a

∫

c

c

∫ c

G G Define v ≡ v ⋅ nz < ψ >2 =

1 ψ dS (Area Average) a ∫a

Area Averaged Mass Equation

ψ =ρ G ψs = 0 ψ g = 0 ∂ ∂t

∫ a

ρ dS +

∂ ∂z

∂ < ρ >2 a ∂t a

a

∫

G G

ρ v ⋅ nz dS = 0

a

(for this problem

∂ ∂ < ρ >2 + ∂t ∂z

∫ ρ vdS = 0

∂ (a) = 0 ) ∂t

ρ v = G is the mass flux

c

∂ ∂ < ρ >2 + < ρ v >2 a = 0 ∂t ∂z

11

Area Averaged Momentum Equation G

ψ = ρv

GG GG ψ s = ( PI − σ ) G ψ g = ρ g G

∂ ∂t

∫ a

∂ G ρ vdS + ∂z

∫ a

G G GG G GG ndC ∂ ∂ GG G GG G GG ndC G ρ vv ⋅ nz dS = − (P Ι ) ⋅ nz dS − P Ι ⋅ G G + σ ⋅nz dS + σ ⋅ G G + ρ gdS ∂z n ⋅ nc ∂z n ⋅ nc

∫

∫

a

∫

c

∫

a

c

∫ a

G Project along the flow direction by dotting against nz ∂ ∂t

∂

∫ ρ vdS + ∂z ∫ ρ vvdS = a

a

G GG nG ⋅ nG GG nG ⋅ nG n ⋅ nz dC ∂ ∂ z − σ zz ⋅dS + σ ⋅ G G dC + σ ⋅ G Gz dC + ρ g z dS PdS − P G G + ∂z ∂z n ⋅ nc n ⋅ nc n ⋅ nc

∫ a

∫

∫

c

a

∫

∫

cN

cM

∫ a

G G ∂ ∂ ∂ n ⋅ nz ∂ a < ρ v > 2 + < ρ vv > 2 a = − < P2 > a − P G G dC + < σ zz > 2 a ∂t ∂z ∂z n ⋅ nc ∂z

∫ c

− < τ w >1 Pw + aΔPp′ + < ρ > 2 ag z

Examine the pressure integral

∫

G G n ⋅ nz P G G dC =< P >1 n ⋅ nc

∫

G G n ⋅ nz P G G dC =< P >1 n ⋅ nc

c

c

∫

G G n ⋅ nz G G dC n ⋅ nc

∫

G G n ⋅ n2 ∂a ∂a ≈ − < P >2 G G = − < P1 > n ⋅ nc ∂z ∂z

c

c

Giving

∂ ∂ ∂ ∂a ∂a < ρ v >2 + < ρ vv > 2 a = −a < P > 2 − < P > 2 + < P >2 + ∂t ∂z ∂z ∂z ∂z ∂ < σ zz >2 a − < τ w >1 Pw + ΔPp′a+ < ρ > 2 ag z ∂z

a

12

or simplifying a

∂ ∂ ∂ ∂ < ρ v >2 + < ρ vv >2 a = −a < P >2 + < σ zz > 2 a− < τ w >1 Pw + ΔPp′a + < ρ > 2 ag z  

∂t ∂z ∂z ∂z −< ρ >2 ag sin θ

Area Averaged Internal Energy Equation G G G G G ∂ G ( ρ u ) + ∇ ⋅( ρ uv ) = −∇ ⋅ q′′ − P∇ ⋅ v + q′′′ + ∂t

∑∑σ j

i

ij

∂vi ∂x j

Integrate over area applying Leibnitz Rule and Gauss’s Theorem G ∂ G G ∂ ∂ G G G G ndC ( ρ u )dS + ρ uv ⋅ nz dS + ρ u( v − vs ) ⋅ G G = − q′′ ⋅ nz dS ∂z ∂t ∂z n ⋅ nc a a c a G G G ∂v G ndC − q′′ ⋅ G G − P∇ ⋅ vdS + q′′′dS + σ ij i dS ∂x j n ⋅ nc j i

∫

∫

∫

∫

∫

c

∫

∫ ∑∑

∫

a

a

a

∂ ∂ < ρ u > 2 + < ρ uv > 2 a = ∂t ∂z G G ∂ ′′ >1 Pw − P∇ ⋅ vdS + < q′′′ > 2 a + − < q′′z > 2 a + < qwall ∂z

a

∫ a

∫ ∑∑σ a

i

j

ij

∂vi dS ∂x j

G G If q′′ = − k ∇T ∂T G G q′′ ⋅ nz = q′′z − k ∂z

∫

G G P∇ ⋅ vdS =< P > 2

a

∫

G G ∇ ⋅ vdS

a

∂ G G G G dS =< P > 2 v ⋅ nz dS + v ⋅ n G G ∂z n ⋅ nc

∫ a

∫ c

Note:

∫ c

G G dS v ⋅n G G = n ⋅ nc

∫ c

G G G ndS (v − vs ) ⋅ G G + n ⋅ nc

∫

cM

G G dS vs ⋅ n G G n ⋅ nc

13

∫

⎧ ⎫ G G G G dS ⎪ ⎪∂ G G P∇ ⋅ vds =< P > 2 ⎨ v ⋅ n z + vs ⋅ n G G ⎬ n ⋅ nc ⎪ ⎪⎩ ∂z a cM ⎭

a

∂ ∂ ∂ ∂ < ρ u > 2 + < ρ uv > 2 a = − < q′′z > a + < qw′′ >1 Pw − < P > 2 < v > 2 a − Wex′ + < q′′′ > 2 a ∂t ∂z ∂z ∂z

a

∫

∫

We make the further assumptions: < P > 2 =< P > 2 < abc >≈< a >< b >< c >

Given: ρ = ρ (u, P)

⇒ < ρ > 2 = ρ (< u > 2 , < P > 2 )

For convenience we will also drop the average notation, such that the area averaged equations are Mass a

∂ρ ∂ + ( ρ va) = 0 ∂t ∂z

Momentum a

∂ρ v ∂ ∂P ∂σ zz a − τ w Pw + ΔPp′a − ρ ga sin θ + ( ρ vva) = −a + z

∂t ∂z ∂z  ∂ ≅0

Internal Energy a

∂ρ u ∂ ∂q′′ ∂va  + ( ρ uva) = − z a + qw′′ Pw − P − Wex′ + q′′′a ∂t ∂z ∂z ∂z N ≅0

State

ρ = ρ (u, P) The 1-D area averaged equations are applicable to any 1-D piping segment. These 1-D segments are assumed to be connected at manifolds. Manifolds have multiple inlets and outlets, and as such the 1-D equations do not apply. To address manifolds we average the general fluids equations over volume.

14

Volume Averaged Equations

The volume V is bounded by a surface S which consists of flow surfaces (Aj), fixed surfaces (SW) and moving surfaces (SM). The volume averaged equations are obtained by integrating over V. Mass Conservation ∂ρ G G + ∇ ⋅ ρv = 0 ∂t

∫ V

G G dρ dV + ∇ ⋅ ρ vdV = 0 dt

∫ V

Apply Leibniz Rule and Gauss’s Theorem G G

G G

d dt

∫ ρdV − ∫ ρ ⋅ v ⋅ ndS + ∫ ρ v ⋅ ndS = 0

d dt

∫

s

V

S

∫

S

G G

G

ρ dV + ρ ( v − vs ) ⋅ ndS

V

S

∫ ρ dV = M =< ρ > V 3

V

G G

G

∫ ∫

ρ v ⋅ ndS = m j =< ρ v j > 2 Aj

S

∫

G G

ρ ( v − vs ) ⋅ ndS =

SW

G G

G

ρ ( v − vs ) ⋅ n dS +

∫

SM

G G

G

ρ ( v − vs ) ⋅ n dS +

∑ ∫ ρ (v − v ) ⋅ ndS G

G

G

s

j

Aj

G G Note: v j = v ⋅ n → (+ ) outgoing (-) incoming

A

15

The volume averaged mass equation is then Mass V

d < ρ >3 + dt

∑ < ρv

j

> 2 Aj = 0

j

or equivalently d M+ dt

∑ m − ∑ m

min

ex

ex

=0

in

It is also common to assume < ρ >3 ≅< ρ > 2 . Internal Energy G G G G G ∂ G ( ρ v) + ∇ ⋅ ( ρ uv ) = −∇ ⋅ q′′ − P∇ ⋅ v + q′′′ + ∂t

∑∑σ i

ij

j

∂vi ∂x j

Integrate over volume applying Leibnitz Rule and Gauss’s Theorem d dt

∫

G G

∫

G

G G

∫

∫

G G

G

∫

ρudV + ρu (v − vs ) ⋅ ndS = − P∇ ⋅ vdV − q′′ ⋅ ndS + q′′′(r )dV +

V

S

V

S

V

∫ ∑∑σ V

i

ij

j

∂vi dV ∂x j

∫ ρu(v − v ) ⋅ ndS = ∫ ρu(v − v ) ⋅ n dS + ∫ ρu(v − v ) ⋅ n dS + ∑ ∫ ρu(v − v ) ⋅ ndS = ∑ < ρ uv > A G G

G

G G

s

G

G G

s

S

G

G

s

SW

G

s

j

SM

j

G

Aj

j

j

∫

G G q′′ ⋅ ndS = qw = − hc As (Tw − T∞ ) averaged over the entire wall surface

∫

G q′′′(r )dV = q g

S

V

giving for the integrated internal energy equation V

d < ρ u >3 + dt

∑ j

G G < ρuv j > 2 Aj = − P∇ ⋅ vdV + qw + q g +

∫ V

∫ ∑∑σ V

j

i

ij

∂vi dV ∂x j

16

∫ ∑∑σ V

∫

j

ij

i

∂vi dV = W Dis (Dissipative work term) ∂x j

G G G G P∇ ⋅ vdV =< P >3 ∇ ⋅ vdV

V

∫ G G ∫ v ⋅ ndS V

=< P >3

S

⎡ G G G G G ⎤ =< P >3 ⎢ ∫ (v − vs ) ⋅ ndS + ∫ vs ⋅ ndS ⎥ S ⎣S ⎦ ⎡ ⎤ G G G G ⎥  ⎢ =< P >3 vs ⋅ ndS + v ⋅ ndS =< P >3 ⎢ ⎥ j A 1 ⎣ Sm ⎦ ⎡ ⎤ = < P >3 ⎢ < v j > 2 Aj ⎥ + Wex ⎢⎣ j ⎥⎦

∑∫

∫

⎡ ⎢ ⎢⎣

∑ j

⎤ < v j > 2 Aj ⎥ + Wex ⎥⎦

∑

Note: < P >3 ≈< P >3 because δ P are negligible V

d < ρ u >3 + dt

∑ < ρuv

j

> 2 Aj = − < P >3

j

∑ ⎡⎣< v

j

> 2 Aj ⎤⎦ − Wex + q w + q g + W DIS

j

Again, for convenience we will drop the average notation and assume < abc >≈< a >< b >< c >

< ρ >3 = ρ (< u >3 , < P >3 )

such that the volume averaged equations are Mass V

d ρ+ dt

∑ ρv A = 0 j

j

j

Internal Energy V

d ρu + dt

∑ (ρuv) A = −P∑ ⎡⎣v A ⎤⎦ − W j

j

j

j

j

ex

+ qw + q g + W DIS

j

17

Momentum is a vector thus volume averaged momentum equations must be analyzed in all directions. For a generic volume this becomes nontrivial. Momentum exchange across a manifold (pressure losses) will be handled through experimentally determined loss coefficients.

18

Systems Analysis

n = 2, ..., N

n=1 Upper Plenum

Reactor

Lower Plenum Pump

The system can be treated as a flow network, where multiple 1-D piping segments are connected at manifolds. The 1-D segments are assumed to be of constant cross sectional area. Area Averaged Mass Equation Ax

∂ρ ∂ + ρ vAx = 0 ∂t ∂z

Area Averaged Internal Energy Equation Ax

∂vAx ∂ ∂ ρ u + ( ρ uvAx ) = − P + qw′′ Pw + q′′′( z ) Ax ∂t ∂z ∂z

Area Averaged Momentum Equation ∂ 1 ∂ ∂P τ w Pw ( ρ vvAx ) = − ρv + − − ρ g sin θ + ΔPp′ ∂t ∂z Ax ∂z Ax

State Equations

ρ = ρ (u, P)

19

u = u (T , P)

Constitutive Relations qw′′ = hc (Tw − T )

τ w Pw Ax

=

f ρ v2 + De 2

∑

K jδ ( z − z j )

ρ v2

j

2

Assumptions: The system is a single phase liquid Global compressibility – Pressure changes around the loop are sufficiently small that the pressure distribution is unimportant and all state equations can be evaluated at system pressure. Because the velocity of the fluid is much greater than the rate of temperature increase, ρ ≅ constant in the mass conservation equation. We employ the Bousinessq Approximation ⇒ density variations only important in the buoyancy terms ( ρ g sin θ ). In addition, we assume that where variations in density are important, the fluid is incompressible with respect to pressure but thermally expandable, i.e. ρ = ρ (T ) . Mass ∂ (GAx ) = 0 ∂z

Internal Energy Ax ρ

∂ρ vAx ⎤ ∂vAx ∂u ∂u ⎡ ∂ + qw′′ Pw + q′′′( z ) Ax + u ⎢ Ax ρ + + ρ vAx = −P ⎥ ∂t  t  ∂z ⎦ ∂z ∂z ⎣ ∂ 0

Ax ρ

∂u ∂u + GAx = qw′′ Pw ∂t ∂z

Momentum ∂P f G 2 ∂G ∂ ⎧ G 2 ⎫ = − − − + ⎨ ⎬ ∂z De 2 ρ ∂t ∂z ⎩ 2 ρ ⎭

∑ j

K jδ ( z − z j )

G2 − ρ g sin θ + ΔPρ′ 2ρ

State

ρ = ρ (u, P) = ρ (u ) or ρ = ρ (T )

20

Loop Momentum Balance Approach Integrate the momentum equation around the flow loops.

∫⇒ ∑∫ i

(break into integrals over sections of constant Ax )

Li

Momentum ∂G −∂P f G 2 = − − ∂t ∂z De 2 ρ

∑ j

G2 − ρ g sin θ + ΔPp′ kiδ ( z − zi ) 2ρ

Loop 1 ∂G dz = − Li ∂t i∈L

∑∫ i∈L1

∑ i

∂P dz − Li ∂z

∑∫ 1

dG Li 1i = − {PL. P. − PU . P.} − dt

∑∫ i∈L1 L i

∑ i

f G2 dz − De 2 ρ

∫∑ L1 j∈L1

f1i G12i − Dei 2 ρ

∑K

f ni Gni2 − Dei 2 ρ

∑

G2 K jδ ( z − z j ) dz − ρ g sin θ dz + ΔPp′dz 2ρ

∫ L1

G12j

L1

∫

− − ρ g sin θ dz + ΔPP1

2ρ

j

j∈L1

∫

L1

Similarly for loops n = 2,…,N

∑ i

Li

dGni = − {PL. P. − PU .P.} − dt

∑ i

Kj

j∈Ln

Gnj2

2ρ

∫

− − ρ g sin θ dz + ΔPPn Ln

and from the lower plenum through the reactor to the upper plenum Lc

dGc f L G2 = −{PU . P. − PL. P.} − c c c − dt Dei 2 ρ

∑ j∈Lc

Kj

Gc2 − ρ gdz 2ρ

∫ Lc

We can further take advantage of mass conservation to eliminate mass flux in terms of mass flow rate m n = Gni Ani

Loops n= 1,…,N

∑ i

∑

Li dm n f ni Li 1 mn2 = −{PL. P. − PU . P.} − − Ani dt Dei Ani2 2 ρ i∈L n

∑ j∈Ln

K j mn2 − ρ g sin θ dz + ΔPpn Anj2 2 ρ

∫

Ln

21

Core Lc dm c f c Lc 1 mc2 = −{PU . P. − PL. P.} − − Ac dt Dec Ac2 2 ρ

∑ j∈Lc

K j mc2 − ρ gdz Ac2 2 ρ

∫ Lc

In general the number of unknowns is m n , m c , PU . P. , PL. P. ⇒ N + 3

We have N+1 momentum equations, and the mass conservation equation m c =

∑ m ⇒ N + 2 Equations implying one addition equation is needed for closure. n

n

Since we are not interested in the pressure distribution, we can reduce the number of unknowns by one by defining ΔPc = PL. P. − PU . P. such that the momentum equations can be written Loops n= 1,…,N

∑ i

Li dm n = −ΔPc − Ani dt

∑ i∈Ln

f ni Li 1 mn2 − Dei Ani2 2 ρ

∑ j∈Ln

K j mn2 − ρ g sin θ dz + ΔPpn Anj2 2 ρ

∫

Ln

Core Lc dm c f L 1 m2 = ΔPc − c c 2 c − Ac dt Dec Ac 2 ρ

∑ j∈Lc

K j mc2 − ρ gdz Ac2 2 ρ

∫ Lc

such that the integrated momentum equations and the mass conservation equation give N+2 equations in the N+2 unknowns m n , m c , ΔPc ⇒ N + 2 Special Case: One asymmetric loop

Assume that only Loop 1 behaves asymmetrically, i.e. m 2 = m 3 = ... = m N . We require then only 3 momentum equations

∑ i

Li dm n = −ΔPc − Ani dt

∑ i∈Ln

f ni Li 1 mn2 − Dei Ani2 2 ρ

∑ j∈Ln

K j mn2 − ρ g sin θ dz + ΔPpn n=1,2 Anj2 2 ρ

∫

Ln

22

Lc dm c f c Lc 1 mc2 = ΔPc − − Ac dt Dec Ac2 2 ρ

∑ j∈Lc

K j mc2 − ρ gdz Ac2 2 ρ

∫ Lc

and the mass conservation equation m 1 + ( N − 1)m 2 = m c

to describe the system. Special Case: All loops behave identically

Solution requires the Loop 1 momentum equation

∑ i

Li dm 1 = −ΔPc − A1i dt

∑ i∈L1

f1i Li 1 m12 − Dei A12i 2 ρ

∑ j∈L1

K j m12 − ρ g sin θ dz + ΔPp1 A12j 2 ρ

∫ L1

the core momentum equation Lc dm c f c Lc 1 mc2 = ΔPc − − Ac dt Dec Ac2 2 ρ

∑ j∈Lc

K j mc2 − ρ gdz Ac2 2 ρ

∫ Lc

and mass conservation m c = Nm 1

In addition, we can use the mass conservation equation to eliminate m 1 in favor of m c and add the momentum equations to yield a single momentum equation in the system flow rate

∑ i

Li dm c =− Ai dt

∑ i

fi Li m c2 1 − Dei 2 ρ Ai2

∑ j

K j m c2 − A2j 2 ρ

v∫ ρ g sin θ dz + ΔP

p1

Analysis of Buoyancy Term Each of the momentum equations contains an integral of the body force term. We can approximate the integral as

∫

Ln

= ∑ ρ k g ΔH k → ΔΗ k is the change in elevation of that node and k

ρ k = ρ (Tk )

23

where Tk (or uk ) can be obtained by solving the energy equation around the loop. Let m 2 = m | m | such that in the event of flow reversal friction forces always oppose the flow. Consider the time advancement scheme

∑ i

Li ⎧ m nt +Δt − m nt ⎫ t +Δt ⎨ ⎬ = −ΔPc − Δt Ani ⎩ ⎭

∑ i

f ni Lni 1 m nt +Δt | m nt +Δt | − 2ρ Deni Ani2

∑ j

K j m nt +Δt | m nt +Δt | − 2ρ Anj2

∑ ρ gΔΗ t k

k

+ ΔPpnt

k

( N − 1)m 2t +Δt + m1t +Δt = m ct +Δt Note: The non derivative terms can be evaluated at either new or old time with no effect on numerical stability. The time level specification shown here is chosen for robustness, with the short time scale phenomenon evaluated at new time and long time scale phenomenon evaluated at past time. Analyze the stability of this time advancement scheme The integrated momentum equations are effectively discretized forms of the differential equations ∂G =0 ∂z ∂G ∂P =− ∂t ∂z

on a mesh

t G t +Δt − G tj+Δ ∂G −1/2 → j +1/ 2 =0 ∂z Δz

∂G j +1/2

⎧ P − Pj ⎫ = − ⎨ j +1 ⎬ ∂t ⎩ Δz ⎭ t t t +Δt t +Δt G tj+Δ ⎪⎧ P − Pj ⎪⎫ +1/2 − G j +1/2 = − ⎨ j +1 ⎬ Δt Δz ⎪⎩ ⎪⎭

24

Expressing each of the parameters as a perturbation about a fixed reference value leads to t t +Δt δ G tj+Δ +1/2 − δ G j −1/ 2

Δz t t δ G tj+Δ +1/2 − δ G j +1/2

Δt

=0 t +Δt t +Δt ⎪⎧ δ Pj +1 − δ Pj ⎪⎫ = −⎨ ⎬ Δz ⎪⎩ ⎪⎭

Express each parameter as a Fourier component of the form

δψ tj = δψ t e

{

δ G t +Δt e δ G t +Δt e

δG

t +Δt

ikz j

ikz j +1/2

−e

ikz j −1/2

− δ Gt e Δt

ikz j +1/2

}=0

ikz j +1/2

=

−δ P t +Δt ikz j+1 − kzj e −e Δz

{

}

Δz Δz − ik ik ⎫ −Δt t +Δt ⎧ 2 −δG = δ P ⎨e − e 2 ⎬ Δz ⎩ ⎭

t

⎛ Δz ⎞ ⎛ Δz ⎞ ⎡ ⎛ Δz ⎞ ⎛ Δz ⎞ ⎤ ⎛ Δz ⎞ cos ⎜ k ⎟ + i sin ⎜ k ⎟ − ⎢ cos ⎜ k ⎟ − i sin ⎜ k ⎟ ⎥ = 2i sin ⎜ k ⎟ = ik ′ ⎝ 2 ⎠ ⎝ 2 ⎠ ⎣ ⎝ 2 ⎠ ⎝ 2 ⎠⎦ ⎝ 2 ⎠

δ G t +Δt = 0 δ G t +Δt − δ G t =

−Δt δ P t +Δt ik ′ Δz

or

δ G t +Δt +

Δt δ P t +Δt ik ′ = δ G t Δz

This can be written in matrix form as

AU t+Δt = BU t U t+Δt = A -1BU t

U t+Δt = λ U t 25

where λ are the eigenvalues of A -1B . Stability requires the spectral radius (largest eigenvalue in absolute value) be less than or equal to one. The eigenvalues of A -1B are found by | Aλ − Β |= 0 0 ⎞ ⎛1 ⎜ ⎟ A= ⎜1 Δt ik ′ ⎟ ⎝ Δz ⎠

⎛0 0⎞ B=⎜ ⎟ ⎝1 0⎠ 0 ⎞ ⎛ λ ⎜ ⎟ | Aλ − Β |= ⎜ λ − 1 λ Δt ik ′ ⎟ ⎝ Δz ⎠

This assumed no perturbations in G via the mass equations, i.e. perturbations do not grow in this system. det( Aλ − Β ) = λ 2

Δt ik ′ = 0 Δz

λ = 0 → Unconditionally stable

26

Numerical Solution of the Loop Momentum Balance Equations For illustration purposes, we consider the one asymmetric loop problem, such that the time discretized loop momentum balance equations can be written as Loop 1

∑ i∈L1

Li ⎧ m 1t +Δt − m 1t ⎫ t +Δt ⎨ ⎬ = −ΔPc − Δt A1i ⎩ ⎭

∑ i∈L1

f1ti +Δt L1i 1 m 1t +Δt | m 1t +Δt | − 2ρ De1i A12i

∑ j∈L1

K j m 1t +Δt | m 1t +Δt | − ρ kt g ΔΗ k + ΔPpt1 2ρ A12j k ∈L

∑ 1

Loop 2 Li ⎧ m 2t +Δt − m 2t ⎫ f 2ti+Δt L2i 1 m 2t +Δt | m 2t +Δt | t +Δt − ⎨ ⎬ = −ΔPc − Δt 2ρ De2i A12i ⎭ 2i ⎩ i∈L2

∑A i∈L2

∑

K j m 2t +Δt | m 2t +Δt | − ρ kt g ΔΗ k + ΔPpt 2 2 2 ρ nj k∈L

∑A j∈L2

∑ 2

Core

Lc Ac

t +Δt t +Δt ⎧ m ct +Δt − m ct ⎫ f ct +Δt Lc 1 m c m c t +Δt − ⎨ ⎬ = ΔPc − Δt Dec Ac2 2ρ ⎩ ⎭

t +Δt t +Δt K j m c m c

∑A j∈Lc

2 c

2ρ

−

∑ ρ g ΔH t k

k

k∈Lc

Mass ( N − 1)m 2t +Δt + m1t +Δt = m ct +Δt If the density distribution is known from the past time step, the three momentum equations and the mass conservation equation are four non linear equations in the four unknowns

m 1t +Δt , m 2t +Δt , m ct +Δt and ΔPct +Δt which can be solved by the Newton-Raphson Technique. Note, the equations are mostly linear, with the nonlinearity appearing in the friction and forms loss terms. Assuming that the friction factor (f) varies slowly with mass flow rate, a Newton-Raphson linearization of the system yields the iteration equation for the new time values Loop 1 Li ⎧ m 1p +1 − m 1t ⎫ p +1 ⎨ ⎬ = −ΔPc − Δ t ⎭ 1i ⎩

∑A i∈L1

f1ip L1i 1 (2m 1p +1 − m 1p ) | m 1p | − 2 2ρ 1i A1i

∑ De i∈L1

K j (2m 1p +1 − m 1p ) | m 1p | − 2 2ρ 1j

∑A j∈L1

∑ ρ g ΔΗ t k

k

+ ΔPpt1

k∈L1

27

Loop 2 Li ⎧ m 2p +1 − m 2t ⎫ p +1 ⎨ ⎬ = −ΔPc − Δt ⎭ 2i ⎩

∑A i∈L2

f 2pi L2i 1 (2m 2p +1 − m 2p ) | m 2p | − 2 2ρ 2 i A2 i

∑ De i∈L2

K j (2m 2p +1 − m 2p ) | m 2p | ρ kt g ΔΗ k + ΔPpt 2 − 2 2ρ 2j k ∈L

∑A j∈L2

∑ 2

Core p +1 p p Lc ⎧ m cp +1 − m ct ⎫ f c p Lc 1 (2m c − m c ) m c p +1 − ⎨ ⎬ = ΔPc − Δt Ac ⎩ Dec Ac2 2ρ ⎭

p +1 p p K j (2m c − m c ) m c

∑A j∈Lc

2 c

2ρ

−

∑ ρ g ΔH t k

k

k∈Lc

Mass ( N − 1)m 2p +1 + m1p +1 = m cp +1 The linearized momentum and mass conservation equations have the form

a1m 1P +1 + ΔPcP +1 = b1 a2 m 2P +1 + ΔPcP +1 = b2 ac m cP +1 − ΔPcP +1 = bc m1P +1 + ( N − 1)m 2P +1 − m cP +1 = 0 which may be written in matrix form 0 0 1 ⎞ ⎛ m 1 ⎞ ⎛ b1 ⎞ ⎛ a1 ⎟ ⎜ ⎟ ⎜ ⎟⎜ a2 0 1 ⎟ ⎜ m 2 ⎟ ⎜ b2 ⎟ ⎜0 = ⎜0 a3 −1⎟ ⎜ m c ⎟ ⎜ bc ⎟ 0 ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎝ 1 (n − 1) −1 0 ⎠ ⎝ ΔPc ⎠ ⎝ 0 ⎠

Due to the sparseness of the coefficient matrix, Cramer’s rule can be used to solve directly for m 1p +1 as

m 1p +1 =

b1a2 + ( N − 1)b1ac − ( N − 1)ac b2 + a2bc a1a2 + ( N − 1)a1ac + a2 ac

The remaining new iterate values can then be easily obtained from a) ΔPc p +1 = b1 − a1m 1p +1

28

p +1 2

b) m

b2 − ΔPc p +1 = a2

c) m cp +1 = m 1p +1 + ( N − 1)m 2p +1 Iteration to convergence yields the new time mass flow rates around the loop.

Fluid Temperature Distribution Proceeding to the next time step requires solution of the energy equation for the temperature distribution around the loop. We have seen that the Internal Energy Equation can be written as

Ax ρ

∂u ∂u + ρ vAx = qw′′ ( z ) Pw ( z ) ∂t ∂z

We can either solve for internal energy, then use a state equation to extract temperature, or solve for temperature directly through the transformation

∂u ∂u ∂T = ∂ξ N ∂T ∂ξ

ξ = z, t

Cp

such that the energy equation may be written

Ax ρ C p

∂T ∂T + ρ vAx C p = qw′′ ( z ) Pw ( z ) ∂t ∂z

Since the wall heat transfer rate is given by

qw′′ = hc (Tw − T ) and the fluid is assumed to be single phase at all times, temperature is the more convenient variable. Previous analysis of time and spatial discretizations for this equation on the mesh below,

showed that for velocity positive, the discretization 29

⎧ T t +Δt − Tkt ⎫ t +Δt t +Δt t* M kCp ⎨ k ⎬ + m k C p {Tk +1/2 − Tk −1/2 } = qk ⎩ Δt ⎭

Tk +1/2 = Tk Tk −1/2 = Tk −1 qk = hc Ak (Twk − Tk ) yielded a solution that was unconditionally stable and produced physical (non negative) results. In a transient, it is possible to have flow reversal such that velocity (mass flow) rate can be negative. An extension of the above to account for negative velocities gives For

m k > 0 Tk +1/2 = Tk Tk −1/2 = Tk −1

m k < 0 Tk +1/2 = Tk +1 Tk −1/2 = Tk

To avoid “if” statements when programming the above, the convected values at the boundaries can be expressed analytically as | m | 1 1 Tk +1/2 = (Tk + Tk +1 ) + k (Tk − Tk +1 ) 2 m k 2 | m | 1 1 Tk −1/2 = (Tk + Tk −1 ) + k (Tk −1 − Tk ) 2 m k 2 Giving for the spatially discretized energy equation M kCp

dTk m k + 2 dt

⎧ m k ⎫ m k ⎨1 − ⎬ C pTk +1 + | m k | C pTk − 2 ⎩ m k ⎭

⎧ | m k | ⎫ ⎨1 + ⎬ C pTk −1 = qk m k ⎭ ⎩

or assuming an implicit time discretization ⎧ T t +Δt − Tkt ⎫ m k M kCp ⎨ k ⎬+ ⎩ Δt ⎭ 2

⎧ m k ⎫ m k ⎧ | m k | ⎫ t +Δt t +Δt t +Δt t* ⎨1 − ⎬ C pTk +1 + | m k | C pTk − ⎨1 + ⎬ C pTk −1 = qk   2 m m k ⎭ k ⎭ ⎩ ⎩

30

Manifold Analysis (Upper Plenum)

The previous discretization holds only for one-dimensional, single inlet, single outlet nodes. The same approach can be extended to manifolds. Consider the example of the upper plenum illustrated above. Mass conservation gives i 2 (Holds even during flow reversal) m c = m 1 + m i 2 = ( N − 1)m m 2

and the energy equation applied to the manifold (assuming no heat transfer) gives

M jcp

dT j dt

 T = m c C pT j −1/ 2 − m 1C pT j +1/2 − mC p k −1/ 2

T j −1/2 =

| m | 1 1 T j + T j −1} + c {T j −1 − T j } { m c 2 2

T j +1/ 2 =

| m | 1 1 T j + T j +1} + 1 {T j − T j +1} { m 1 2 2

Tk −1/2 =

| m | 1 1 T j + Tk } + 2 {T j − Tk } { m 2 2 2

Substituting above gives M jC p

dT j dt

=

m c ⎧ | m c | ⎫ m 1 ⎧ | m 1 | ⎫ m 2 ⎧ | m 2 | ⎫ 1 i 2 | C T ⎨1 + ⎬ C pT j −1 − ⎨1 − ⎬ C pT j +1 − ⎨1 −  ⎬ C pTk − | m c | + | m 1 | + | m p j m c ⎭ m 1 ⎭ m 2 ⎭ 2 ⎩ 2 ⎩ 2 ⎩ 2

{

}

31

Manifold Analysis (Lower Plenum)

A similar treatment of the lower plenum gives

M jC p

dT j dt

= m 2C pTk +1/2 + m 1C pT j −1/2 − m c C pT j +1/ 2

Tk +1/2 =

| m | 1 1 T j + Tk } + 2 {Tk − T j } { m 2 2 2

T j −1/2 =

| m | 1 1 T j + T j −1} + 1 {T j −1 − T j } { m 1 2 2

T j +1/ 2 =

| m | 1 1 T j + T j +1} + c {T j − T j +1} { m c 2 2

M jC p

dT j dt

=

m c m 1 ⎧ | m 1 | ⎫ ⎨1 + ⎬ C pT j −1 − m1 ⎭ 2 ⎩ 2

⎧ | m c | ⎫ m 2 ⎨1 − ⎬ C pT j +1 + mc ⎭ 2 ⎩

⎧ | m 2 | ⎫ 1 ⎨1 + ⎬ C pTk − | m c | + | m 1 | + | m 2 | C pT j  m 2 2 ⎭ ⎩

{

}

32

Numerical Solution of the Internal Energy Equation

Consider the simple flow network illustrated above. The numbering scheme chosen for the temperature solution will affect the matrix structure (non zero distribution). Choose a numbering scheme such that the 1-D segments are numbered first, and the manifolds are numbered last. This numbering produces the matrix structure illustrated below.

While the matrix is 11 x 11, it is primarily tridiagonal with nonzero columns and rows corresponding to connections at the manifolds. We can think of the number of manifolds as a measure of the “complexity” of the problem. Increasing the number of nodes in the 1-D segments increases the length of the tridiagonal segment only. To increase the number of off-tridiagonal rows and columns requires the addition of additional manifolds.

33

Matrix Solution by Block Elimination The matrix equation for the temperature solution has the form

The matrix F is a square matrix with the same number of rows and columns as the number of manifolds, or the “complexity”. The matrix T is tridiagonal. Rewrite the matrix equation as the system of matrix equations a) Tx1 + Rx2 = S1 b) Bx1 + Fx2 = S 2 Solve Equation a) for x1

x1 + T −1 Rx2 = T −1S1 c) x1 = T −1S1 − T −1 Rx2 Construction of T −1S1 involves application of the Thomas Algorithm for tridiagonal systems to S1 . Construction of T −1 R involves application of the Thomas Algorithm to each column in R. Substitute

x1 = T −1S − T −1 Rx2 into b)

B[T −1S − T −1 Rx2 ] + Fx2 = S2 [ F − BT −1 R]x2 = S 2 − BT −1S1 which is a 2 x 2 linear system for x2 . Given x2 , x1 may be obtained directly from Equation c). 34