landfill_gas_flow

<< Click to Display Table of Contents >>

Navigation:  Sample Problems > Applications > Fluids >

landfill_gas_flow

Previous pageReturn to chapter overviewNext page

{ 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

 

 Kp

 Cp

 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)))

 a0

 b0

 a1

 b1

 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