Sample Problems > Applications > Fluids > landfill_gas

{ LANDFILL_GAS_FLOW.PDE

This script solves 2D flow of ideal gas through a porous medium combined of

2 contiguous layers of distinct permeability in elliptic coordinates. The

outermost of the two layers generates gas that flows inwards and is

collected by a pipe in the centre. When the cross-sections are ellipses of

a common focal length f, there are exact solutions available that are used

for verification. Note: the domain appears almost circular, but if the

coordinates are checked carefully, all boundaries are ellipses.

Includes usage of:

- redefinition of Cartesian coordinates to custom curvilinear system

- definition of inverse hyperbolic functions

- discontinuous material properties

- definition of features via a repeat loop

- file export via a repeat loop

Written at:

Department of Mathematics, Thompson Rivers University (British

Columbia, Canada) by Damian Halvorsen and Yana Nec.

Reference for the geometry and all formulae: appendix B in

DOI: https://doi.org/10.1016/j.pce.2018.10.003 ;

http://faculty.tru.ca/ynec/index_papers.html

}

Title 'Planar flow, elliptic geometry, 2 laminae'

COORDINATES cartesian2

VARIABLES P

DEFINITIONS

C2K = 273.15 { 0 deg C in Kelvin }

g = 9.8067 { gravity constant m/s^2 }

ps = 101325 { standard atmospheric pressure at sea level Pa }

Ts = C2K+15 { standard atmospheric temperature at sea level K }

Tc = C2K+20 { standard chemical temperature K }

Latm = 0.0065 { standard atmospheric cooling rate deg/m }

R_air = 286.9 { air gas constant p = rho R T }

Ro = 8.3145 { universal gas constant, specific R = Ro/M }

i2m_p = ps/406.8 { inWc to Pa conversion factor }

elev = 0; { surface elevation m }

p_bar = ps { barometer pressure Pa }

p_atm = p_bar*(1-Latm*elev/Ts)^(g/(R_air*Latm)) { atmospheric pressure at given elevation}

pB = p_atm-0.5*i2m_p { pressure under the cover (if applicable), converted from inWc }

p_out = p_atm-15*i2m_p { outlet pressure, converted from inWc }

d_P = 0.5*0.3048 { pipe diameter m }

r_P = d_P/2 {pipe radius}

{ pack, waste and cover laminae thickness m }

h = array(1,8,3)

r_A = r_P+h[1] {gravel pack radius}

r_B = r_A+h[2] {landfill cavity radius}

{ porosity }

p_eps = array(0.6,0.4,0.7)

{ effective packing sphere diameter to represent grain size m }

sphere_diam = array(0.05,0.1,0.01) {k_b>k_a}

{ tortuosity }

tau = 100

k = array for i (1 by 1 to 3) : p_eps[i]^3*sphere_diam[i]^2*(72*tau*(1-p_eps[i])^2)^(-1)

{ temperature deg C }

Tlfg = 15+C2K

{ hole diameter m }

dh = 3/8*0.0254

{ # of holes in each perforated cross-section }

nh = 2

dl = nh*(dh/2)^2/d_P

{molar weights}

M_CH4 = 0.016044

M_CO2 = 0.04401

M_O2 = 0.0319988

M_N2 = 0.0280134

Mw = array(M_CH4,M_CO2,M_O2,M_N2)

{ gas composition CH4, CO2 and O2, balance N2 }

Xlfg0 = array(0.5,0.4,0.01);

Xlfg = array(Xlfg0[1],Xlfg0[2],Xlfg0[3],1-sum(i,1,3,Xlfg0[i]))

Rlfg = Ro/sum(i,1,4,Xlfg[i]*Mw[i])

{ Sutherland's formula }

{ base values for Sutherland's formula }

To_CH4 = 273.15; s_CH4 = 197.8; muo_CH4 = 12.01*10^(-6);

To_CO2 = 293.15; s_CO2 = 240; muo_CO2 = 14.8*10^(-6);

!To_CO2 = 273.15; s_CO2 = 222.2; muo_CO2 = 13.7*10^(-6);

To_O2 = 292.25; s_O2 = 127; muo_O2 = 20.18*10^(-6);

To_N2 = 300.55; s_N2 = 111; muo_N2 = 17.81*10^(-6);

mu_CH4 = muo_CH4*(Tlfg/To_CH4)^1.5*(To_CH4+s_CH4)/(Tlfg+s_CH4);

mu_CO2 = muo_CO2*(Tlfg/To_CO2)^1.5*(To_CO2+s_CO2)/(Tlfg+s_CO2);

mu_O2 = muo_O2*(Tlfg/To_O2)^1.5*(To_O2+s_O2)/(Tlfg+s_O2);

mu_N2 = muo_N2*(Tlfg/To_N2)^1.5*(To_N2+s_N2)/(Tlfg+s_N2);

{ exponential formula }

mu = array(mu_CH4, mu_CO2, mu_O2, mu_N2)

a = 3/8 { empirical constant from Thomas A. Davidson's paper - 3/8 or 1/3

a = 3/8 appears to work negligibly better }

Ylfg = sum(i,1,4,sqrt(Mw[i])*Xlfg[i])

fl = 2^a/Ylfg^2*sum(i,1,4,sum(j,1,4,

Xlfg[i]*Xlfg[j]*Mw[i]^((a+1)/2)*Mw[j]^((a+1)/2)/

(sqrt(mu[i]*mu[j])*(Mw[i]+Mw[j])^a)))

mu_lfg = 1/fl

{ number of perforated cross-sections }

n = 28

{ generation rate m^3/hr }

C_init = 375

C = C_init*ps/(3600*Rlfg*Tc*pi*(r_B^2-r_A^2)*dl*n)

RT = Rlfg*Tlfg

muRT = mu_lfg*RT

muRTC = muRT*C

ang_jump = pi/8

Px ! constant, to be used as a normalisation pressure value

Px2 = Px^2

f = 0.5*r_P ! focal length, common to all ellipses

fb = f^2*muRTC/(4*Px2*k[2]) ! combined constants

fc = f^2*muRT*Cp/(4*Px2*Kp)

muRTCpx = muRT*Cp/Px2

grav = 0*vector(0,-g,0) ! toggle the factor to turn gravity on and off

! define inverse hyperbolic functions not available as built in functions

atanh(q) = 0.5*ln((1+q)/(1-q))

acosh(q) = ln(q+sqrt(q^2 -1))

asinh(q) = ln(q+sqrt(q^2 +1))

{ define the elliptic coordinates: xi runs along the hyperbolae (normal

coordinate) and phi is the elliptic arc length (tangential coordinate)

Note: radial arc length = angle from horizontal,

elliptic arc length != angle from horizontal except 0, pi

}

xyf = f^2+y^2-x^2

tanhxi = if(x=0) then sqrt(y^2/xyf)

else sqrt(0.5*(-xyf+sqrt(xyf^2+(2*x*y)^2))/x^2)

xi = atanh(tanhxi)

phi = if(x=0) then pi/2*sign(y)

else if(x>0) then arctan(y/(x*tanhxi))

else arctan(y/(x*tanhxi))+pi

! common expressions needed in exact solutions

xip = acosh(r_P/f)

xia = acosh(r_A/f)

xib = acosh(r_B/f)

xiPB = xip-xib

xiBA = xib-xia

xiAP = xia-xip

r_outx = f*cosh(xib)

r_outy = f*sinh(xib)

r_inx = f*cosh(xia)

r_iny = f*sinh(xia)

r_px = f*cosh(xip)

r_py = f*sinh(xip)

! coefficients needed in exact solutions

! block/unblock one type of boundary condition

! pressure boundary condition

Px = pB

rhs = -fb*(1+k[2]/(k[1]-k[2])*cosh(2*xiBA)+

k[1]/(k[1]-k[2])*(sinh(2*xiPB)+sinh(2*xiBA)*cosh(2*xiAP))/sinh(2*xiAP))

theta_P = 0

theta_B = rhs-theta_P*k[1]/k[2]*sinh(2*xiBA)/sinh(2*xiAP)

a0a = k[2]/(k[2]*xiAP+k[1]*xiBA)*((pB/Px)^2-(p_out/Px)^2+

fb*(cosh(2*xib)-cosh(2*xia))-2*fb*sinh(2*xia)*xiBA)

a0b = 1/(k[2]*xiAP+k[1]*xiBA)*(((pB/Px)^2-(p_out/Px)^2+

fb*(cosh(2*xib)-cosh(2*xia)))*k[1]+2*k[2]*fb*sinh(2*xia)*xiAP)

b0a = (p_out/Px)^2-a0a*xip

b0b = (pB/Px)^2-a0b*xib+fb*cosh(2*xib)

a1a = 1/(sinh(2*xiPB))*(theta_P*cosh(2*xib)-theta_B*k[2]/k[1]*cosh(2*xip)+

fb*k[2]/k[1]*cosh(2*xip)*(cosh(2*xiBA)-1))

b1b = 1/(sinh(2*xiPB))*(-theta_P*k[1]/k[2]*sinh(2*xib)+theta_B*sinh(2*xip)+

fb*(sinh(2*xip)-sinh(2*xib)*cosh(2*xiAP)))

a1b = 1/(sinh(2*xiPB))*(theta_P*k[1]/k[2]*cosh(2*xib)-theta_B*cosh(2*xip)+

fb*(cosh(2*xib)*cosh(2*xiAP)-cosh(2*xip)))

b1a = 1/(sinh(2*xiPB))*(-theta_P*sinh(2*xib)+theta_B*k[2]/k[1]*sinh(2*xip)+

fb*k[2]/k[1]*sinh(2*xip)*(1-cosh(2*xiBA)))

P_ell = sqrt(a0*xi+b0-fc*cosh(2*xi)+

cos(2*phi)*(a1*sinh(2*xi)+b1*cosh(2*xi)-fc))

MATERIALS

'ACmat' :

Kp = k[2]

Cp = C

a0 = a0b

b0 = b0b

a1 = a1b

b1 = b1b

'ABmat' :

Kp = k[1]

Cp = 0

a0 = a0a

b0 = b0a

a1 = a1a

b1 = b1a

INITIAL VALUES

P = P_atm/Px

EQUATIONS

div(Kp*p*(grad(p)-p/RT*grav)) = -muRTCpx

BOUNDARIES

REGION 'AC'

use material 'ACmat'

start (r_outx,0)

value(P) = sqrt((pB/Px)^2 +theta_B*cos(2*phi)) ! pressure boundary condition

arc(center=0,0) to (0,r_outy) to (-r_outx,0)

to (0,-r_outy) to (r_outx,0)

start (r_px,0)

value(P) = sqrt((P_out/Px)^2 + theta_P*cos(2*phi))

arc(center=0,0) to (0,r_py) to (-r_px,0)

to (0,-r_py) to (r_px,0)

REGION 'AB'

use material 'ABmat'

start (r_inx,0)

arc(center=0,0) to (0,r_iny) to (-r_inx,0)

to (0,-r_iny) to (r_inx,0)

start (r_px,0)

arc(center=0,0) to (0,r_py) to (-r_px,0)

to (0,-r_py) to (r_px,0)

repeat th = 0 by ang_jump to 3/8*pi

Feature 'Ray_'+$th

start(r_px*cos(th),r_py*sin(th)) line to (r_outx*cos(th), r_outy*sin(th))

endrepeat

PLOTS

contour(P) painted zoom(-r_px,-r_px,2*r_px)

repeat th = 0 by ang_jump to 3/8*pi

Elevation (P, P_ell) on 'Ray_'+$th as 'P on Ray_'+$th

export points=30 file = 'qpe2L_'+$round(th*10000)+'.tbl'

Elevation (P-P_ell) on 'Ray_'+$th as 'P error on Ray_'+$th

endrepeat

END

landfill_gas_flow

landfill_gas_flow