ACES PSC Module:  Equations File   (Last modified:  30/08/12)

 This file contains a list of all equations used by the analysis routine.
 When performing an analysis calculations are performed in the strict
 order specified within this file.   New parameters can be defined and
 incorporated into the this Equations file but they must first be defined
 in the Parameters List file.   

 EDITING INSTRUCTIONS

 The rules for editing this file and creating and defining equations are
 given below.

 GENERAL RULES

 Since this module mimics a "spreadsheet" equations must be inserted
 into the list in the strict order in which you require them to be
 executed when performing a design.  Consequently they must conform
 to the following rules:

  (1)  Variables must first be defined in the Parameters List file.
  (2)  Blank lines may be used anywhere to delineate equation blocks.
  (3)  Comment lines must be prefaced with the exclamation symbol (!)
  (4)  Each equation may be up to 120 characters long.
  (5)  Allowable operands are ‘+’, ‘-‘, ‘*’, ‘/’ and ‘^’ (exponential)
  (6)  Variables must appear in equations without curly braces.
  (7)  Undefined variables will be allocated the value ‘1’.
  (8)  Any number of bracket groups may be used in an equation.
  (9)  Nested brackets are NOT permitted.
 (10) Spaces in the equation will be ignored.
 (11) Constants may be specified as decimal, integer or exponential
 (12) String (text) variables must contain the '$' symbol (e.g. Note$)
 (13) String variables can only be assigned a maximum of 44 characters
 (14) Comments on an IF statement line must not include the ( ) symbols
 (15) Only equations that lie within the **** boundaries will be used

 ________________________________________________________________________

 DISPLAYING EQUATIONS in the CALCULATION LOG

    .  Equations have been listed sequentially in the equation list to
       follow the logical order in which the design is usually performed.
       Blocks of equations therefore correspond to each tab in the 
       module.  The tabs have been internally numbered 1-9 and the 
       first block of equations in any tab is also given that same 
       number  e.g.   !4   PRESTRESS LOSSES TAB

    .  When the <Check calcs> button at the bottom of a tab is clicked 
       ACES will display the calculation log of all equations associated
       with that tab.

    .  Sub-groups of equations represent various frames and calculation
       buttons within a given tab.  They have been given a label and, 
       where applicable, a letter designating that sub-group  e.g.,
                !7b  SHEAR DESIGN - Flexure shear cracking

    .  Only the equations relating to this block will be displayed in the
       calculation log when the 'View calcs' button on that dialog box 
       (or in that frame) is clicked.

    .  When the <Full calcs> button at the bottom of a tab is clicked
       ACES will display the full calculation log i.e., all equations
       associated with the entire analysis.

    .  A block or sub-block is ended when another tab group number or
       sub-group number is encountered.

    .  !0 signifies that the equations in that block are not to be 
       associated with any button.

    .  !+ signifies that all equations in that block will always be 
       displayed in the log.

    .  !- signifies that equations in that block will not be displayed 
       in the log for subgroups.

 ________________________________________________________________________

 GENERAL FUNCTIONS

    .  The MAX and MIN functions 

       The MAXimising and MINimising functions operate on variables 
       and numbers only - no expressions are allowed.  Variables and
       numbers may be signed  (e.g.  -ve) and up to 20 values only are
       allowed within the brackets.  Note that the ABS and TAN 
       functions CANNOT be included in a MAX/MIN list.

    .  TAN (tangent) function.  This takes the form:
                              var = TAN (Angle)
       'Angle' must be specified in radians.  If it is not already in radians
       then convert it using one of the nine temporary work variables
       (WorkVarX) described in the Parameters List file and the two
       expressions:
                              WorkVarX = (Angle*3.14159/180.)
                              var = TAN (WorkVarX)
       Note that the ABS and TAN functions CANNOT be included together
       on the same line nor in a MAX/MIN list.  In such cases a
       temporary parameter (WorkVarX) must be used to calculate an
       intermediate value.

    .  ABS (absolute value) function.  This takes the form:
                              var = ABS (xxx)
       where 'var' is any parameter in the parameter list and 'xxx' is either
       an optionally signed variable or a number.  It MUST be included
       within normal brackets (as shown).   Note that the ABS and TAN
       functions CANNOT be included together on the same line nor in a
       MAX/MIN list.  In such cases a temporary parameter (WorkVarX)
       must be used to calculate an intermediate value.  Note also that in
       a statement such as var = -ABS (xxx)  var will always be negative.

    .  OUTPUT function

       The function is specified as:  OUTPUT a, b, c, d, .... etc  where a,
       b, c, etc  are parameters that have already been defined and used
       in the equation list.  ACES will simply display the specified values,
       separated by commas.   This can be very useful in LOOP and
       ITERATE functions to track the behaviour of one or more variables.
       A maximum of 20 parameters can be listed in OUTPUT ......

       Note that the display of values will occur even if the NOTRACK
       flag has been invoked (refer to the LOOP function below for a
       description of the NOTRACK flag).


    .  IF function.  There are two types of IF function:

       (1) Single line form:  The function can be expressed either as:
                       IF (xxx > yyy)  expression
                                       or as:
                       IF (xxx > yyy) THEN expression

       (2) Multiple line form.  This must be expressed as:
                       IF (xxx > yyy) THEN
                            expressions1
                       ELSE
                            expressions2
                       END IF

            Both expressions1 and expressions2 may contain up to 20
            statements (or none at all) and the operands can only be:
            ">", "<" or "=".   Values 'xxx' and 'yyy' can only be variables
            or numbers and must be included within normal brackets.
            Comments on an IF statement line must not include embedded
            quotation ( ) symbols.

       NOTE that only the Single line form can be used within a loop


    .  LOOP function.  This takes the form:

            LOOP  i = a, b NOTRACK
                  expression 1
                  function 2
	   ...... etc ......
            END LOOP

       NOTRACK is an optional attribute - it suppresses the tracking and
       display of intermediate results during the looping process.

       If any parameter name within 'expression' contains the character 'i'
       then the 'i' itself will be replaced with the number corresponding 
       to the loop variable 'i'.

       For example, in the calculation of tension forces in PS strands the
       general equation for setting the bar depth (Dpbari) of bars in row
       'i' is:  “Dpbari = D – Ybarri”.  Using the LOOP function, this can
       be performed for bar rows 1 to 8 using three lines of code  viz:

             LOOP  i = 1,8
                   Dpbarri = D - Ybarri
             END LOOP

       Internally ACES will create 8 equations of the form:
       Dpbarr1 = D - Ybarr1, Dpbarr2 = D - Ybarr2, etc.

       A maximum of 20 expressions and/or functions can be included in
       the loop and the loop variables 'a' and 'b' can only lie betweeb 0 
       and 20

      NOTE that only the Single line form of IF can be used within a loop

    .  ITERATE funtion.  This must take the form:

        ITERATE
          .........
        END ITERATE

              or

        ITERATE ON param IN STEPS OF xx  NOTRACK
            ......  (up to twenty statements allowed) .....
        END ITERATION WHEN (param ><= number) OR AFTER n CYCLES

        where:  'param' = any parameter in the parameter list

                         'xx' = the number to be added to 'param' at the
                                 beginning of each cycle. This may be zero 
                                 if the iteration is naturally convergent.

                         ' n' = maximum number of cycles to perform if the
                                 iteration is not converging quickly enough.

       The quantity evaluated, (param ><= number) must be enclosed in
       normal brackets.  Note that the tolerance for equality has been
       set to 1%.  That is, if two numbers are within 1% of each other,
       they are considered equal. 'number' may be a parameter or a 
       number.

       NOTRACK is optional and suppresses the tracking and display of
       intermediate results during the iteration process.

       Example:

         ITERATE ON WorkVar1 IN STEPS OF 5
                WorkVar2 = WorkVar1 * 2.
         END ITERATION WHEN (WorkVar1 = 25) OR AFTER 30 CYCLES

	For notes (TnoteXX) involving the Torsion equations, prepending 
        the string with a ? will display the note in red, prepending 
        with ^ will display in green.

 ________________________________________________________________________

SPECIAL FUNCTIONS

    .  SORT functions:
                       SORT_LOW_HIGH (p1i,p2i,r1i,r2i) NOTRACK
                       SORT_HIGH_LOW (p1i,p2i,r1i,r2i) NOTRACK

       These two functions sort a set of 8 values 'p1i' in ascending and/or
       descending order respectively and place the sorted results into the
       set labelled 'r1i'.  Values in set 'p2i' correspond to those in set 
       'p1i' and are placed into 'r2i' during the sort process.  In effect
       the function sorts the set 'p1i' and its corresponding values 'p2i' 
       and stores the results into 'r1i' and 'r2i'.  Note, however, that 
       'p1i' and 'p2i' must contain 8 values each.

       NOTRACK suppresses the tracking and display of sorted output.
       The functions have been set up specifically for sorting strand and
       passive reinforcement bars and their corresponding forces in order
       to determine if they lie in the tension zone of the section.

       Example:      SORT_LOW_HIGH (Ybarri, Fbari, Dsorti, Fsorti)


    .  COMBINE and SORT functions:

              COMBINE_SORT_LOW_HIGH (p1i, p2i, p3i, p4i, r1i, r2i, r3i)
              COMBINE_SORT_HIGH_LOW (p1i, p2i, p3i, p4i, r1i, r2i, r3i)

       This is similar to the SORT function described above with the
       exception that the 'p1i' and 'p3i' data sets and their corresponding
      'p2i' and 'p4i' sets are first combined before sorting.  The 8 values 
       in each of the individual 'p1i' and 'p3i' sets then appear as 16 
       sorted values in the final 'r1i' set.  

       This function is used to combine the strand table ('p1i'/'p3i') and 
       the passive reinforcement tables into one consolidated table which 
       is then sorted by bar distance ('r1i') and the corresponding bar 
       forces ('r2i').

       The values in 'r3i' are flags that indicate which set the 'r1i' values
       originally came from (either 'p1i' or 'p3i').  A value of 1 indicates 
       'p1i' and 2 indicates 'p3i'.  This information is required by the 
       strain compatibility routine to differentiate a PS strand from a 
       reinforcing bar when calculating the yield strain.

       Example:

       COMBINE_SORT_HIGH_LOW (Dbari,Fbari,Dpfi,Fbfi,Dsti,Fsti,Typei)


    .  LOCATE Ultimate Neutral Axis Depth (dn) Functions:

                     LOCATE_DN_USING_SIMPLIFIED
                                           and
                     LOCATE_DN_USING_STRAIN_COMPATIBILITY

       Both functions use the values returned by the COMBINE_SORT
       functions.  They locate the neutral axis depth (Dcb) and populate
       the bar force table automatically.  Although these functions could
       have been replaced by inserting relevant equations directly into the
       equations file it has been found that the hard-coded routines
       implemented here-in execute much faster.

       The '_SIMPLIFIED' method assumes:  (1) a concrete stress block of
       0.85f'c; and (2) that all bars and strands have yielded (except 
       those in the compression zone, which are excluded from the analysis).

       The '_STRAIN _COMPATIBILITY' method derives the concrete 
       ultimate strain from the parameter 'ucu' which is editable on the
       'Ult M Check' tab but which defaults to 0.003.  The function 
       incrementally adjusts the depth Dcb until the concrete force 
       balances the steel forces.  Bars and strands in the compression 
       zone are excluded and the bar/strand forces are reduced if the
       bar/strand strain is below yield.  (This is commonly refered to as
       'Partial Prestressing').


***************************** Start of Equations *****************************

!-  CALCULATION OF DEFAULT VALUES

IF (f`cg = 0.) THEN f`cg = 50.
IF (Eg = 0.) THEN Eg = 0.043 * (2500.^1.5) * (f`cg^0.5)
IF (f`cs = 0.) THEN f`cs = 40.
IF (Es = 0.) THEN Es = 0.043 * (2400.^1.5) * (f`cs^0.5)
IF (f`cmi = 0.) THEN f`cmi = 40.
IF (Egmi = 0.) THEN Egmi = 0.043 * (2500.^1.5) * (f`cmi^0.5)
IF (f`cmt = 0.) THEN f`cmt = 30.
IF (u5 = 0.) THEN u5 = us
IF (u6 = 0.) THEN u6 = u1
IF (fcrack = 0.) THEN fcrack = 200.

IF (Egmt = 0.) THEN 
    WorkVar1 = 0.8 * f`cg
    Egmt = 0.043 * (2500.^1.5) * (WorkVar1^0.5)
END IF

IF (Esr = 0.) THEN Esr = Ep
IF (ø = 0.) THEN ø = 0.9
IF (O = 0.) THEN O = 0.9
IF (øm= 0.) THEN øm= 0.9
IF (Om= 0.) THEN Om= 0.9
IF (ucu = 0.) THEN ucu = 0.003


!2   SECTION PROPERTIES TAB - Section property calculations

Zt = Ig / (Dg - Yb)
Zb = Ig / Yb

IF (Ws = 0.) THEN Ws = Wtf
IF (bef = 0.) THEN bef = Ws
IF (Es = 0.) THEN Es = n*Eg

n  =  Es/Eg
Bw2 = 2*bw
bv = Bw2

IF (CalcProp = 1) THEN
   CalcPrp$ = Composite section props recalculated by ACES
   D = Dg + Ts
   As = Ws * Ts
   Ys = Dg + Ts/2.
   Ac = Ag + As*n
   Yc = (Ag * Yb + As*n*Ys) / (Ag + As*n)
   WorkVar1 = Yc - Yb
   WorkVar2 = Ys - Yc
   Ic = Ig + Ag*WorkVar1^2 + n*Ws*Ts^3/12. + n*As*WorkVar2^2
   WorkVar1 = D - Yc
   Zst = Ic / (WorkVar1*n)
   WorkVar1 = Dg - Yc
   Zsb = Ic / (WorkVar1*n)
   Zgt = Ic / (WorkVar1)
   Zgb = Ic / Yc
ELSE
   CalcPrp$ = Composite props NOT recalculated by ACES 
END IF

!  Calculation of shear flow at centroid and web/flange

IF (CalcProp = 1) THEN
   Ybflan = (Yc - Ybf)                                   ! Bottom flange to composite centroid
   WorkVar1 = (Yc - Dwf)                                 ! Height of web block
   Areac = Abf + WorkVar1 * Bw2                          ! Shear area below composite centroid
   Ymomsc = Abf * Ybflan + WorkVar1 * Bw2* WorkVar1 / 2. ! MOI below composite centroid
   Ybarc = Ymomsc / Areac                                ! Centroidal offset of section below centroid
   Qna = Areac * Ybarc                                   ! Shear flow at centroid     
!    Calculations for Web/flange joint
   Hwebfl = Dwf - Tbf                                    ! Height of web block
   Ywebfl = Yc - Tbf - Hwebfl / 2.                       ! Centroidal offset of web block
   Areac = Abf + Hwebfl*Bw2                              ! Area of bottom flange + web block
   Ymomsc = Abf*Ybflan + Hwebfl*Bw2*Ywebfl               ! Moment of area of section below centroid
   Ybarc = Ymomsc / Areac                                ! Centroidal offset of section below centroid
   Qfw = Areac * (Yc-Dwf)                                ! Shear flow at web/flange interface    
END IF

IF (Qna < 0.) THEN Qna = -Qna
IF (Qfw < 0.) THEN Qfw = -Qfw


!2a   SECTION PROPERTIES TAB - Non-prestressed reinforcement

Nrfbarst = 0
  Arft     = 0
Ynbarsrf = 0

LOOP  i = 1,20
  Nrfbarst = Nrfbarst + Nrfbrsi
END LOOP

LOOP  i = 1,20
  WorkVar1 = 3.14159 * BarDtri * BarDtri / 4
  Arfi = WorkVar1 * Nrfbrsi
END LOOP

LOOP  i = 1,20
  Arft = Arft + Arfi
END LOOP

Arl = Arft                                   ! Total area of passive reinforcement

LOOP  i = 1,20
  YbxNrfi = Yrfi * Arfi
END LOOP

LOOP i = 1,20
  Ynbarsrf = Ynbarsrf + YbxNrfi
END LOOP


!3a PS ESTIMATES TAB - Differential Shrinkage

u = u2 - u1
n = Es/Eg
ec = Ys-Yc

Fshr = u*Es*As*1E-9*(1 - 2.71828^-Rcf)/Rcf   ! Calculation of shrinkage force
Mshr = Fshr*ec/1000
SumNA = n*As + Ag
Fsus = Fshr*1000*(1/As - n/SumNA)
Fsug = -Fshr*1000/SumNA

fts = Fsus - Mshr*1E6/Zst
fbs = Fsus - Mshr*1E6/Zsb
ftg = Fsug - Mshr*1E6/Zgt
fbg = Fsug + Mshr*1E6/Zgb


!3b PS ESTIMATES TAB - Prestress Forces
! Previous option to included debonded bars has been removed

Nbarb1 = Nbart1 - Nbard1
Nbarb2 = Nbart2 - Nbard2
Nbarb3 = Nbart3 - Nbard3
Nbarb4 = Nbart4 - Nbard4
Nbarb5 = Nbart5 - Nbard5
Nbarb6 = Nbart6 - Nbard6
Nbarb7 = Nbart7 - Nbard7
Nbarb8 = Nbart8 - Nbard8

IncDBar$ = Debonded bars are NOT included in analysis


YbxNb1 = Ybarr1 * Nbarb1
YbxNb2 = Ybarr2 * Nbarb2
YbxNb3 = Ybarr3 * Nbarb3
YbxNb4 = Ybarr4 * Nbarb4
YbxNb5 = Ybarr5 * Nbarb5
YbxNb6 = Ybarr6 * Nbarb6
YbxNb7 = Ybarr7 * Nbarb7
YbxNb8 = Ybarr8 * Nbarb8

!   Find total number of bars in layers 1 to 8 (Nbart1 to Nbart8)

Nbars = 0

LOOP  i = 1,8
  Nbars = Nbars + Nbarti
END LOOP

!   Find total number of debonded bars in layers 1 to 8

Ndbars = 0

LOOP  i = 1,8
  Ndbars= Ndbars+ Nbardi
END LOOP

IF (Nbars = 0) THEN
   Ntubbars = 0.
ELSE
   Ntubbars = 100. * Ndbars / Nbars
END IF

!   Find total number of bars used in the analysis (may/may not include debonded bars)

Nbbars = 0

LOOP  i = 1,8
  Nbbars= Nbbars+ Nbarbi
END LOOP

!   Find total moment of area of bars in each layer

Ynbars = 0

LOOP  i = 1,8
  Ynbars= Ynbars+ YbxNbi
END LOOP

!   Calculate Pj & Centroid of bar group & Eccentricity of group

Pj = Nbbars*Pult*Jf
Ycgs = Ynbars / Nbbars
e = Yb - Ycgs


!4  PRESTRESS LOSSES TAB - Steam Relaxation Losses

k5a = 1 + (Jf-0.7)*0.5/0.1
k5b = (Jf-0.4)/0.3
WorkVar3 = 0.
k5 = MAX (k5a,k5b,WorkVar3)

Lsrl = 0.1*k5/1.5
Lsr = 100*Lsrl
Prl = -Lsrl*Pj
Pjr = Pj + Prl


!   PRESTRESS LOSSES TAB - Elastic Deformation Losses

Ap = Nbbars * Aps
Apt= Ap
fcgs = -Pjr*1000*(1/Ag + e*e/Ig) + Msw*1E6*e/Ig
Pelastic = fcgs*Ep*Ap/(Egmt*1000)
Ledl = -Pelastic*100/Pj
Pt = Pjr + Pelastic
Ltr = Pt*100/Pj

!   PRESTRESS LOSSES TAB - Shrinkage Losses

Ac= n*As + Ag
Pshr = -us*Ep*Ap*1E-9/(1 + 15*Arl/Ac)
Lshr = -Pshr*100/Pj
Prs = Pt + Pshr

!   PRESTRESS LOSSES TAB - Creep loss parameters

At = Ag + Ws * Ts
th = 2*At/(Gp + 0.5*Vp)
Fratio = f`cmt/f`cg


!4e PRESTRESS LOSSES TAB - Creep losses due to PS + SW

øcc = øccb*k2*k3
Occ = øcc                                            ! HTML Reports cannot handle ø symbol
Occb = øccb                                          ! HTML Reports cannot handle ø symbol
O = ø                                                ! HTML Reports cannot handle ø symbol
fcscgs = -Pt*1000/Ag - Pt*1000*e*e/Ig + Msw*1E6*e/Ig
ucc1 = 1.E6*fcscgs*øcc/(Eg)


!4f PRESTRESS LOSSES TAB - Creep losses due to Deck + SDL

Fdeck = Mslab*1E6*(Yb-Ycgs)/Ig
Fsdl  = Msdl * 1E6*(Yc-Ycgs)/Ic
øcc2  = øccb * k2s*k3s
Occ2  = øcc2                                         ! HTML Reports cannot handle ø symbol
ucc2  = øcc2*1E6*(Fdeck+Fsdl)/Eg

!   PRESTRESS LOSSES TAB - Total creep losses

ucc    =  ucc1 + ucc2
Pcreep = ucc*Ep*Ap/1E9
Lcr    =  -Pcreep*100/Pj
P      =  Pt+Pshr+Pcreep
Ltt    =  P*100/Pj


!5  SERVICEABILITY CHECK - Girder stresses at transfer  (Tension +ve  Compression -ve)

f`csat = 0.5*f`cmt^0.5
f`csac = 0.6*f`cmt
f`at   = 0.25*f`cg^0.5
f`ac   = 0.4*f`cg
f`atsv   = 0.5*f`cg^0.5
f`acsv   = 0.6*f`cg

Mpte   = -Pt*e/1000

fgt1   = -Pt*1000/Ag    
fgb1   =  fgt1
fgt2   = -Mpte*1E6/Zt
fgb2   =  Mpte*1E6/Zb
fgt3   = -Msw *1E6/Zt
fgb3   =  Msw*1E6/Zb
fgt4   =  fgt1+fgt2+fgt3
fgb4   =  fgb1+fgb2+fgb3

!5b SERVICEABILITY CHECK - Design moment stresses (Tension +ve;  Compressn -ve)

Mpe    = -P*e/1000
fgtss1 = -P*1000/Ag
fgbss1 =  fgtss1
fgtss2 = -Mpe*1E6/Zt
fgbss2 =  Mpe*1E6/Zb
fgtss3 = -Msw*1E6/Zt
fgbss3 =  Msw*1E6/Zb
fgtss4 = -Mslab*1E6/Zt
fgbss4 =  Mslab*1E6/Zb
fstss5 = -Msdl*1E6*SDLf/Zst
fsbss5 = -Msdl*1E6*SDLf/Zsb
fgtss5 = -Msdl*1E6*SDLf/Zgt
fgbss5 =  Msdl*1E6*SDLf/Zgb
fstss6 =  fts
fsbss6 =  fbs
fgtss6 =  ftg
fgbss6 =  fbg
fstss7 = -Mll*1E6/Zst
fsbss7 = -Mll*1E6/Zsb
fgtss7 = -Mll*1E6/Zgt
fgbss7 =  Mll*1E6/Zgb
fstss8 = -Mhvl*1E6/Zst
fsbss8 = -Mhvl*1E6/Zsb
fgtss8 = -Mhvl*1E6/Zgt
fgbss8 =  Mhvl*1E6/Zgb

!   SERVICEABILITY CHECK - Final design DL+SDL+LL stresses

fstll  =  fstss5+fstss6+fstss7+ftempst
fsbll  =  fsbss5+fsbss6+fsbss7+ftempsb
fgtll  =  fgtss1+fgtss2+fgtss3+fgtss4+fgtss5+fgtss6+fgtss7+ftempgt
fgbll  =  fgbss1+fgbss2+fgbss3+fgbss4+fgbss5+fgbss6+fgbss7+ftempgb

!   SERVICEABILITY CHECK - Final design DL+SDL+HLV stresses

fstsv  = fstss5+fstss6+fstss8+ftempst
fsbsv  = fsbss5+fsbss6+fsbss8+ftempsb
fgtsv  = fgtss1+fgtss2+fgtss3+fgtss4+fgtss5+fgtss6+fgtss8+ftempgt
fgbsv  = fgbss1+fgbss2+fgbss3+fgbss4+fgbss5+fgbss6+fgbss8+ftempgb

!   SERVICEABILITY CHECK - Determine maximum values

Msdlf = Msdl*SDLf 

fmaxsllt = MAX (fstll,fsbll,fgtll,fgbll)
IF (fmaxsllt < 0.) THEN fmaxsllt = 0.

fmaxsllc = MIN (fstll,fsbll,fgtll,fgbll)
IF (fmaxsllc > 0.) THEN fmaxsllc = 0.

fmaxssvt = MAX (fstsv,fsbsv,fgtsv,fgbsv)
IF (fmaxssvt < 0.) THEN fmaxssvt = 0.

fmaxssvc = MIN (fstsv,fsbsv,fgtsv,fgbsv)
IF (fmaxssvc > 0.) THEN fmaxssvc = 0.


!5c SERVICEABILITY CHECK - Cracking Check

!   Perform the SLS check using the linear-elastic strain compatibility method
!  
!   Four conditions are analysed:
!
!   - cracking of PC girder only
!   - PC girder under self weight + slab
!   - cracking of composite girder (incl selfwt+slab strains in PC girder)
!   - composite girder serviceability design moment (incl strains from non-comp PC girder)
!
!  This module creates the following parameters
!
!  curvco  , Composite girder curvature
!  curvpc  , Precast curvature
!  fcocrps , Composite cracking strand stress
!  fcocrreo, Composite cracking reinf. stress
!  fcodlps , Composite moment strand stress
!  fcodlreo, Composite moment reinf. stress
!  fcoincps, Composite strand stress increment
!  fcoincre, Composite reinf. stress increment
!  fpccrps , Precast cracking strand stress
!  fpccrreo, Precast cracking reinf. stress
!  fpcdlps , Precast strand stress
!  fpcdlreo, Precast reinf. stress
!  fpcincps, Precast strand stress increment
!  fpcincrf, Precast reinf. stress increment
!  Mcocrack, Composite cracking (decompression) moment
!  Mcodesig, Composite moment
!  Mcoincr , Composite moment increment,
!  Mpccrack, Precast cracking (decompression) moment
!  Mpcdlsla, Precast moment - SW + deck
!  Mpcincr , Precast moment increment
!  ucocrbot, Composite cracking strain at bottom
!  ucocrtop, Composite cracking strain at top
!  ucodlbot, Composite strain at bottom
!  ucodltop, Composite strain at top
!  uconeutr, Composite neutral axis depth
!  upccrbot, Precast cracking strain at bottom
!  upccrtop, Precast cracking strain at top
!  upcdlbot, Precast strain at bottom
!  upcdltop, Precast strain at top
!  upcneutr, Precast neutral axis depth
!
!  Call the routine that performs the SLS analysis using strain compatibility

   ANALYSE_SLS_BY_LINEAR_ELASTIC_STRAIN

 fslscrck = MAX (fcoincps,fcoincre)                 ! Determine max stress increment

IF (fslscrck > fcrack) THEN
   CrkNote$ = Crack limit exceeded (fslscrck > fcrack)
ELSE
   CrkNote$ = Crack check OK (fslscrck < fcrack)
END IF
                    
IF (Mcodesig > Mcocrack) THEN
   MslNote$ = SLS design moment exceeds capacity
ELSE
   MslNote$ = SLS moment capacity OK
END IF

IF (SLScase = 0) THEN
   SLScase$ = SLS Crack Check for Live Load
ELSE
   SLScase$ = SLS Crack Check for Special Vehicle Load
END IF


!6b ULTIMATE MOMENT CHECK - Input data & constants 

Dp = D - Ycgs
fp = 1000.*Pult/Aps

!   Set parameters to generate stress-strain curve for strain compatibility

IF (fpus = 0.) fpus = 1750.
IF (epus = 0.) epus = 0.035
IF (epyf = 0.) epyf = 0.85
IF (epy  = 0.) epy = epyf*fpus/Ep

fpy = Ep*epy
epe = P/(Nbbars*Aps*Ep)

WorkVar1 = Ap * fp
WorkVar2 = (Ast-Asc) * fsy
k2u = (WorkVar1 + WorkVar2) / (Ws*Dp*f`cs)

v = Srf - 0.007 * (f`cs-28)
fpu = fp * (1-k1u*k2u/v)
fuc = Srf * f`cs

!   Tension forces in steel

LOOP i = 1,8                       ! Strand depths
   Dpbari = D - Ybarri
END LOOP

LOOP i = 1,20                      ! Passive reinforcement depths
   Dprfi = D - Yrfi
END LOOP

LOOP i = 1,8                       ! Strand areas
   Apbari = Nbarbi * Aps
END LOOP

!   Passive reinforcement areas (Arfi) are already calculated

LOOP i = 1,8    NOTRACK            ! Forces in strand (for simplified method)
   Fbari = Apbari * fpu/1000.
   IF (Fbari = 0.) THEN Dpbari = 0.
END LOOP

LOOP i = 1,20  NOTRACK             ! Forces in passive R/F (for simplified method)
   Fbrfi = Arfi * fsy / 1000.
   IF (Fbrfi = 0.) THEN Dprfi = 0.
END LOOP

!   LOCATE POSITION of NA (dn) & STEEL/CONCRETE CAPACITIES
!
!   Combine and sort prestressed and non-prestressed R/F
!   in descending order of distance from top fibre

    COMBINE_SORT_HIGH_LOW (Dpbari,Fbari,Dprfi,Fbrfi,Dsorti,Fsorti,Typei) NOTRACK

!   Calculate depth to neutral axis (dn) & steel/concrete capacities.
!   Use values stored in Dsorti and Fsorti and base capacity
!   calculations on the type of method selected (see below).
!
!  SIMPLIFIED METHOD - Assumptions
!
!  1. The ultimate concrete stress state is assumed to be
!     rectangular down to a distance of:  gamma*neutral axis depth
!     (viz  gammaDn = v*dn  where v = gamma  and dn = depth to NA)
!
!  2. The ultimate concrete compressive stress Fuc = 0.85F’c
!
!  3. All steel (including passive reinforcement) below the neutral
!     axis is at yield.  Steel in compression is ignored.
!
!  4. The ultimate condition occurs when the compressive force in
!     the concrete balances the tensile forces in the steel.  
!
!
!  STRAIN-COMPATIBILITY METHOD - Assumptions
!
!  1. The ultimate concrete stress state is assumed to be
!     rectangular down to the depth of the equivalent concrete
!     compression block (gamma * Dcb).  Dcb is increased in increments
!     of 0.1mm from an initial value of 0.1mm.
!
!  2. The ultimate concrete strain at the top surface is given by
!     'ucu' (usually 0.003 microstrain).  The strain at each steel level
!     is determined by linear interpolation of a straight line running
!     from 'ucu' and through the neutral axis.  The strain in the compression
!     reinforcement is limited to the maximum allowed for the concrete.
!
!  3. If the interpolated strain exceeds the rupture strain the strand is
!     assumed to be ruptured.  If the interpolated strain exceeds the yield
!     strain for its particular grade of steel, the full yield force for
!     that bar is included.
!
!  4. If the interpolated strain is less that the yield strain, the
!     elastic force at that strain is used.
!
!  5. The total concrete force may include contributions from the slab, 
!     flanges and web depending on the location of gamma*dn.  Ultimate
!     strand stresses are based on the stress-strain interaction diagram.
!
!  6. The ultimate condition occurs when the compressive force in
!     the concrete balances the tensile forces in the steel.  If not,
!     increase the value of Dcb and iterate.
!
!  7. For legacy reasons, two parameters are used to represent the
!     neutral axis depth - Dcb and dn.
!   
!  8. Note that the concrete forces (Fultcs, Fultcf, Fultcw) and their
!     correpsonding moments (Multcs, Multcf, Multcw) are calculated
!     within the LOCATE_DN_.... functions, as is the value of 'gammaDn'.

IF (CalcDcb = 0) THEN
   DcbNote$ = Based on simplified method
   CalcDcb$ = dn calculated using simplified method
   LOCATE_DN_USING_SIMPLIFIED_METHOD
ELSE
   DcbNote$ = Based on strain compatibility method
   CalcDcb$ = dn calc'd using strain compatibility method
   LOCATE_DN_USING_STRAIN_COMPATIBILITY
END IF

WorkVar1 = Ts + Ttf
DnNote$ = NA in slab

IF (dn > Ts)       THEN DnNote$ = Na in flange
IF (dn > WorkVar1) THEN DnNote$ = NA in Web
IF (CalcDcb = 0)   THEN DnNote$ = Simplified method assumes dn in slab

!   Calculate area of passive reinforcement in tension zone

Ast = 0.

LOOP i = 1,20
   WorkVar1 = Arfi
   IF (Dprfi > Dcb) Ast = Ast + WorkVar1
END LOOP

Ftstrand = 0.

LOOP i = 1,8
   IF (Fbari = 0.)  Yelds$i =
   IF (Ybarri > Yc) Ftstrand = Ftstrand + Fbari
END LOOP

LOOP i = 1,8
   ubari = 0.
   IF (Apbari > 0.) ubari = 1E9 * Fbari / (Ep*Apbari)
END LOOP

LOOP i = 1,20
   IF (Fbrfi = 0.)  Yeldr$i =
END LOOP

LOOP i = 1,20
   ubrfi = 0.
   IF (Arfi > 0.) ubrfi = 1E9 * Fbrfi / (Esr*Arfi)
END LOOP

!   Calculate area of non-P/S steel in compression

Asc = Arft - Ast
Dcbnok  = Dcb
Dcbnok1 = Dcb
Dcbnok2 = Dcb

!   Sum strand and R/F forces in WorkVar1 and Moment in WorkVar2

Fsteel = 0.
Msteel = 0.

LOOP i = 1,8 NOTRACK
   Fsteel = Fsteel + Fbari
   Msteel = Msteel + Mbari
END LOOP

LOOP i = 1,20 NOTRACK
   Fsteel = Fsteel + Fbrfi
   Msteel = Msteel + Mbrfi
END LOOP

!   Calculate Ku

IF (Fsteel = 0) THEN
   Ku = 0.
   WorkVar3 = 0.
ELSE
   WorkVar3 = 1000.*Msteel/Fsteel
END IF

IF (WorkVar3 = 0) THEN
   Ku = 0
ELSE
   Ku = dn/WorkVar3
END IF

IF (Ku > 0.4) THEN
   KuNote$ = ? Section is over-reinforced (Ku > 0.4)
ELSE
   KuNote$ = ^ Passive reinforcement OK (Ku < 0.4)
END IF


!6c ULTIMATE MOMENT CHECK - Final Ultimate Design Moments

Musw = Msw * LFsw
Muslab = Mslab * LFslab
Musdl = Msdl * LFsdl
Mull = Mll * LFll
Muhvl = Mhvl * LFhvl
Mu1 = Musw + Muslab + Musdl + Mull
Mu2 = Musw + Muslab + Musdl + Muhvl
Mu = Mcs + Mc

øMu =  øm * Mu
Om = øm                          ! HTML report can`t handle ø symbol
OMu  = øMu                       ! HTML report can`t handle ø symbol
Mmax = MAX (Mu1,Mu2)
Mllcap = øMu - Musw - Muslab - Musdl


!7  SHEAR DESIGN - Ultimate Design Parameters

WorkVar1 = LFsw*Vsw + LFslab*Vslab + LFsdl*Vsdl + LFll*Vll

IF (Vult = 0) THEN Vult = WorkVar1

st = 0.33*f`cg^0.5
f`cf = 0.6*f`cg^0.5


!7f   Calculate Cracking Moment (Mcr) and Muomin

Mcr = 1E-6*Zgb*(f`cf + 1E3*P/Ag) + P*e*1E-3 
Muomin = 1.2 * Mcr 

! Suppress calculation of cracked moment of inertia & effective moment of inertia
!
! WorkVar1 = D - Ycgs
! Iccr = 0.3333*Ws*dn^3 + (Ep/Eg)*Ap*WorkVar1^2
! 
! WorkVar1 = MAX(Mll,Mhvl)
! Msv = Msw + Mslab + Msdl + WorkVar1
! IF (Msv = 0.) Msv = Mcr
!
! WorkVar1 = (Mcr/Msv)
! Ief = Iccr + (Ic - Iccr)*WorkVar1^3


!7b  SHEAR DESIGN - Flexure shear cracking

! Decompression moment & corresponding shear

scpf = -1000*P/Ag - 1000*P*e*Yb/Ig
Mo = -scpf*Zgb/10^6                        

IF (Mult = 0.) THEN
   Vult = 0.
ELSE
   Vo = Mo*Vult/Mult
END IF

! Calculate distance from the extreme compressn fibre to the
! centroid of the outermost layer of tensile R/F or tendons

SORT_LOW_HIGH (Ybarri, Fbari, Dsorti, Fsorti)

Ybarmin = Dsort1
do = D - Ybarmin
Def = D
doMin = 0.8*D

Beta1 = 1.1*(1.6 - do/1000)
IF (Beta1 < 1.1) THEN Beta1 = 1.1
Beta2 = 1.0     !  For Super-Ts (AS5100 Clause 8.2.7.1)

WorkVar1 = Span/2.

IF (x > WorkVar1) THEN
   WorkVar2 = Span - x
ELSE
   WorkVar2 = x
END IF

IF (Beta3 = 0) THEN Beta3 = 2*do*1000/WorkVar2
IF (Beta3 < 1) THEN Beta3 = 1.0
IF (Beta3 > 2) THEN Beta3 = 2.0

WorkVar1 = (Ast+Ap)*f`cg/(bv*do)
WorkVar2 = WorkVar1 ^ 0.333

Vuct = Beta1*Beta2*Beta3*bv*do*WorkVar2/1000.
Vucs = Vuct + Vo


!7c  SHEAR DESIGN - Web shear cracking at centroid of composite section (direct solution)

sc = -P*1000/Ag + P*e*1000*(Yc-Yb)/Ig
Qlbwc = Qna/(Ic*Bw2)

IF (sc < st) THEN
   tc = (st*st-st*sc)^0.5              ! Shear stress tau
ELSE
   tc = 0.
END IF

WorkVar1 = 0.5*sc
WorkVar2 = WorkVar1*WorkVar1
s1c = (WorkVar2 + tc*tc)^0.5 + WorkVar1
Vtc = 1E-3*tc/Qlbwc

Mtc = Mult*Vtc/Vult                    ! Moment capacity corresponding to shear capacity Vtc

IF (Mtc > Mcr) THEN                    ! Section is cracked if Mtc > cracking moment Mcr
    VuccNte$ = Cracked
    Vucc = 0
ELSE
    VuccNte$ = OK
    Vucc = Vtc
END IF


!7d  SHEAR DESIGN - Web shear cracking at the bottom of the web-flange interface
!
!  The algorithm to determine the ultimate concrete shear capacity,
!  Vtf, is based on RF Warner, BV Rangan, AS Hall and KA Faulkes,
!  "Concrete Structures" (Addison Wesley Longman, Aust).  Assume:
!
!   y2 = Dist from web/flange interface to composite centroid
!   y3 = Dist from web/flange interface to centroid of precast girder
!  Dwf = Distance from bottom of girder to web/flange interface
!    P = Final prestress force
!    e = Eccentricity (CG girder-strand group)
!   Ag = Area of precast girder
!   Ig = Moment of Inertia of precast girder
!   Ic = Composite moment of inertia
!   sf = Flexural stress at web-flange interface for uncracked section
!   st = Allowable principal tensile strength
!   tf = Ultimate shear stress capacity of concrete on its own (tau)
!   Bw2= Sum of widths of both webs
!  Qfw = Shear flow constant at flange/web junction
!  Mcf = Corresponding moment factor (= Mult/Vult)
!        (used to calculate the corresponding flexural stress 
!         when the shear capacity is reached)
!
!  The solution for Vtf is based on the following relationships:
!
!  sf = (-P/Ag) - (P*e*y3/Ig) + Vtf*Mcf*y2/Ic .......... (1) Flexural stress
!
!  st = 0.33SQRT(f'cg) = SQRT((sf/2)^2 + tf^2) + sf/2 .. (2) Tensile stress capacity
!
!  Vtf = tf*Ic*Bw2/Qfw                      ............. (3) Shear force capacity
!
!  Equations (1), (2) and (3) can be combined to form a quadratic with
!  the only unknown being Vtf.  The equations may be simplified by
!  grouping the constants  viz:
!
!  Equation (1):  Let
!    s1 = (-P*1000/Ag) - (P*e*y3*1000/Ig)
!    Mcf = Mult/Vult
!    s2 =  y2*10^6/Ic
!    s3 =  Mcf*s2
!        
!  Substituting, Equation (1) becomes:  sf = s1 + Vtf*s3
!
!  Equation (2):  Let  st = 0.33SQRT(f'cg)
!    then:               SQRT((sf/2)^2 + tf^2) + sf/2 = st
!    reorganise:         SQRT((sf/2)^2 + tf^2) = st - sf/2
!    square both sides:  (sf/2)^2 + tf^2 = st^2 + (sf/2)^2 - 2*st*sf/2
!    simplify:           tf^2 = st^2 - st*sf
!    reorganise:         tf^2 + st*sf - st^2 = 0.
!
!  Equation (3):  If  Qlbwf = Qfw/(Ic*Bw2)  then  tf = Vtf*Qlbwf
!
!  Substitute for tf and sf in equation (2):
!               (Vtf*Qlbwf)^2 + st*(s1 + Vtf*s3) - st^2 = 0
!  Expand: Qlbwf^2 * Vtf^2 + st*s3*Vtf + st*s1 - st^2 = 0
!
!  To simplify the above expression let: 
!        a = Qlbwf^2
!        b = st*s3 = st*Mcf*s2
!        c = st*s1 - st*st
!
!  Then:   a*Vtf*Vtf + b*Vtf + c = 0  and the solution for Vtf is:
!             Vtf = (-b +- SQRT(b*b - 4*a*c))/2*a
!
!  Setting the determinant of this quadratic equation to:
!        f = b*b - 4*a*c
!  there is no solution if 'f' is negative.  If 'f' is positive then:
!      Vtf = (-b + SQRT(f))/2*a

y2 = Yc - Dwf
y3 = Yb - Dwf

Qlbwf = Qfw/(Ic*Bw2)
Mcf = Mult/Vult

s1 = -P*1000/Ag + P*e*1000*y3/Ig        ! Convert kN to N
s2 = 10^6*y2/Ic                         ! Convert kN.m to N.mm

a = Qlbwf*Qlbwf*1E6
b = st*Mcf*s2
c = st*s1 - st*st

f = (b*b - 4*a*c)

IF (f < 0.) THEN
    s1f = 0.
    Vtf = 0.
    Mtf = 999999.
ELSE
    Vtf = (-b + f^0.5 )/(2*a)                             
    sf = -P*1000/Ag + P*e*1000*y3/Ig + Vtf*Mcf*y2*1E6/Ic
    tf = 1000*Vtf*Qlbwf
    WorkVar1 = (0.5*sf)
    WorkVar2 = WorkVar1*WorkVar1
    WorkVar3 = WorkVar2 + tf*tf
    WorkVar4 = WorkVar3^0.5
    s1f = WorkVar4 + 0.5*sf
    Mtf = Mult*Vtf/Vult          ! Moment capacity corresponding to shear capacity (Vtf)
END IF

!   Section is cracked if Mtf > cracking moment Mcr

IF (Mtf > Mcr) THEN
    Vucf = 0.
    VucfNte$ = Cracked
ELSE
    Vucf = Vtf
    VucfNte$ = OK
END IF


!   SHEAR DESIGN - Determine Vuc
!
!   As per AS5100.5-2004, Sectn 8.2.7.2 and 
!   Warner et al.,'Concrete Structures' 1998, Sectn 12.4

WorkVar1 = MAX (Vucs, Vucc, Vucf)

WorkVar2 = Vucs
IF (Vucs = 0.) THEN Workvar2 = WorkVar1
WorkVar3 = Vucc
IF (Vucc = 0.) THEN WorkVar3 = WorkVar1
WorkVar4 = Vucf
IF (Vucf = 0.) THEN WorkVar4 = WorkVar1

Vuc = MIN (WorkVar2, WorkVar3, WorkVar4)

!  If section is cracked in flexure, only flexural-shear 
!  cracking Vuc applies (i.e. Vucs)

IF (Mult > Mcr) THEN Vuc = Vucs


!7e  SHEAR DESIGN - Design of web shear reinforcement

IF (Os = 0.) THEN Os = 0.7         ! HTML report can`t handle `ø` symbol

Vumin = (Vuc + 0.6*Bw2*do/1000)
øVumin = Os*Vumin
Vumax = 0.2*f`cg*Bw2*do/1000 + Pv
øVumax = Os*Vumax
OVumin = øVumin                    ! HTML report can`t handle `ø` symbol
OVumax = øVumax                    ! HTML report can`t handle `ø` symbol

IF (Sv = 0.) Sv = Svmax

Svr1 = 0.
Svr2 = 0.
Asvsr = 0.
Asvsr1 = 0.
Asvsr2 = 0.

IF (øVumin > Vult) THEN
   Ov = 30.
   Svr = Svmax                     ! Required shear R/F spacing
   Vus = 0.
   AsNote1$ = Asv-min
   AsNote2$ = N/A
ELSE
   Ov = 0.
   Svr = 0.
   Vus = (Vult - 0.7*Vuc)/0.7
   AsNote1$ = N/A
   AsNote2$ = Asv Reqd
END IF

Ov1 = 30.
Ov2 = 30 + 15*(Vult - øVumin)/ (øVumax - øVumin)

IF (Ov2 > Ov1) THEN Ov1 = Ov2
IF (Ov1 > 45) THEN Ov1 = 45

IF (øVumin > Vult) THEN Ov1 = 0
IF (øVumin > Vult) THEN Ov2 = 0

IF (Dsr = 0) THEN Dsr = 16

Asv = 2*3.14*Dsr*Dsr/4

IF (fsysr = 0.) THEN fsysr = 250

IF (øVumin > Vult) THEN
    Ovc = Ov
    Asvsr1 = 0.35*Bw2/fsysr       ! Required area of shear R/F per unit length
    Svr1 = Asv/Asvsr1
ELSE
    Ovc = Ov2
    WorkVar1 = Ov1*3.14159/180.   ! Convert to radians [for TAN function]
    WorkVar2 = TAN(WorkVar1)
    WorkVar3 =  ABS(WorkVar2)
    Asvsr2 = Vus*1000*WorkVar3/(do*fsysr)
    Svr2 = Asv/Asvsr2
END IF

Asvsr = Asvsr1
IF (Asvsr < Asvsr2) THEN Asvsr = Asvsr2
IF (Svr2 > 300) THEN Svr2 = Svmax

IF (øVumin > Vult) THEN
    Svr = Svmax
ELSE
    Svr = ABS(Svr2)
END IF

Asvs = Asv/Sv


!8a  DEFORMATIONS - Axial shortening and transfer hog

!    Deflection & hog calculations are only valid at mid-span

IF (Lg = 0.) THEN Lg = Span
IF (Ptm = 0.) Ptm = Pt

Xe = Ptm*Lg*1E6/(Ag*Egmt)        ! Elastic chortening
u3 = k1s*850
Xs = u3*Lg/1000                  ! Shortening due to shrinkage

øccsc = 2*k2sc*k3sc
u4 = fcscgs*øccsc*1E6/Egmt
Xc = -u4*Lg/1000                 ! Shortening due to creep

Xa = Xe + Xs + Xc                ! Total axial shortening
Occsc = øccsc                    ! HTML reports can't handle ø symbol

Dsw = 40.*Mswx*Lg*Lg*1E12/(384.*Egmt*Ig)
Hps = -Ptm*1000*e*(Lg*1000)*(Lg*1000)/(8.*Egmt*Ig)
Htr = Dsw + Hps

!8b DEFORMATIONS - Deflections

Ddeck = 40. * Mslabx * Lg*Lg * 1E12 / (384. * Eg * Ig)
Dsdl  = 40. * Msdlx * Lg*Lg * 1E12 / (384. * Eg * Ic)
Ddl = Ddeck + Dsdl
Dratio = Lg*1000/Dll
Dperm = Lg*1000/Drperm

!8c DEFORMATIONS - Hog due to creep

WorkVar5 = Ti                    ! Number of days to installation of girder
th1 = 2*Ag/Per
Fratio3 = f`cmi/f`cg

øcc1 = øccb1*k2c*k3c
Occ1 = øcc1                      ! HTML report can`t handle ø symbol
Occb1= øccb1                     ! HTML report can`t handle ø symbol
Ecc1 = fcgs * øcc1 * 1E6 / Egmi
Pc1 = Ecc1 * Ep * Ap / 1E9

!8d DEFORMATIONS - Girder hog at installation

WorkVar5 = Ty                    ! Estimated life of girder after installation
Hgic = øcc1*Htr
Hgips = -(Ptm + Pc1)*e*Lg*Lg*1E9/(8*Egmi*Ig)
Hgi = Hgic + Hgips + Dsw

!8e DEFORMATIONS - Final girder hog after installation

th2 = 2 * At / (Gp + 0.5 * Vp)
Fratio4 = f`cg / f`cg
øcc3 = øccb3 * (k2f-k2d) * k3f
Occ3 = øcc3                      ! HTML report can`t handle ø symbol
Occb3 = øccb3                    ! HTML report can`t handle ø symbol

Hel = Hgips + Dsw + Ddl
Hgfps = -P*e*Lg*Lg*1E9/(8*Egmt*Ig)
Hgfc = øcc3*Hel
Hgf = Hgfps + Hgfc + Dsw + Ddl


!9   SHEAR & TORSION CHECK - Design information to check strength of beam in torsion
!
!  The design check is for:
!     (a) Combined torsion & shear (re sufficiency of web steel)
!     (b) Combined torsion & bending (re sufficiency of longitudinal steel)
!
!  Two load cases are considered - Tult with correpsonding V and M;  and
!                                - Vult with corresponding T and M
!
!  References: (1) Section 8.3 AS5100.5-2004 
!              (2) "Concrete Structures" Warner et al. 1998 
!              (3) "Design of Prestressed Concrete", Ch 5, Gilbert & Mickleborough, 1990
!              (4) Clause 8.3.3 AS5100.5 (to check for Web Crushing failure)
!
!  Value of Jt may be assumed to be one of the following (Clause 8.3.3(3) AS5100.5):
!      0.4b2D      for solid rectangles
!      0.4SUM(x2y) for solid T or I sections
!      Am*2*bw     for thin wall hollow sections

IF (Mvcorr = 0) THEN Mvcorr = Mult

DTRflag = 0         ! Used to indicate if torsional reinforcement is inadequate


!9a  Calculate ultimate torsional strength limited by web crushing

Jt = Am*2*bw
Tumax = 0.2*f`cg*Jt/1000000         ! Eq 8.3.3(2) AS5100.5

!  To prevent web crushing beams must be proportioned so that ratio Rtv < 1 
!  (Eq 8.3.3(1) AS5100.5)   Vumax is calculated in the Shear Check section.

IF (Ot = 0.) THEN Ot = 0.7

Rtvt = Tult  / (Ot * Tumax) + Vtcorr / (Os * Vumax )
Rtvv = Tvcorr / (Ot * Tumax) + Vult  / (Os * Vumax )

Rtv = MAX(Rtvt,Rtvv)
Rtv = MIN(Rtv,2.0)                  ! Limit ratio to a maximum of 2.0

Tnote1$ = ^Web will not crush
Tnote18$ = Rtv must be less than 1.0 (Eq 8.3.3(1))

IF (Rtv > 1.) THEN
   Tnote1$ = ?Web will crush
   DTRflag = 1
END IF 


!9b  Check if torsional reinforcement is needed  (Clause 8.3.4 AS5100.5)
!
!    Calculate ultimate strength in pure torsion (Tuc)

Scp = P*1000/Ag                                                
WorkVar1 = (1 + 10 * Scp / f`cg) ^ 0.5
Tuc = 0.3 * f`cg ^ 0.5* Jt * WorkVar1/1000000   ! Eq 8.3.5(1) AS5100.5

!  Check if torsional reinforcement is needed  (Clause 8.3.4(a) AS5100.5)

Tucmax = 0.25*Ot*Tuc 
CloseTie = 0

Tnote2$ = ^Torsion reinforcement not required

IF (Tucmax < Tult) THEN
   Tnote2$ = ?Torsional reinforcement (closed ties) required
   CloseTie = 1
   DTRflag = 1
END IF

Rtv1tmax = Tult/(Ot*Tuc) + Vtcorr/(Os*Vuc)
Rtv1tmax = MIN(Rtv1tmax,2.0)                 ! Limit ratio to a maximum of 2.0
Rtv1vmax = Tvcorr/(Ot*Tuc) + Vult/(Os*Vuc)
Rtv1vmax = MIN(Rtv1vmax,2.0)                 ! Limit ratio to a maximum of 2.0
Rtv1 = MAX(Rtv1tmax,Rtv1vmax)

Tnote20$ = Note: Tucmax must be > Tult and Rtv1 < 0.5

IF (Rtv1 > 0.5) THEN
   Tnote2$ = ?Torsional reinforcement (closed ties) required
   CloseTie = 1
   DTRflag = 1
END IF

Tnote3$ = ^Closed ties are not required
IF (CloseTie > 0) THEN Tnote3$ = ?Closed ties are required


!9c  Check for combined torsion and shear

!   (a) Determine required Asw/s ratio
  
Vumin = Vuc + 0.6*bv*do/1000         ! Eq 8.2.9 AS5100
Vdesign = Vtcorr/Os                                      

Tnote4$ = ^Vumin > Vdesign - No extra shear RF needed

IF (Vumin < Vdesign) THEN 
   Tnote4$ = ?Vumin < Vdesign - Extra shear RF needed
   DTRflag = 1
END IF

Asvmin$ = Use Asvmin

IF (Vumin < Vdesign) Asvmin$ = Asvmin insuficient

!   Calculate thetav and thetat for both Tult/Vtcorr and Vult/Tvcorr

thetav1 = 30 + 15*(Vtcorr - Os*Vumin)/(Os*Vumax - Os*Vumin)
thetav = thetav1
IF (thetav < 30) THEN thetav = 30
IF (thetav > 45) THEN thetav = 45
thetavc = thetav                                        ! thetav due to corresp Vt

thetav1 = 30 + 15*(Vult - Os*Vumin)/(Os*Vumax - Os*Vumin)
thetav = thetav1
IF (thetav < 30) THEN thetav = 30
IF (thetav > 45) THEN thetav = 45
thetavm = thetav                                        ! thetav due to max Vt

thetat1 = 30 + 15*(Tult - Ot*Tuc)/(Ot*Tumax - Os*Tuc)   ! Warner, Eq 13.17
thetat = thetat1
IF (thetat < 30) THEN thetat = 30
IF (thetat > 45) THEN thetat = 45
thetatm = thetat                     ! thetat due to max T

thetat1 = 30 + 15*(Tvcorr - Ot*Tuc)/(Ot*Tumax - Os*Tuc)  ! Warner, Eq 13.17
thetat = thetat1
IF (thetat < 30) THEN thetat = 30
IF (thetat > 45) THEN thetat = 45
thetatc = thetat                                         ! Use Thetat due to corresp T

!   Atg = area of polygon with vertices at the centre of the bars

nlegsh = 2 * nclty

WorkVar8 = thetatm * 3.14159 / 180.                      ! Use Theta corresp to max Torsion
WorkVar8 = TAN(WorkVar8)
Aswstm = Tult*1000000 / (Ot*2*fsysr*Atg/WorkVar8)        ! Warner, Eq 13.21

WorkVar8 = thetatc * 3.14159 / 180.                      ! Use Theta corresp to correp Torsion
WorkVar8 = TAN(WorkVar8)
Aswstc = Tvcorr*1000000 / (Ot*2*fsysr*Atg/WorkVar8)      ! Warner, Eq 13.21

Aswst = MAX(Aswstm, Aswstc)

WorkVar8 = thetavm * 3.14159 / 180.                      ! Use Theta corresp to max Vt
WorkVar8 = TAN(WorkVar8)
Aswsvm = Vult*1000 / (Os*nlegsh*fsysr*do/WorkVar8)       ! Warner, Eq 13.23

WorkVar8 = thetavc * 3.14159 / 180.                      ! Use Theta corresp to corresp Vt
WorkVar8 = TAN(WorkVar8)
Aswsvc = Vtcorr*1000 / (Os*nlegsh*fsysr*do/WorkVar8)     ! Warner, Eq 13.23

Aswsv = MAX(Aswsvm, Aswsvc)

Aswons = Aswst + Aswsv                                   ! Min total required R/F (Warner, Eq 13.26)

IF (Aswons > 1E9) THEN Aswons = 1E9

IF (Acltie = 0) THEN 
   Acltie = 3.14*TieDiam^2 / 4                        
END IF

Asws = Acltie/Sdesign                            
 
IF (Asws > 1E9) THEN Asws = 1E9

!   (b) Check minimum torsional reinforcement  (Clause 8.3.7 AS5100.5)

Aswsmin1 = 0.35*bv / (nlegsh*fsysr)                      ! Warner,  Eq 13.32
Aswsmin2 = 0.2*yt1 / fsysr                               ! AS5100 - Clause 8.3.7

IF (Aswsmin1 > Asws) THEN
   Tnote6$ = ?Use minimum torsional reinforcement
   DTRflag = 1
ELSE
   Tnote6$ = ^Satisfies minimum ligs test (Warner, Eq 13.32)
END IF

IF (Aswsmin2 > Asws) THEN 
   Tnote7$ = ?Does not satisfy minimum ligs test
   DTRflag = 1
ELSE
   Tnote7$ = ^Torsional RF OK (Clause 8.3.7 AS5100)
END IF

WorkVar8 = thetatm * 3.14159 / 180. 
WorkVar8 = TAN(WorkVar8)
Tus1 = 1E-6*2*fsysr*Asws*Atg / WorkVar8                  ! Eq 8.3.5(2) AS5100

IF (Tus1 > 1E9) THEN Tus1 = 0

WorkVar8 = thetatc * 3.14159 / 180.
WorkVar8 = TAN(WorkVar8)
Tus2 = 1E-6*2*fsysr*Asws*Atg / WorkVar8                   ! Eq 8.3.5(2) AS5100

IF (Tus2 > 1E9) THEN Tus2 = 0

Tus = MAX(Tus1, Tus2)

IF (Tus > Tuc) THEN 
   Tnote8$ = ^Asw/s satisfies Tus > Tuc
ELSE
   Tnote8$ = ?May need to increase Asw/s to satisfy Tus>Tuc   
   DTRflag = 1
END IF

!   Check if ratio T*/0t*Tus + V*/(0s*Vs+Pv) is less than 1

WorkVar8 = thetavm*3.14159/180.
WorkVar8 = TAN(WorkVar8)
WorkVar8 = WorkVar8 * 1000.
Vusteel1 = nlegsh*Asws*fsysr*do/WorkVar8           ! Eq 8.2.10 AS5100

Rtv2t = Tult/(Ot*Tus) + Vtcorr/(Os*Vusteel1 + Pv)  ! Check max Torsion

WorkVar8 = thetavc*3.14159/180.
WorkVar8 = TAN(WorkVar8)
WorkVar8 = WorkVar8 * 1000.
Vusteel2 = nlegsh*Asws*fsysr*do/WorkVar8           ! Eq 8.2.10 AS5100

Rtv2c = Tvcorr/(Ot*Tus) + Vult/(Os*Vusteel2 + Pv)  ! Check corrresp torsion
Rtv2 = MAX(Rtv2c, Rtv2t)
Vusteel = MAX(Vusteel1, Vusteel2)                  ! Ult shear strength of stirrups

IF (Rtv2 < 1) THEN				   ! Cl 8.3.4(b) AS5100
   Tnote9$ = ^Combined shear/torsion ratio is OK
ELSE
   Tnote9$ = ?Combined shear/torsion ratio > 1. Adjust Asws  
   DTflag = 1
END IF

Tnote22$ = Rtv2 must be < 1. If not, adjust Asws

!   Design for closed ties for combined shear and torsion

Stiecalc = Acltie / Asws                          ! Raswsv is input by user
Sclostie = Stiecalc

Sclmax1 = 0.12 * ut
Sclmax  = Sclmax1

IF (Sclmax > 300) THEN Sclmax = 300

IF (Sclostie > Sclmax) THEN                        ! Need to reduce spacing to maximum allowed
   Sclostie = Sclmax
END IF

Raswsy = Acltie / Sclostie

IF (DTRflag = 0) THEN
   DTRnote$ = ^Torsional R/F is OK
ELSE
   DTRnote$ = ?Torsional R/F may be inadequate
END IF


!9d  Design for additional longitudinal reinforcement  (Clause 8.3.6 AS5100)
!
!    (a) Flexural tensile zone
!
!    Ftadd  =  Design tensile forces to be taken in addition to any design tensile forces due to flexure
!    Ttcapbs = Tensile capacity of bars in tension zone including both strand and passive RF
!    Ftmcorr = Tensile force due to moment corresponding to design Torsion
!    Ftbdif  =  Residual capacity of bars in tension zone at the prescribed corresponding moment 
!
!    Ftadd = 0.5*Aswst*ut*(1/TAN(thetat/180*3.14))^2*fsysr/1000  ........Eqn 8.3.6(1)
!    Atadd = 0.5*Aswst*ut*(1/TAN(thetat/180*3.14))^2             ........derived from 8.3.6(1)
!
!  Where  Aswst = Area of torsional R/F;  thetat = corresponding angle of torsional RF
!  Need to check the case for both max T as well as corresp T

WorkVar8 = thetatm * 3.14159 / 180.
WorkVar8 = TAN(WorkVar8)
WorkVar8 = 1.0 / WorkVar8
WorkVar8 = WorkVar8 ^ 2

ALRflag = 0

Ftadd1 = 0.5* Aswstm * ut * WorkVar8 * fsysr / 1000
Atadd1 = 0.5* Aswstm * ut * WorkVar8

WorkVar8 = thetatc * 3.14159 / 180.
WorkVar8 = TAN(WorkVar8)
WorkVar8 = 1.0 / WorkVar8
WorkVar8 = WorkVar8 ^ 2

Ftadd2 = 0.5* Aswstc * ut * WorkVar8 * fsysr / 1000
Atadd2 = 0.5* Aswstc * ut * WorkVar8

Ftadd = MAX(Ftadd1, Ftadd2)
Atadd = MAX(Atadd1, Atadd2)

!   Note:  Ftcapbs = tensile capacity of bars & strands in tensile zone

Ftcapbs = 0.

LOOP i = 1,8
   IF (Ybarri < Yc) Ftcapbs = Ftcapbs + Fbari
END LOOP

LOOP i = 1,20
   IF (Yrfi < Yc) Ftcapbs = Ftcapbs + Fbrfi
END LOOP

Ftmcorr = Mtcorr*1000 / (0.8*0.9*do)                   ! Tensile force due to corresp M
Ftbdif = Ftcapbs - Ftmcorr

LongReo = 0
Aaddltr = 0                                            ! Additional long'l torsional RF

Tnote21$ = Ftbdif must be > Ftadd and Abotreo > Aadltr

IF (Ftbdif < Ftadd) THEN
   Tnote10$ = ?Additional longitudinal RF is required
   LongReo = 1
   Aadltr = (Ftadd - Ftbdif) * 1000 / fsysr
   ALRflag = ALRflag + 1
ELSE
   Tnote10$ = ^No extra longitudinal RF is required
END IF

!   Calculate area of passive RF in tension zone

Abotreo = 3.14159*(BarDtens^2/4)*Ntbars                ! NOTE: User sets BarDtens & Ntbars
WorkVar1 = 0

IF( Abotreo > Aadltr) THEN
   Tnote11$ = ^Longitudinal RF in flexural tensile zone is OK
ELSE
   Tnote11$ = ?Need more longitudinal RF in flexural tensile zone
   ALRflag = ALRflag + 1
   WorkVar1 = 1
END IF

IF( Abotreo = Aadltr) THEN
   Tnote11$ = ^Longitudinal RF in flexural tensile zone is OK
   ALRflag = ALRflag - WorkVar1
END IF

!   Check flexural compression zone

!   Fcadd    = Additional force required in the compression zone
!   Ftstrand = Tensile capacity of strand in compression zone
!
!   Usually strands in the flange area are located within the flange/web region
!   and hence enclosed by stirups.  They can therefore be considered to be able 
!   to contribute to the torsional strength and have consequently been included 
!   in the calculations.

Fcadd0 = Ftadd - Mtcorr*1000000/(0.8*fsysr*0.85*do)
Fcadd  = Fcadd0

IF (Fcadd < 0) Fcadd = 0

!  Note Ftstrand = Tensile capacity of strands in compression zone 
!  (must be calculated here)

IF (Fcadd < Ftstrand) THEN
   Aaddlcr = 0.
ELSE
   Aaddlcr = (Fcadd - Ftstrand)*1000000 / (fsysr*Wtf)
END IF

IF (Scbars = 0) THEN
   Atopreo = 0.
ELSE
   Atopreo = 3.14159*(BarDcomp^2/4)*(1000/Scbars)                 ! User sets Scbars and BarDcomp
END IF

WorkVar1 = 0

IF (Atopreo > Aaddlcr) THEN
   Tnote12$ = ^Long. RF in flexural compression zone OK
ELSE
   Tnote12$ = ?Need more longitudinal RF in flexural compression zone
   ALRflag = ALRflag + 1
   WorkVar1 = 1
END IF

IF (Atopreo = Aaddlcr) THEN
   Tnote12$ = ^Long. RF in flexural compression zone OK
   ALRflag = ALRflag - WorkVar1
END IF

IF (ALRflag = 0) THEN
   ALRnote$ = ^Additional longitudinal R/F is not required
ELSE
   ALRnote$ = ?Additional longitudinal R/F is required
END IF