fire tube

a.Rmd

# Water Heater

For backpacking weight is an important constraint. Often boiling water is
enough for food and drink -- the food can be precooked and/or dehydrated etc.
such that it doesn't need to simmer. Therefore a stove and a pot could be
unnecessary, if a tube containing burning fuel is put inside a water bottle.

Methanol stoves where the pot has a shield ("caldera cone", jetboil etc.)

The fuel is tentatively methanol. Using compressed gas (propane isobutane etc.)
or liquid fuel under pressure would allow control of the fuel rate, but those
are harder to make. Esbit hasn't been rejected.

People carry phones and possibly chargers for other reasons. Therefore a fan
rather than natural draft is under consideration. The fan reduces thermal
resistance on the tube side, and there may be the possibility of an optimal fan
speed profile. The outlet temperature could be measured, but the relation between
the two is unclear:

        increasing fan speed increases the heat transfer coefficient

        increasing fan speed increases the mass transfer coefficient for fuel evaporation

Using empirical correlations below we see that thermal resistance between the
gas and the internal wall is the highest, though natural convection between the water
and the outside of the tube is not negligible.

Ideally there would be very little excess air. The increase in mass
transfer coefficient following an increase in fan speed should result in no
change in the excess air. This part is unknown. In one sense it is a
humidification (of methanol), and the temperature depends on the reaction too.

CHNOSZ didn't work. Perhaps because it is used in geochemistry and biochemistry.

NASA-CEA is from http://akrmys.com/links/softw.html.en For now I just read the
output file. I changed the EFMT subroutine. In the future it could be
worthwhile to adapt the code to compile as a shared object and then do a
dyn.load for R?  The fortran seems to open thermo.lib and trans.lib files not
read-only ie. it's not reentrant?

Rate: http://akrmys.com/KUCRS/distrib.htm.en unless the methanol-specific paper is of interest?

TODO:
        methanol: diffusivity in air(?), hvap, hcombust(CHNOSZ is unreliable)

        qs = (NaMaCpa + NbMbCpb)/ (1- exp(- (NaMaCpa+NbMbCpb)/h ))

        but 

        Using R Treybal Mass Transfer Operations

        chilton-colburn (3.57)


        We have the case "steady-state diffusion of A through nondiffusing B"
        Nb = 0
        Na / (Na + Nb) = 1 (2.27)
        
        (3.1)
        Na = F ln [ ( 1 - yBulk) / (1 - ySurf) ]


```{r}
library(glue)
library(tidyverse)

cea2 <- function(T.K=c(400, 500), P=c(0.9, 1)) {
        inp <- tempfile(fileext='.inp')
        ceainp <- strsplit(inp, '\\.') %>% `[[`(1) %>% head(-1) %>% paste(collapse='.')
        out <- glue("{ceainp}.out")
        # TODO  reactants input
        # this means tp might not be the right -- want transport properties
        # before equilibrium too, but SUBROUTINE TRANP is not that easy?
        system(glue("
          (
           cat <<-CEAINP
                problem
                  tp
                  p(bar) = {paste(P, collapse=' ')}
                  t(k)   = {paste(T.K, collapse=' ')}
                reactant
                  name = CH3OH  moles = 0.11
                  name = O2   moles = 0.21
                  name = N2   moles = 0.79
                end
                outp transport
                end
CEAINP
            ) > {inp}
          fcea2m.sh {ceainp}
        "))

        outct <- readLines(out)

        system(glue("rm -f {inp} {out}"))

        bp <- match( c(" THERMODYNAMIC PROPERTIES", " MOLE FRACTIONS", "    PRODUCTS WHICH WERE CONSIDERED BUT WHOSE MOLE FRACTIONS", "   CONDUCTIVITY IN UNITS OF MILLIWATTS/(CM)(K)", "  WITH EQUILIBRIUM REACTIONS", "  WITH FROZEN REACTIONS", " MOLE FRACTIONS"), outct)

        ceaBlock <- function(ix, lw=16) {
                x <- outct[ix]
                xD <- str_sub(x, lw) %>% read_delim(' ', col_names=F, trim_ws=T) %>% t()

                xN <- str_sub(x, 2, lw)
                colnames(xD) <- xN[ xN != ''] %>%
                        strsplit(',') %>% map_chr(`[`(1)) %>%
                        strsplit(' ') %>% map_chr(`[`(1))
                xD
        }

        # THERMODYNAMIC PROPERTIES
        thermoP <- ceaBlock((bp[1]+2):(bp[2]-2))

        # MOLE FRACTIONS
        # TODO drop the leading * ?
        molFrac <- ceaBlock((bp[2]+2):(bp[3]-4))

        visc <- ceaBlock( (bp[4]+2) : (bp[4]+1 + bp[6]-bp[5] - 4), 17  )
        transpo.eq     <- ceaBlock( (bp[5]+2) : (bp[6]-2) )
        transpo.frozen <- ceaBlock( (bp[6]+2) : (2*bp[6]-bp[5] - 2) )

        colnames(transpo.frozen) <- paste('frozen', colnames(transpo.frozen), sep='.')
        colnames(transpo.eq)[1] <- glue('eq.{colnames(transpo.eq)[1]}')
        
        cbind( thermoP, visc, transpo.eq, transpo.frozen, molFrac) %>%
                as_tibble(.name_repair='universal') %>%
                mutate_at('CONDUCTIVITY', ~ . *  100/1000 ) # 
}
cea2 ()
```


## experiment

How do I design the wick?

Is there an experiment using water evaporation instead?

I can vary things like input flow rate, outer water temperature, and measure temperature and humidity.

That can yield a Kga for water?

We have claims about heat-mass transfer analogy holds better than the momentum-{mass,heat} analogy. Therefore stoichiometry can be assured?




```{r}
library(units)
library(ggforce)
library(CHNOSZ)
library(tidyverse)
library(deSolve)
library(rootSolve)

T.units("K")
P.units("bar")
E.units("J")


info('methanol', 'liq')
info('methanol', 'gas')
info('oxygen')
info('nitrogen')
info('CO2', 'gas')
info('H2O', 'liq')
# H2O vapor is missing,
# but it looks like we can use "liquid" water at the right pressure
# and get vapor densities, H is good at 373K P=1 -> P=1.1

# but we also have methanol evaporating
# use the meox1 or meox2
# how to write the equation when T cannot be solved explicitly
# how to determine meox2 vs meox1? There is CO +1/2O2 <=> CO2
# at a specific temperature.   coxox$out$logK = log10( 
#  run.wjd should be fine for that
run.wjd.example <-
 run.wjd( c('carbon dioxide', 'carbon monoxide', 'oxygen') %>% map(info),
        Y=c(0.1, 0.2, 0.1),
        T= 400, P=0.1)



# https://www.engineeringtoolbox.com/water-properties-temperature-equilibrium-pressure-d_2099.html
water <- read_csv('data/liquidWater.csv')


# negative coefficient means reactant
meox2 <-  subcrt(  c(  info('methanol', 'gas') , 'oxygen', 'carbon dioxide', 'water'), c(-1, -1.5, 1, 2 ) , P=0.2)
meox1 <-  subcrt(  c(  info('methanol', 'gas') , 'oxygen', 'carbon monoxide', 'water'), c(-1, -1, 1, 2 ) , P=0.2)

coxox <- subcrt(  c('carbon dioxide', 'carbon monoxide', 'oxygen'), c( 1, -1, -0.5), P=0.1 )
# water-gas shift?

# I want this except the composition changes
air <- within(read_csv('data/engToolboxSIair.csv'), {
        T <- T # * K
        Cp <- Cp # * kJ / kg / K
        Cv <- Cv # * kJ / kg / K
        mu <- mu * 10^-5 # * kg / m / s
        nu <- nu # * m^2 / s
        rho<- rho # * kg / m^3
        alpha<-alpha # *m^2/s
})

cps <- data.frame(
           name=c('argon', 'methanol', 'carbon monoxide', 'carbon dioxide', 'water', 'oxygen', 'nitrogen' ) ,
           mw = c(39.948, 32.042, 28.010, 44.010, 18.015, 31.999, 28.014), # kg / kmol...
           a0 = c(2.5, 4.714, 3.912, 3.259, 4.395, 3.630, 3.539),
           a1 = c(0, -6.986, -3.913, 1.356, -4.186,  -1.794,  -0.261) / 10^3,
           a2 = c(0, 4.211, 1.182, 1.502, 1.405,  0.658,  0.007) / 10^5,
           a3 = c(0, -4.443, -1.302, -2.374, -1.564,  -0.601,  0.157)/ 10^8,
           a4 = c(0, 1.535, 0.515, 1.056, 0.632,  0.179,  -0.099) / 10^11 )


airfrac <- c(nitrogen=78.084 , oxygen=20.9476 , argon=0.934 , 'carbon dioxide'=0.0314 )
airfrac <- airfrac / sum(airfrac)

gasR <- 8.3145 # * J / mol / K )

cpsF <- function(n, x) with(cps[ cps$name == n, ],  gasR * ( a0 + x*(a1 + x*(a2 + x*(a3 + x*a4))) ))
cpsFderiv <- function(n, x) with(cps[ cps$name == n, ],  gasR * ( a1 + x*(2*a2 + x*(3*a3 + 4*x*a4)))) 

mw <- function(n) cps[ cps$name == n, 'mw' ]


# kJ/ kg.K
# (.x / mw(.y)) is [g of species .y]   in ( sum( airfrac ) == 1 [mol of air]  )
cpsAir <-  . %>%
        map(.,  function(temp)  imap( airfrac, ~  cpsF(.y, temp) * .x/mw(.y) ) %>%
                        reduce(`+`)) %>%
        do.call(c, .)


# cross_df...
air2 <- air %>%
        mutate( cpsmix = cpsAir(.data$T),
                cpRat = (cpsmix - Cp) / cpsmix)


# tube properties
astmB88 <- within( read_csv('data/astmB88cu.csv'), {
        weight <- weight * drop_units(make_units(kg/m) / make_units( lb / ft ))
})
astmB88[ , 3:7]  <- lapply( astmB88[ , 3:7], function(x) x * drop_units( ud_units[['in']] / ud_units$m ) )



# lennard jones viscosity
ljv <- merge( read_csv('data/lennard-jones-visco.csv', comment='#'), cps, by='name')


neufeld1972.coefs <- list(A=1.16145, B=0.14874, C=0.52487, D=0.77320, E=2.16178, F=2.43787)
```
```{r}

neufeld1972.omega <- function(x) with(neufeld1972.coefs, A*x^B + C*exp(-D*x) + E*exp(-F*x) )

# atomic/molecular gas viscosity
mu.chapmanEnskog <- function(n, T) {
        ljvi <- ljv%>%filter(name==n)
        26.69e-7 * sqrt( mw(n) * T )  /
                        ljvi$sigma^2 / neufeld1972.omega( T / ljvi$epsilon.by.k )
                        
}

# gas mixture viscosity
mu.reichenberg1979 <- function(ys=airfrac, T=298,
                            mu_i = imap_dbl( ys, ~ mu.chapmanEnskog( .y , T) )) {
        
        ljvi <- ljv %>% filter(name %in% names(ys))

        Mi <- ljvi$mw # g/mol used on pg 9.20


        # 9-4.17
        mu_ri <- with(ljvi, 52.46 * mu.debye * Pc.bar / Tc.K^2)


        T_ri <- T / ljvi$Tc.K
        # 9-5.5
        F_ri <- (T_ri ^ 3.5 + (10*mu_ri)^7 ) / T_ri^3.5 / (1 + (10*mu_ri)^7)


        U_i <- (1+0.36*T_ri*(T_ri- 1))^(1/6)* F_ri / sqrt(T_ri)
        C_i <- Mi^(1/4) / (mu_i * U_i)^0.5

        CCterm <- outer(C_i, C_i , FUN='+')
        T_rij <- T / with(ljvi, outer(Tc.K, Tc.K))

        mu_rij <- sqrt( mu_ri %o% mu_ri )
        F_rij <- (T_rij ^ 3.5 + (10*mu_rij)^7 ) / T_rij^3.5 / (1 + (10*mu_rij)^7)

        MMterm <- sqrt(   outer(Mi, Mi, FUN=function(i, j) (i*j) / 32 / (i+j)) )


        H_ij <- MMterm * CCterm * (1+0.36*T_rij*(T_rij-1))^(1/6) * F_rij / sqrt(T_rij)

        HMM <- H_ij * outer( Mi, Mi, FUN= function(i, j) 3 + 2*i/j)
        diag(HMM) <- 0

        stopifnot(  abs(sum(ys) - 1) < 1e-8 ) # or just fix it?

        Hterm <- HMM %*% ys

        K_i <- ys * mu_i / (ys + mu_i * Hterm )

        diag(H_ij) <- 0

        mu_m <- t(1 + 2*H_ij%*%K_i + (H_ij%*%K_i) * (H_ij%*%K_i)) %*% K_i

        c(mu_m)
}

# conductivity
# 10-3.14 Poling etc.
# it is quite bad in comparison to tabulated data.
lambda.chung88 <- function(sp='oxygen', T=298)  { 
        M <- mw(sp) / 1000 # kg/mol
        mu <- mu.chapmanEnskog(sp, T)
        Cv <- cpsF(sp, T) - gasR # ideal gas

        ljvi <- ljv[ sp == ljv$name, ]

        Tr <- T / ljvi$Tc
        omega <- ljvi$omega

        alpha <- Cv/ gasR - 3/2
        beta <- 0.7862 - 0.7109*omega + 1.3168*omega^2
        Z <- 2 + 10.5 * Tr^2

        psi <- 1 + alpha * ((0.215 + 0.28288*alpha - 1.061*beta + 0.26665*Z) /
                            (0.6366 + beta * Z + 1.061*alpha*beta))

        3.75* psi * mu * Cv / (alpha+3/2) / M
}

cmw <- read_csv('data/conductivitymW.csv') %>%
        pivot_longer(matches('.00'), names_to='T', values_to='lambda', values_drop_na=T) %>%
        modify_at('T', as.numeric) %>%
        modify_at('extrapp0', ~ c('Y'=T, 'N'=F)[.] ) %>%
        modify_at('lambda', ~ . / 1000) #  * with(ud_units, W / m / K))


cmw.approx <- function(name='oxygen', T=298) {
        cmi <- cmw[ cmw$name==name, ]
        spline( cmi$T, cmi$lambda, xout=T)$y
}

# 10-6.2
lambda.masonSaxena1958 <- function(ys=airfrac, T=298, epsilon=1) {

        ljvi <- ljv %>% filter(name %in% names(ys) )
        M <- ljvi$mw
        Trs <- T/ljvi$Tc.K

        # 10-3.9
        Gammas <- 210 * (ljvi$Tc.K * M^3 / ljvi$Pc.bar^4)^(1/6)

        lambdas <- map_dbl(names(ys), ~ cmw.approx(., T))


        Aijyj <- map_dbl( seq_along(ys), function(i)
                         sum( map_dbl(seq_along(ys), function(j)  {

              #10-6.5
              lamrat <- Gammas[j] * (exp(0.0464*Trs[i]) - exp(-0.2412*Trs[i])) /
                      ( Gammas[i] * (exp(0.0464*Trs[j]) - exp(-0.2412*Trs[j]))  )
              #lamrat <- lambdas[i] / lambdas[j]

              # 10-6.3
              #lamrat <- mu[i]* M[j] / (mu[j] * M[i])

              # 10-6.2
              ys[j] * epsilon * (1 + lamrat^0.5 * (M[i]/M[j])^0.25)^2 /
                      ( 8 * (1 + M[i]/M[j]))^0.5
        } )))

        sum( ys * lambdas / Aijyj )
}

# if we just take the Aijyj=1,
lambda.lerp <- function(ys=airfrac, T=298) {
        lambdas <- map_dbl(names(ys), ~ cmw.approx(., T))
        sum(ys * lambdas)
}

# reuse the reichenberg1979 as suggested by 10-6.4
# but this one is too low by about 10x
lambda.reichenberg1979 <- function(ys=airfrac, T) 10*mu.reichenberg1979(ys, T,
                                                          mu_i = imap_dbl(ys,~ cmw.approx(.y, T)) ) 



# compare lambda.masonSaxena1958 with  cmw for air
# it is not too wrong
cmw.cmp <- cmw %>% filter(name == 'air') %>% within(., {
        lambda.measured <- lambda
        lambda.mason <- Vectorize(function(x) lambda.masonSaxena1958(T=x))(T)
        lambda.lerp <- Vectorize(function(x) lambda.lerp(T=x))(T)
        lambda.reichenberg.times.10 <- Vectorize(function(x) lambda.reichenberg1979(T=x))(T)
}) %>% pivot_longer(matches('lambda.'), names_to='method', names_prefix='lambda.')


h.nconvect.tab <- rbind(
        tibble( orient='v', ymax = c(1e4, 1e9, Inf),
                            ymin = c(-Inf, 1e4, 1e9),
                                a=c(1.36, 0.59, 0.13),
                                m=c(1/5, 1/4, 1/3)),
        tibble( orient='h',
                ymax = c(1e-5, 1e-3,1, 1e4,1e9, Inf),
                ymin = c(-Inf, 1e-5, 1e-3,1, 1e4,1e9),
                a=c(0.49, 0.71,1.09,1.09,0.53,0.13),
                m=c(0, 1/25,0.1,1/5,1/4,1/3)))

# XXX boiling?
h.nconvect.water <- function(Tw, Tinf, L, orientation='v') {

        if (Tw > 373) warning('h.nconvect.water needs a boiling correlation')
        Tf <- (Tw + Tinf)/2
        Pr <- with(water, approxfun(T, Pr))(Tf)
        k <- with(water, approxfun(T, k))(Tf)
        nu <- with(water, approxfun(T, mu/rho))(Tf)
        Beta <- with(water, splinefun(T, 1/log(rho)))(Tf, deriv=1)
        Gr <- 9.8 * Beta * (Tw - Tinf) * L^3 / nu^2

        Y <- Gr * Pr

        am <- h.nconvect.tab %>% filter(ymax > Y & ymin < Y & orientation == orient)

        as.numeric( k/L* am[ 1, 'a'] * Y^am[1, 'm'] )
}

```

Unclear why the units package above prevents the axis ticks from being labeled.
Also reichenberg1979 is off by ~ 8 - 10 when measured thermal conductivities
are substituted for Chapman Enskog theoretical viscosities. Anyways just use
the `lambda.masonSaxena1958`



## reaction thermodynamics

Something is wrong because T2=690K is the maximum which is far below the
tabulated Methanol-Air adiabatic flame temperature of 1949C. Maybe because we
have (wet). Maybe because property='G' should be 'H' should this leads to no
root. Hydrogen is not produced? I fixed a problem with H1 being calculated at
T2.

It should be straightforward:

    Hrxn(T0) = integral Cp nproduct from T0 to T2

```{r}
meohox <- function(phi=2, RH=0.5, T0=298, atmBar = 1) { # phi is the https://en.wikipedia.org/wiki/Equivalence_ratio
        nw_by_nair <- RH * with(water, splinefun(T, Pbar / (atmBar - Pbar)))(T0)
        dryAir <- 1.5 * phi * airfrac / airfrac['oxygen']

        ispecies <- c(info('methanol', 'gas'), map_int(c(names(airfrac), 'water', 'carbon monoxide'), info) )



        Y0 <- c(methanol=1, dryAir,  water=nw_by_nair * sum(dryAir), 'carbon monoxide'=1e-15)
        run.wjdT <- function(T2) run.wjd( ispecies, Y=Y0, T = T2 )

        # sink(file=file(open='meohox.log'), type='message')


        H1 <- sum( map2_dbl( ispecies, Y0, ~ subcrt( .x, P= atmBar * .y / sum(Y0),
                                          T=T0, property='H' )$out[[1]][[1]] * .y))
        rootF <- function(T2) {
                r <- run.wjdT(T2)
                H2 <-sum( map2_dbl( ispecies, r$Y, ~ subcrt( .x, P= atmBar * .y / sum(r$Y),
                                                             T=T2, property='H' )$out[[1]][[1]] * .y))
                H2 - H1
        }
        browser()
        rootR <- multiroot(rootF, 600)
        final <- run.wjdT(rootR$root)
        # sink()

        list( T=rootR$root, final)
}
meohox(phi=1)
```

```{r}
inputT <- crossing( phi=c(0.8, seq(0.95, 1.05, length.out=5), 1.2),
                RH = c(0.1, 0.5, 0.9)) %>%
         bind_cols(. , T2 = pmap_dbl(., ~ {
                                             print(.x)
                                             meohox(phi=.x, RH=.y)$T}) )
plot(ggplot(inputT, aes(phi, T2, col=RH, group=RH)) + geom_line())
# T2 is the initial temperature
```



```{r}
# friction factor
f.churchill <- function(Re, epsilonByD=0.01) {
        A <- (2.457 * log(((7/Re)^0.9 + 0.27*epsilonByD)^(-1)))^16
        B <- (37530/Re)^16
        2 * ( (8/Re)^12 + (A+B)^(-1.5))^(1/12)
}
# coleburn 1933 analogy
# but maybe also do a laminar flow?
# St * Pr^(2/3) = f/2
# St = h / rho Cp U
# Pr = nu / alpha = 0.7 no need to estimate I guess? alpha from lambda from chung 1988 I guess
st.coleburn1933 <- function(f, Pr=0.7) Pr^(-2/3) * f/2


ps0 <- list(tube=astmB88[1, ], Pr=0.7, kWall = 380, # W / m / k
            L = 0.3, Twater=298, ns = airfrac / 100 )

odefun <- function(x, ys, ps=ps0) with(ps, {
        T <- ys[1]
        P <- ys[2]
        # ns <- ys[ 3 : (nns +2) ]
        # when we decide to have the reaction which involves many more species
        # if following lindstedt2002

        # ns should be mol/s

        ID <- with(tube, OD - 2*wall)
        Ac <- pi * ID^2 / 4


        ljvi <- ljv[ ljv$name %in% names(ns), ]



        kgps <- sum(imap_dbl(ns, ~ mw(.y) * .x )) / 1000 # [ kg / s ]
        molps <- sum(ns) # * make_units(mol / s)

        # [kg/mol]
        kgramPerMol <- kgps / molps


        # [ J / K s]
        cPmdot <- sum(imap_dbl( ns, ~ .x * cpsF(.y, T))) # [ J / K  / s]


        # [kg/m^3] = M n / V = P / RT  * M
        rho <- P / (gasR * T) * kgramPerMol

        # [kg / m^2 /s]
        G <- kgps / Ac


        mu <- mu.reichenberg1979(ns / sum(ns), T) # [ kg / m / s]

        # [1]
        Re <- G * ID / mu

        # [1]
        epsilon <- 0.0015e-3 # in m
        f <- f.churchill(Re, epsilon / ID)


        # the usual variables
        #
        # [ J/K/kg ]
        # cP <- cPmdot / kgps
        #
        # [m/s]
        # U <- G / rho
        #
        # h <- rho * U * cP * st.coleburn1933(f)

        # but we can divide out a rho and a kgps
        h <- cPmdot / Ac * st.coleburn1933(f)


        # R * Q is [K]
        # Q [W / m], since we have d/dx

        # R [K m / W]
        # k [W / K m]
        RwallA <- with(tube, log(OD/ID)/(2*pi*kWall))


        Two <- multiroot(   function(Two) {
                        A <- 1/(h * pi * ID)
                        B <- RwallA
                        C <- 1 / (pi * tube$OD * h.nconvect.water( Two, Twater, L ))
                        Two * (A+B+C) - ((A + B) * Twater + C* T ) },
                  (T+Twater)/2 )$root

        ABC <- c(A =  1/ (h*pi*ID),
                 B=RwallA ,
                 C= 1/ (pi * tube$OD * h.nconvect.water( Two, Twater, L )))

        Tws <- with(as.list(ABC), c( A * Twater + (C+B) * T, (A+B)*Twater + C*T) / (A+B+C))
        if (Tws[2] > 373) print(Tws[2])

        # really dQ/dx
        # [W / m]
        Q <- with(as.list(ABC), (T - Twater) / (A+B+C))

        # d(Cp * T)/dt = -Q
        # Cp * dT/dt + T dCp/dT * dT/dT = -Q
        # dT/dt = -Q / [ Cp + T dCp/dT ]
        # or the Kirchoff transformation meaning Cp*T is the variable...
        # what is T dCp/dT, for water it's 1.7 while Cp is 33 (J/molK), which doesn't matter

        # (%i6) eqs : [ q * A = T1 - T2, q*B = T2 - T3, q*C = T3 - T4]$
        # (%i7) solve(eqs, [T2, T3, q]);
        #              A T4 + (C + B) T1       B T4 + A T4 + C T1       T1 - T4
        # (%o7) [[T2 = -----------------, T3 = ------------------, q = ---------]]
        #                  C + B + A               C + B + A           C + B + A
        # T1 = T
        # T2 = Twi
        # T3 = Two
        # T4 = Twater

        
        # how do T and P interact?
        # as the gas is cooled, the gas takes up less volume
        # can kinetic energy be neglected?
        list(# dT/dx, dP/dx
             c( -Q / cPmdot , # [W/m] / [W/K] => K/m
                - 2*f * G^2/rho / ID
             )

             )

})


# input:
#   % excess air
# initial temperature from energy balance (run.wjd)
oderesult <- ode( c(T=1949/2+273, P=100100), seq(0, ps0$L, length.out=100), odefun, ps0 )
```



```{r}
jpeg(width=300, height=300)
plot(ggplot(as.data.frame(oderesult) %>% mutate_at('T', ~ (T - ps0$Twater) / (T[1] - ps0$Twater) )
            %>% mutate_at('P', ~ (P - 100000) / (P[1] - 100000) )
            %>% pivot_longer(!"time"), aes(time, value, col=name))
        + geom_line()
        + ylab('m'))
dev.off()

```




```{r}

#
# dP/dL = 2*f/D * rho * u^2
# u = 4 V / pi D^2
# V = nRT / P
#
# n constant / variable with reaction
#
# reaction rate...
#
# ode( c(T=298, P=1.01, airfrac, methanol=0.2), c(0, 1)



#
# with(meox2$out, approxfun(T, H))(T)
#
#
# c(airfrac, methanol = airfrac['oxygen'] / excess )

# water
# fuel
# air (humid?)
# tubing dimensions
# fan characteristics
# fluid mechanical correlation
# excessive generalization leading to CAPE-OPEN
# use CHNOSZ

```



## unused laminar natural convection similarity solution

```{r}
# laminar natural convection from
#
# Bejan Kraus Heat transfer handbook 2003
# 
# cites Sparrow&Gregg Trans ASME 1958, 78, 435-440
# for the following:
#
# Tw-Tinf = N*x^n   (7.34)
#
# f( eta = (y/x) (Gr_x / 4)^0.25 )
# theta( eta ) = (T - Tinf) / (Tw - Tinf)
# Gr_x = g*beta * x^3 * (Tw - Tinf) / nu^2
#
# (7.35)
#
# f''' + (n+3)f f'' - 2(n+1)(f')^2 + theta = 0
# theta'' / Pr + (n+3)f theta' - 4 n f' theta = 0
#
# and boundary conditions are  0 = f'(0) = f(0) = 1-th(0) = f'(inf) = th(inf)
#
# 7.35 is solved numerically. But there is a problem
#
# q = [ - k* th1(0)/sqrt2 * (g*Beta/nu^2)^0.25 * g(x,n,N) ]|x=0 to L
# where g = 4*N^1.25*x^((5n+3)/4) / (5n + 3)
#
# N and n are unknown. They have to be solved along with the Tw.
# Probably the q should be split into two(?) sections, one from x=0 to L/2 and L/2 to L.
# then there is q(Two, Tinf, n, N) = hA (Tg - Twi) = ... (Twi - Two)
# at different locations.
#
# I don't have the original paper.  Rahman&Satish don't get 7.35 in their 20 or 22,
# but it seems like their perturbation method is not the same as (7.34)
#
# the other options are openfoam and a Nusselt equation
convect.shoot <- function(y, etas =  c(0,1), n=1, Pr=8)
        ode(  c( f2 = y[1], f1=0, f0=0, th1 = y[2], th0=1  ),
                etas,
            function(nu, fth, pars) with(c(fth,pars),
                        list(c( -(n+3)*f0*f2 + 2*(n+1)*f1^2 - th0,
                                f2,
                                f1,
                                -Pr*(n+3)*f0*th1 + 4*Pr*n*f1*th0,
                                th1))),
            list(n=n, Pr=Pr))

convect.shoot.f <- function(etaEnd, n, Pr) function(y) convect.shoot(y, etas=c(0, etaEnd), n, Pr) %>%
                                                                `[`(nrow(.), c('f1', 'th0'))

# works faster if we start with etas c(0,1) then increase the 1

check_converged <- function(xs) {
        p <- sum( abs( xs[[ length(xs) ]] - xs[[ length(xs) - 1 ]]  )) < 1e-5

        if (!p) warning("convect shooting method needs a longer etaEndSeq")

        xs[[ length(xs) ]]
}

# shooting method
convect <- function(n=0.5, Pr=8, etas = c(0, 1, 2, 3, 8, 2^10), etaEndSeq=2^(0:10))
        accumulate( etaEndSeq, ~ multiroot( convect.shoot.f(.y, n, Pr), .x)$root
                                , .init = c(1,1) ) %>%
                check_converged() %>%
                convect.shoot( etas, n, Pr )
# Nu_x = -th1(0)/ sqrt(2) * Gr_x^(1/4)
# Gr_x includes Tw-Tinf
#
# Nu_x -> Nu or otherwise

```


```{r}
plot( ggplot( air2 %>% subset( T < 1100 ), aes(T, cpRat) ) + geom_line() +
        ggtitle('[engineering toolbox air Cp - Cp]  / Cp,\n where Cp uses a 4th degree polynomial for\n Cp of N2,O2,Ar,CO2 \n with coefficients from Poling Prausnitz OConnell 2001') )



```
```{r}
plot(ggplot( cmw.cmp, aes(T,value/lambda, col= method)) + geom_point() +
     ylab(expression(lambda[method] / lambda[measured])) + ggtitle('air conductivities') )
```