Schreiner equations for CCR?

[huwporter]

I'm adding CCR functionality to my homemade deco software and I've hit a hurdle;

The tissue loading equations are straightforward enough for a stage at a constant depth. They are the same Buhlmann equations as for open circuit, you just need to work out your PN2 and PHe after taking account of the PPO2 and the proportions of inert gasses in the diluent (code is Python):

TissueN2[cntr] = TissueN2[cntr]+(PN2-TissueN2[cntr])*(1-math.pow(2,(-1*(TimeInterval/HalftimesN2[cntr]))))
TissueHe[cntr] = TissueHe[cntr]+(PHe-TissueHe[cntr])*(1-math.pow(2,(-1*(TimeInterval/HalftimesHe[cntr]))))

Equally, ascent and descent phases are simple enough on OC, using the Schreiner equations together with the fractions of inert gas inspired:

TissueN2[cntr] = (FN2*(CurrentPressure-WaterVapour))+(DepthRate*FN2)*(MoveTime-(1/(math.log(2)/HalftimesN2[cntr])))-((FN2*(CurrentPressure-WaterVapour))-TissueN2[cntr]-((DepthRate*FN2)/(math.log(2)/HalftimesN2[cntr])))*math.exp(-1*(math.log(2)/HalftimesN2[cntr])*MoveTime)
TissueHe[cntr] = (FHe*(CurrentPressure-WaterVapour))+(DepthRate*FHe)*(MoveTime-(1/(math.log(2)/HalftimesHe[cntr])))-((FHe*(CurrentPressure-WaterVapour))-TissueHe[cntr]-((DepthRate*FHe)/(math.log(2)/HalftimesHe[cntr])))*math.exp(-1*(math.log(2)/HalftimesHe[cntr])*MoveTime)

However this assumes FN2 and FHe remain the same throughout the stage - which would only be true OC.

What I can't find are the equivalent equations for ascent/descent phases on CCR where not only is the pressure changing but so also are the FN2 and FHe - which are changing at a different rate to the pressure change...!

Every description of Buhlmann I've come across so far (multiple online including the Baker decolessons.pdf, Deco For Divers etc) only describe the OC equations. Yes, I could approximate it by dividing it in to small discrete steps but I'd prefer to use the exact equations if possible.

To save me trying to sweat it out with a pencil myself, does anyone happen to have the CCR equivalent for Schreiner?

Thanks for any leads...

 

[elmo]

Disclaimer: this hasn't been tested at all, so use at your own risk. I use Subsurface for deco planning, and am reasonably familiar with the deco code. Subsurface uses the Buhlmann equation with small time intervals for both OC and CC. It is computationally less intense than the Schreiner equation, even with the smaller intervals required to get a similar level of precision. The maths below assumes I haven't forgotten calculus completely.

The term DepthRate*FN2 is rate of change (per time) of pN2, which for constant setpoint is equal to DepthRate by FN2InInert (fraction of N2 in inert gas), assuming DepthRate is change in ambient pressure per time. This should be the only part of the equation that is different from OC.

E.g. for 15/55 diluent,
FN2InInert = (1 - 0.15 - 0.55) / (1 - 0.15) = 0.353
FHeInInert = 0.55 / (1 - 0.15) = 0.647

Assuming constant setpoint...

FN2 = (CurrentPressure - SetPoint) / CurrentPressure*FN2InInert
FHe = (CurrentPressure - SetPoint) / CurrentPressure*FHeInInert
TissueN2[cntr] = (FN2*(CurrentPressure-WaterVapour))+(DepthRate*FN2InInert)*(MoveTime-(1/(math.log(2)/HalftimesN2[cntr])))-((FN2*(CurrentPressure-WaterVapour))-TissueN2[cntr]-((DepthRate*FN2InInert)/(math.log(2)/HalftimesN2[cntr])))*math.exp(-1*(math.log(2)/HalftimesN2[cntr])*MoveTime)
TissueHe[cntr] = (FHe*(CurrentPressure-WaterVapour))+(DepthRate*FHeInInert)*(MoveTime-(1/(math.log(2)/HalftimesHe[cntr])))-((FHe*(CurrentPressure-WaterVapour))-TissueHe[cntr]-((DepthRate*FHeInInert)/(math.log(2)/HalftimesHe[cntr])))*math.exp(-1*(math.log(2)/HalftimesHe[cntr])*MoveTime)

Note, you'll have to set a limit on the setpoint in the shallows so that it can't be greater than current pressure.

 

[atdotde]

I derived the equation you are looking for in my blog: Please follow me to CCR Schreiner Equation – The Theoretical Diver

 

[huwporter]

Huge thanks to both for your excellent replies - I'll implement both your solutions tomorrow and let you know the results. I was definitely feeling it has been too many decades since my calculus classes.

In the meantime I realised there is another approach that I was about to post here: absolute pressure is only relevant in the Buhlmann/Schreiner equations for calculating the inspired pressure of your inert gasses. if PO2 is constant, then the total inspired pressure of inert gasses still changes linearly with depth, it just doesn't cross the pressure axis at zero bar, instead at the setpoint; so you can recalculate an 'effective pressure' as far as the inerts are concerned with zero at the absolute pressure of the setpoint and feed that in to the equations instead, and then multiply by the fraction of the ratio of inert gases only. I have implemented that and I'll compare with the above tomorrow.

(And yes, you need to apply a limit in the shallows to avoid calculating a negative pressure on the inerts, cool though it would be in real life to be able to actively suck inert gasses out of the tissues )

 

[huwporter]

As promised, results of comparing the various solutions to the Schreiner equation for CCR:

These are the compartment pressures for example compartments 1,4,7 and 16 at the end of a descent to 90m from surface at 18m/min with 10/70 diluent at setpoint 1.0 (to dodge having to code around the situation where the setpoint is greater than ambient):

OK, I'm scratching my head a little here. The Elmo and Helling solutions are identical, indicating they are mathematically equivalent. My solution is similar enough to the stepwise iterative limit of (my) single depth function, probably within computational error. But these two groups are noticeably different from each other.

My first guess that the difference is down to the way water vapour is treated?

Edit: Yes, the difference is down to water vapour. If I set water vapour to zero in my versions then they give the same values as the elmo/Helling solutions. (But elmo explicitly includes water vapour?)

Reference for checking if I have made any mistakes:

elmo version:

FN2InInert = DilN2/(DilN2+DilHe)
FHeInInert = DilHe/(DilN2+DilHe)

FN2 = (CurrentPressure - SetPoint) / CurrentPressure*FN2InInert
FHe = (CurrentPressure - SetPoint) / CurrentPressure*FHeInInert

for cntr in Compartments:
   TissueN2[cntr] = (FN2*(CurrentPressure-WaterVapour))+(DepthRate*FN2InInert)*(MoveTime-(1/(math.log(2)/HalftimesN2[cntr])))-((FN2*(CurrentPressure-WaterVapour))-TissueN2[cntr]-((DepthRate*FN2InInert)/(math.log(2)/HalftimesN2[cntr])))*math.exp(-1*(math.log(2)/HalftimesN2[cntr])*MoveTime)
   TissueHe[cntr] = (FHe*(CurrentPressure-WaterVapour))+(DepthRate*FHeInInert)*(MoveTime-(1/(math.log(2)/HalftimesHe[cntr])))-((FHe*(CurrentPressure-WaterVapour))-TissueHe[cntr]-((DepthRate*FHeInInert)/(math.log(2)/HalftimesHe[cntr])))*math.exp(-1*(math.log(2)/HalftimesHe[cntr])*MoveTime)

Helling version:

for cntr in Compartments:
   TissueN2[cntr] = ((DilN2/(DilN2+DilHe))*((math.exp(-1*(math.log(2)/HalftimesN2[cntr])*MoveTime))-1)*(SetPoint+DepthRate/(math.log(2)/HalftimesN2[cntr])))+((math.exp(-1*(math.log(2)/HalftimesN2[cntr])*MoveTime))*(TissueN2[cntr]-((DilN2/(DilN2+DilHe)) *CurrentPressure)))+ ((DilN2/(DilN2+DilHe))* (CurrentPressure+(MoveTime*DepthRate)))
   TissueHe[cntr] = ((DilHe/(DilN2+DilHe))*((math.exp(-1*(math.log(2)/HalftimesHe[cntr])*MoveTime))-1)*(SetPoint+DepthRate/(math.log(2)/HalftimesHe[cntr])))+((math.exp(-1*(math.log(2)/HalftimesHe[cntr])*MoveTime))*(TissueHe[cntr]-((DilHe/(DilN2+DilHe)) *CurrentPressure)))+ ((DilHe/(DilN2+DilHe))* (CurrentPressure+(MoveTime*DepthRate)))

My version:

EffectiveCurrentPressure = CurrentPressure - SetPoint
EffectiveFN2 = DilN2/(DilN2+DilHe)
EffectiveFHe = DilHe/(DilN2+DilHe)

for cntr in Compartments:
   TissueN2[cntr] = (EffectiveFN2*(EffectiveCurrentPressure-WaterVapour))+(DepthRate*EffectiveFN2)*(MoveTime-(1/(math.log(2)/HalftimesN2[cntr])))-((EffectiveFN2*(EffectiveCurrentPressure-WaterVapour))-TissueN2[cntr]-((DepthRate*EffectiveFN2)/(math.log(2)/HalftimesN2[cntr])))*math.exp(-1*(math.log(2)/HalftimesN2[cntr])*MoveTime)
   TissueHe[cntr] = (EffectiveFHe*(EffectiveCurrentPressure-WaterVapour))+(DepthRate*EffectiveFHe)*(MoveTime-(1/(math.log(2)/HalftimesHe[cntr])))-((EffectiveFHe*(EffectiveCurrentPressure-WaterVapour))-TissueHe[cntr]-((DepthRate*EffectiveFHe)/(math.log(2)/HalftimesHe[cntr])))*math.exp(-1*(math.log(2)/HalftimesHe[cntr])*MoveTime)

Iterative version (e.g. 1m):

Increment = 0.1
CurrentPressure = 1.0 + (Increment/2)

while CurrentPressure < (MaxDepth/10)+1:

   TimeInterval = ((MaxDepth/10)/DepthRate)/(MaxDepth/(Increment*10))
   PO2 = SetPoint
 
   PN2 = (CurrentPressure-(PO2+WaterVapour))*(DilN2/(DilN2+DilHe))
   PHe = (CurrentPressure-(PO2+WaterVapour))*(DilHe/(DilN2+DilHe))
 
   FN2 = PN2/CurrentPressure
   FHe = PHe/CurrentPressure

   for cntr in Compartments:
       TissueN2[cntr] = TissueN2[cntr]+(PN2-TissueN2[cntr])*(1-math.pow(2,(-1*(TimeInterval/HalftimesN2[cntr]))))
       TissueHe[cntr] = TissueHe[cntr]+(PHe-TissueHe[cntr])*(1-math.pow(2,(-1*(TimeInterval/HalftimesHe[cntr]))))

   CurrentPressure = CurrentPressure + Increment

Thoughts?

(It's definitely worth saying that in the "measure with a micrometer, cut by a madman wildly swinging a blunt axe" of real world decompression, pretty any of these would work just fine)