Hurmet module with functions for strength of steel sections per AISC 360-16

steelStrengthPerAISC360-16.hrms

# Hurmet functions to find strength of steel beams and columns per AISC 360-16.
# Copyright 2020 - 2022 Ron Kok
# Released under terms of the MIT License, https://opensource.org/licenses/MIT

function Ps(section, Fy, kLx, kLy)
  # LRFD axial strength of a steel member, per AISC 360-16 sections B & E
  Pcr = criticalAxialStrength(section, Fy, kLx, kLy)
  return 0.9 × Pcr
end

function Pas(section, Fy, kLx, kLy)
  # Service level allowable axial strength of a steel member, per AISC 360-16 sections B & E
  Pcr = criticalAxialStrength(section, Fy, kLx, kLy)
  return Pcr / 1.67
end

function Ms(section, Fy, Lb, Cb, axis)
  # LRFD bending strength of a steel member, per AISC 360-16 sections B & F
  if axis = undefined
    axis = "x"
  end
  if Cb = undefined
    Cb = 1
  end
  Mcr = criticalBendingStrength(section, Fy, Lb, Cb, axis)
  return 0.9 × Mcr
end

function Mas(section, Fy, Lb, Cb, axis)
  # Service level allowable bendingstrength of a steel member, per AISC 360-16 sections B & F
  if axis = undefined
    axis = "x"
  end
  if Cb = undefined
    Cb = 1
  end
  Mcr = criticalBendingStrength(section, Fy, Lb, Cb, axis)
  return 0.9 × Mcr
end

private function criticalAxialStrength(section, Fy, kLx, kLy)
  # Find critical axial strength

  E = 29000 'ksi'  # steel modulus of elasticity

  # What kind of section? e.g. "I", "channel", "HSS", etc.
  shape = shapeOf(section)
  if shape = "L"
    raise "Error. This function doesn’t do single angles."
  end
  if shape = "N/A"
    raise "Error. Unrecognized section: " & name
  end

  if shape = "round"
    A, rx = section["A", "r"]
    ry = rx
  else
    A, rx, ry = section["A", "rx", "ry"]
  end

  if kLy = undefined
    if shape = "round" or section.Ix = section.Iy
      kKy = kLx
    else
      raise "Error. Undefined kLy"
    end
  end

  # Get critical compressive stress, Fcr, per Section E3
  S_r = max(kLx / rx, kLy / ry)
  Fe = π² E / S_r                   # Eqn E3-4, Euler buckling
  Fcr = {
    0.658^(Fy/Fe) × Fy  if S_r ≤ 4.71  √(E / Fy) ;  # Eqn E3-2
    0.877 Fe            otherwise                   # Eqn E3-3
  }  

  # Get the effective area after checking for slender elements per AISC Table B4-1a
  A = {
    AeffRound(section, Fy, E)    if shape = "round" ;
    AeffHSS(section, Fy, E, Fcr)   if shape = "HSS" ;
    AeffAngles(section, Fy, E, Fcr)  if shape = "2L" ;
    Aeff(section, Fy, E, Fcr, shape) otherwise    # I, channel, or tee
  }

  # Check the Section E3 critical stress against Section E4
  if A = section.A and shape ∉ ["round", "HSS"]
    # Section E4 applies only to elements without slender elements
    Cw, J, Ix, Iy = section["Cw", "J", "Ix", "Iy"]
    G = 12000 'ksi'
    if shape = "I"
      # The section is doubly symmetric.
      Lcz = max(kLx, kLy)  # conservative
      Fe = ((π² E Cw) / Lcz² + G J) × 1 / (Ix + Iy)  # Eqn E4-2
    else if shape ∈ ["channel", "tee", "2L"]
      # Singly symmetric.
      rx, ry, y = section["rx", "ry", "y"]
      xo = { section.eo + section.x if shape = "channel"; 0 'in' otherwise }
      yo = {
        0 'in'       if shape = "channel" ;
        y - section.t / 2  otherwise
      }
      ro = xo² + yo² + (Ix + Iy) / A
      Fey = π² E / (kLy / ry)²
      Fez = {
        ((π² E Cw)/min(rx, ry) + G J) × 1/(A ro²) if shape = "channel" ;
        G J / (A ro²)               otherwise
      }
      H = 1 - (xo² + yo²) / ro²
      Fe = ((Fey + Fez)/(2 H))(1 -√(1 - (4 Fey Fez H)/(Fey + Fez)²))  # Eqn E4-3
    end
    Fcr = min(Fcr, Fe)
  end
  return A × Fcr
end

private function shapeOf(section)
  name = section.name
  s1 = name[1]
  s2 = name[1..2]
  return {
    "round"   if s2 = "Pi" or (s2 = "HS" and "x" ∉ name) ;
    "HSS"     if s2 = "HS" ;
    "tee"     if s2 ∈ ["WT", "MT", "ST"] ;
    "channel" if s1 = "C" or s2 = "MC" ;
    "I"       if s1 ∈ ["W", "M", "S"] or s2 = "HP" ;
    "2L"      if s2 = "2L" ;
    "N/A"     otherwise
  }
end

private function AeffRound(section, Fy, E)
  # Get effective area for round HSS and pipes, taking into account slender walls.
  A, D, t = section["A", "OD", "tdes"]
  if D / t > 0.45 × E / Fy 
    raise "AISC doesn’t give a value for a D/t ratio this big."
  else if D / t < 0.11 × E / Fy
    return A                                  # Table B4.1a case 9
  else
    return A × min(1, 0.38 E / (Fy × D/t) + 2/3)  # Eqn E7-19
  end
end

private function AeffHSS(section, Fy, E, Fcr)
  # Get effective area of rectangular or square HSS section
  A, b, h, t = section["A", "b", "h", "tdes"]
  λf = b / t             # width to thickness ratio of flat part of "flange"
  λw = h / t             # ditto for "web"
  λr = 1.40 × √(E / Fy)  # limiting ratio per Table B4.1a

  if λf < λr and λw < λr
    return A
  end
  
  sqrFyFc = √(Fy/Fcr)
  if λf > sqrFyFc × λr
    Fel = (1.38 λr / λf)² Fy                      # Eqn E7-5
    be = b (1 - 0.2 √(Fel / Fcr)) × √(Fel / Fcr)  # Eqn E7-3
    A = A - 2 (b - be) t
  else if λw > sqrFyFc × λr
    Fel = (1.38 λr / λw)² Fy
    he = h (1 - 0.2 √(Fel / Fcr)) × √(Fel / Fcr)
    A = A - 2 (h - he) t
  end
  return A
end

private function Aeff(section, Fy, E, Fcr, shape)
  # Get effective area for I sections, channels, and tees
  A, d, bf, tf, tw, kdes = section["A", "d", "bf", "tf", "tw", "kdes"]
  b = { bf if shape = "channel"; bf /2 otherwise }  # per section B4.1a(a)
  h = { d if shape = "tee"; d - 2 kdes otherwise }

  λf = b / tf              # width to thickness ratio of flange
  λf_r = 0.56 × √(E / Fy)  # limiting ratio per Table B4.1a
  λw = h / tw        # web
  λw_r = { λf_r if shape = "tee"; 1.49 × √(E / Fy) otherwise }

  if λf < λf_r and λw < λw_r
    return A
  end
  
  sqrFyFc = √(Fy/Fcr)
  if λf > sqrFyFc × λf_r
    Fel = (1.49 λf_r / λf)² Fy                     # Eqn E7-5
    be = b (1 - 0.22 √(Fel / Fcr)) × √(Fel / Fcr)  # Eqn E7-3
    A = A - 2 (b - be) tf
  end

  if λw > sqrFyFc × λw_r
    if shape = "tee"
      Fel = (1.49 λf_r / λf)² Fy
      de = d (1 - 0.22 √(Fel / Fcr)) × √(Fel / Fcr)
    else
      Fel = (1.31 λf_r / λf)² Fy
      de = d (1 - 0.18 √(Fel / Fcr)) × √(Fel / Fcr)
    end
    A = A - (d - de) tw
  end
  return A
end

private function AeffAngles(section, Fy, E, Fcr)
  # Get effective area for double angles
  A, b, d, t = section["A", "b", "d", "t"]
  λf = b / t             # width to thickness ratio of flat part of "flange"
  λw = h / t             # ditto for "web"
  λr = 0.45 × √(E / Fy)  # limiting ratio per Table B4.1a

  if λf < λr and λw < λr
    return A
  end
  
  sqrFyFc = √(Fy/Fcr)
  if λf > sqrFyFc × λr
    Fel = (1.49 λr / λf)² Fy             # Eqn E7-5
    be = b (1 - 0.22 √(Fel / Fcr)) × √(Fel / Fcr)  # Eqn E7-3
    A = A - 2 (b - be) t
  end
  if λw > sqrFyFc × λr
    Fel = (1.49 λr / λw)² Fy
    he = h (1 - 0.22 √(Fel / Fcr)) × √(Fel / Fcr)
    A = A - 2 (h - he) t
  end

  return A
end

private function criticalBendingStrength(section, Fy, Lb, Cb, axis)
  # Find critical bending strength

  E = 29000 'ksi'  # steel modulus of elasticity
  S_EFy = √(E / Fy)

  # What kind of section?
  shape = shapeOf(section)
  if shape = "L"
    raise "Error. This function doesn’t do single angles."
  end
  if shape = "N/A"
    raise "Error. Unrecognized section: " & name
  end
  if shape ∈ ["tee", "2L"] and axis ∉ "xX"
    raise "Error. This function does not do tees or double angles bent on their y-axis."
  end

  if shape = "round"
    # AISC section F8
    D, S, Z, t = section["OD", "S", "Z", "tdes"]
    if D / t > 0.45 E / Fy
      raise "AISC doesn't give values for D/t this large."
    end
    Mp = Fy × Z
    if D / t < 0.07 E / Fy
      return Mp
    else if D / t < 0.31 E / Fy
      Mb = (0.021 E / (D / t) + Fy) × S      # EQ F8-2
    else
      Fcr = 0.33 E / (D / t)       # EQ F8-4
      Mb = Fcr × S                 # EQ F8-3
    end
    return min(Mp, Mb)
  end if

  if shape = "HSS"
    # Section F7 for rectangular HSS
    if "X" ∈ name
      if axis ∈ "xX"
        H, S, Z, b, t = section["Ht", "Sx", "Zx", "b", "tdes"]
      else
        H, S, Z, b, t = section["B", "Sy", "Zy", "h", "tdes"]
      end
    else
      H, S, Z, b, t = section["Ht", "S", "Z", "b", "tdes"]
    end
    Mp = Fy Z                 # Eqn F7-1
    λ = b / t                 # "flange" slenderness ratio
    if λ ≤ 1.12 × √(E / Fy)   # Table B4-1.b case 17, compact limit
      return Mp
    end
    if λ ≤ 1.4 × √(E / Fy)    # ditto, slender limit
      return min(Mp, Mp - (Mp - Fy S) (3.57 b/t √(Fy/E) - 4.0))
    end
    be = min(b, 1.92 t √(E/Fy) (1 - 0.38/(b/t) √(E/Fy)))
    ΔI = 2((b - be)t³)/12 + 2(b - be)t(1/2(H-t/2))²
    S = S - ΔI / (H/2)
    return Fy S
  end

  if shape = "tee"
    # Section F9
    d, Sx, t = section["d", "Sx", "t"]
    # This function assumes that stems are in compression, which is conservative.
    Mp = Fy Sx
    # Section F9.4
    Fcr = {
      Fy                           if d/t ≤ 0.84 S_EFy ;
      (1.43 - 0.515 d/t S_EFy) Fy  if d/t ≤ 1.52 S_EFy ;
      (1.52 E) / (d/t)²            otherwise
    }
    return Fcr Sx
  end

  if shape = "2L"
    #Section F9
    b, Sx, t = section["d", "Sx", "t"]
    # This function assumes that stems are in compression, which is conservative.
    # Section F9.4 & F10.3 & Table B4.1b
    Sc = 0.8 Sx
    if b/t ≤ 0.84 S_EFy
      return Fy Sc
    else if b/t ≤ 1.52 S_EFy
      return Fy Sc (2.43 - 1.72 × (b/t) × S_EFy)  # Eqn F10-6
    else
      return (0.71 E) / (b/t)² × Sc               # Eqn F10-7 & F10-8
    end
  end

  # The rest of this function deals with I-shapes and channels
  
  λpf = 0.38 × √(E / Fy)    # Table B4-1.b case 10, compact limit for flanges
  λrf = 0.56 × √(E / Fy)    # ditto, slender limit

  if axis ∉ "xX"
    # Minor axis. Use section F6
    Sy, Zy, bf, tf = section["Sy", "Zy", "bf", "tf"]
    Mp = min(Fy Zy, 1.6 Fy Sy)  # Eqn F6-1
    b = bf / 2          # ref Table  B4-1
    λ = b / tf
    if λ ≤ λpf     # compact flange
      return Mp
    end
    if λ ≤ λrf     # non-compact flange
      return Mp - (Mp - 0.7 Fy Sy) ((λ - λpf) / (λrf - λpf))   # Eqn F6-2
    end
    # slender flange
    Fcr = 0.69 E / (b / tf)²
    return Fcr Sy
  end

  G = 12000 'ksi'
  Sx, Zx, Iy, d, bf, tw, tf, k, J, Cw, rt = section["Sx", "Zx", "Iy", "d", "bf", "tw", "tf", "kdes", "J", "Cw", "rts"]
  b = { bf if shape = "channel"; bf / 2 otherwise }
  λf = b / tf
  h = d - 2 k
  ho = d - tf
  λw = h / tw
  λpw = 3.76 × √(E / Fy)    # Table B4-1.b case 15, compact limit
  λrw = 5.79 × √(E / Fy)    # ditto, slender limit
  Mp = Fy Zx
  
  if λf ≤ λpf and λw ≤ λpw
    # Compact section. Use section F2.
    Mcr = Cb π/Lb √(E Iy G J + ((π E)/Lb)² Iy Cw)
    return min(Mp, Mcr)
  end
  
  if shape = "channel"
    raise "This function does not handle non-compact channels. All current channels are compact per user note in AISC 360 section F2."
  end
  
  # The remaining cases are non-compact I-sections bent on their major axis.
  # This function uses AISC 360 section F5 instead of F4. So we're conservative.

  # Check compression flange yielding. F5.1
  aw = (h tw) / (bf tf)                     # Eqn F4-12
  Rpg = min(1, 1 - aw / (1200 + 300 aw) × (h/tw - 5.7 S_EFy))   # Eqn F5-6
  Mn1 = Rpg Fy Sx

  # Check lateral-torsional buckling. F5.2
  Lp = (Cb π² E) / (Lb/rt)² × √(1 + 0.078 × J / (Sx ho) × (Lb/rt)²)   # Eqn F4-7
  Lr = π rt × √(E / (0.7 Fy))
  Mn2 = {
    Mp                                                        if Lb < Lp ;
    min(Fy, Cb (Fy - (0.3 Fy)((Lb - Lp) / (Lr - Lp)))) × Sx   if Lb < Lr ;   # Eqn F5-3
    min(Fy, (Cb π² E) / (Lb/rt)²) × Sx                        otherwise ;    # Eqn F5-4
  }

  # Check compression flange local buckling. F5.3
  kc = 4 / √(h / tw)
  Mn3 = {
    Mp                                               if λf ≤ λpf ;
    (Fy - (0.3 Fy) ((λf - λpf) / (λrf - λpf))) × Sx  if λf ≤ λrf ;   # Eqn F5-8
    ((0.9 E kc) / (bf / (2 tf))²) × Sx               otherwise     # Eqn F5-9
  }
  
  return min(Mp, Mn1, Mn2, Mn3)
end