• Welcome to ScubaBoard


  1. Welcome to ScubaBoard, the world's largest scuba diving community. Registration is not required to read the forums, but we encourage you to join. Joining has its benefits and enables you to participate in the discussions.

    Benefits of registering include

    • Ability to post and comment on topics and discussions.
    • A Free photo gallery to share your dive photos with the world.
    • You can make this box go away

    Joining is quick and easy. Login or Register now by clicking on the button

Schreiner equations for CCR?

Discussion in 'Dive Software' started by huwporter, Nov 29, 2017.

  1. huwporter

    huwporter DIR Practitioner

    # of Dives: 1,000 - 2,499
    Location: Sydney, Australia.
    174
    154
    43
    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):

    Code:
    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:

    Code:
    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...
     
  2. elmo

    elmo DIR Practitioner

    # of Dives: 200 - 499
    Location: Melbourne
    112
    55
    28
    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.
     
    huwporter likes this.
  3. atdotde

    atdotde Solo Diver

    # of Dives: 200 - 499
    Location: Munich, Germany
    93
    70
    18
    huwporter and elmo like this.
  4. huwporter

    huwporter DIR Practitioner

    # of Dives: 1,000 - 2,499
    Location: Sydney, Australia.
    174
    154
    43
    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 :D)
     
  5. huwporter

    huwporter DIR Practitioner

    # of Dives: 1,000 - 2,499
    Location: Sydney, Australia.
    174
    154
    43
    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):

    CCRSchreinerTesting.png


    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:
    Code:
    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:
    Code:
    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:
    Code:
    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):
    Code:
    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)
     
    Last edited: Dec 1, 2017

Share This Page