Hurmet module to create a Concrete Column Interaction Diagram

concreteColumnInteraction.txt

# ConcreteColumnInteraction
#
# By Ron Kok. Copyright 2023.
# Released under terms of the MIT license, https://mit-license.org/
#
# This is a Hurmet module that draws a concrete column interaction diagram
# for recangular columns per ACI 318.
# The first input, `column`, must be a single-row dataframe, like this example:
#
# column = 
# ``f_c′	f_y	width	depth	barSize	bars	cover	tieSize
# psi	psi	in	in			in	
# 4500	60000	24	30	#7	4x5	2	#4``
#
# To create a diagram in Hurmet.app, first create a cell that fetches this code, e.g.:
# colDiagram = import("https://gist.githubusercontent.com/ronkok/6ea53b79cd49a0ab5c6c60e3f9e8c874/raw/concreteColumnInteraction.txt") = !
#
# Then create a cell that uses the module to draw a diagram, e.g.:
# colDiagram.draw(column, Pu, Mu) = @
#
# Pu and Mu are optional arguments that contain the axial and bending forces on the column.
#
# This module is not unit-aware. The units in the example are mandatory.
# If the number of bars is evenly divisible by 4, you can write "bars" as a single number.
# Otherwise, the bar pattern must be written in the form: 𝐦x𝐧
#  where 𝐦 and 𝐧 are integers ≥ 2.
#
# This diagram is not the last word on column capacity. Slenderness must be
# addressed elsewhere. Also, ACI 318 has several prescriptive requirements that
# must be met in order for this diagram to be valid. Such as: maximum bar
# spacing, minimum ρ values, and minimum cover. You
# are responsible to check those prescriptive requirements.
#
# Only qualified engineers should use this diagram.

draw(column; Pu = 0, Mu = 0)
    title "Concrete Interaction Diagram"
    frame 640, 180
    view 0, 800, 0, 225
    fontsize 12

    text [80, 215], "**Concrete Column**"
    text [102, 185], "f~c~′", "left"
    text [102, 163], "f~y~", "left"
    text [102, 141], "ACI 318-", "left"
    text [102, 119], "Width", "left"
    text [102, 97], "Depth", "left"
    text [102, 75], "Bar size", "left"
    if isnan(column.bars)
        text [102, 53], "Bar pattern", "left"
    else
        text [102, 53], "Bar qty", "left"
    end
    text [102, 31], "Clear cover", "left"
    text [102, 9], "Tie size", "left"

    if isnan(column["f_c′"]) throw "Error. fc′ must be numeric"
    if isnan(column["f_y"]) throw "Error. fy must be numeric"
    if isnan(column.width) throw "Error. Width must be numeric"
    if isnan(column.depth)
        throw "Error. Depth must be numeric"
    end

    fc′ = column["f_c′"] / 1000
    fy = column["f_y"] / 1000
    b = column.width
    H = column.depth
    barSize = column.barSize
    cover = column.cover
    tieSize = column.tieSize
    barPattern = column.bars

    rebar = ``#size	diameter	area
	in	in²
#3	0.375	0.11
#4	0.50	0.20
#5	0.625	0.31
#6	0.75	0.44
#7	0.875	0.60
#8	1.00	0.79
#9	1.125	1.00
#10	1.27	1.27
#11	1.40	1.56
#14	1.69	2.25
#18	2.26	4.00``

    barDia = rebar[barSize].diameter
    barRadius = barDia / 2
    areaPerBar = rebar[barSize].area
    tieDia = rebar[tieSize].diameter
    tieIR = 2 * tieDia  # inside radius of tie bending curve
    tieOR = 3 * tieDia
    Ag = b * H
    if isnan(barPattern)
        pos = findfirst("x", barPattern)
        numTopBars = number(barPattern[1:pos-1])
        numLayers = number(barPattern[pos+1:end])
        numBars = 2 * numTopBars + 2 * (numLayers - 2)
    else
        numBars = barPattern
        if mod(numBars, 4) ≠ 0 throw "“bars” should be divisible by 4, or written with two integers: 𝐦x𝐧."
        numTopBars = (numBars / 4) + 1
        numLayers = numTopBars
    end
    if numBars < 4  throw "ACI 318 section 10.7 requires at least 4 vertical bars."
    As = numBars * areaPerBar  # steel area
    ρ = As / Ag
    if ρ < 0.005 or ρ > 0.08
        throw "ACI 318 sections 10.6.1 and 10.3.1.2 require ρ to be between 0.5% and 8%."
    end
    # Draw the column in a box 100px wide.
    scale = 100 / max(b, H) / 0.8

    text [102, 185], string(fc′, "r2") , "right"
    text [140, 185], "ksi", "right"
    text [102, 163], string(fy, "f0"), "right"
    text [140, 163], "ksi", "right"
    text [102, 141], "19", "right"
    text [102, 119], string(b, "f1"), "right"
    text [140, 119], "in", "right"
    text [102, 97], string(H, "f1"), "right"
    text [140, 97], "in", "right"
    text [102, 75], barSize, "right"
    if isnan(barPattern)
        text [102, 53], barPattern, "right"
    else
        text [102, 53], string(numBars, "f0"), "right"
    end
    text [102, 31], string(cover, "f2"), "right"
    text [140, 31], "in", "right"
    text [102, 9], tieSize, "right"
    # Draw the column outline
    bot = 95
    rect [210, bot; 210 + scale * b, bot + scale * H]
    # Draw the ties
    fill "black"
    strokewidth 0
    rect [210 + scale * cover, bot + scale * cover;  210 + scale * (b - cover), bot + scale * (H - cover)], scale * tieOR
    fill "white"
    strokewidth 0
    rect [210 + scale * (cover + tieDia), bot + scale * (cover + tieDia);  210 + scale * (b - cover - tieDia), bot + scale * (H - cover - tieDia)], scale * tieIR
    # Draw the vertical bars
    fill "black"
    
    brp = barRadius * scale # bar radius in px
    topSpacing = 0
    sideSpacing = 0
    margin = {
        cover + tieDia + (1 - 0.7071) * tieIR - (1 - 0.7071) * barRadius + barRadius if (1 - 0.7071) * tieIR - (1 - 0.7071) * barRadius > 0;
        cover + tieDia + barRadius + 0.0625 otherwise
    }
    topSpacing = (b - 2 * margin) / (numTopBars - 1)
    sideSpacing = (H - 2 * margin) / (numLayers - 1)
    for i in 0:(numTopBars - 1)
        circle [210 + (margin + i * topSpacing) * scale, bot + margin * scale], brp
        circle [210 + (margin + i * topSpacing) * scale, bot + (H - margin) * scale], brp
    end
    for i in 1:numLayers - 2
        circle [210 + margin * scale, bot + (margin + i * sideSpacing) * scale], brp
        circle [210 + (b - margin) * scale, bot + (margin + i * sideSpacing) * scale], brp
    end
    # Write in the Ag, As, and ρ values
    strokewidth 1
    text [210, bot - 20], "A~g~ = " & string(Ag, "f1") & " in²", "right"
    text [210, bot - 42], "A~s~ = " & string(As, "f1") & " in²", "right"
    text [210, bot - 64], "ρ = " & string(ρ * 100, "f2") & "%", "right"
    if topSpacing == sideSpacing
        text [210, bot - 86], "bars at " & string(topSpacing, "f2") & "″ c/c", "right"
    else
       text [210, bot - 86], "bars at " & string(topSpacing, "f2") & "″ × " & string(sideSpacing, "f2") & "″ c/c", "right"
    end
  
    # Draw the x & y axis for the diagram
    baseX = 460
    baseY = 40
    topOfGraph = baseY + 150
    line [baseX, baseY; baseX, baseY + 180]
    line [baseX, baseY; baseX + 200, baseY]
    
    # øPn label on the y-axis
    text [baseX, 210], "ϕP~n~ (kips)", "left"
    # 0 at the origin
    text [baseX, 40], "0", "belowleft"
    # øMn label on the x-axis
    text [630, 10], "ϕM~n~ (kip·ft)"

    # Do the math to find the curve points
    Px, Mx, Mxb, Pmax, MxMax, PatMxMax, MxAtPmax = curveArray(fc′, fy, b, H, Ag, numBars, numTopBars, numLayers, areaPerBar, barRadius, barPattern, cover, tieDia, margin)
  
    needTwoCurves = false
    Py = 0
    My = 0
    Myb = 0
    MyMax = 0
    MyatPmax = 0
    PatMyMax = 0
    MyAtPmax = 0
    if b ≠ H or numTopBars ≠ numLayers
        needTwoCurves = true
        if isnan(barPattern)
          # The bar pattern is written as `m x n`, where m and n are integers.
          # So we reverse the order to look at moments about the y axis.
          barPattern = string(numTopBars, "f0") & "x" & string(numLayers, "f0")
        end
        Py, My, Myb, Pmax, MyMax, PatMyMax, MyAtPmax = curveArray(fc′, fy, H, b, Ag, numBars, numLayers, numTopBars, areaPerBar, barRadius, barPattern, cover, tieDia, margin)
    end

    # find scales and find the max M
    vertScale = 150 / Pmax
    
    Mb = 0
    Mmax = 0
    PatMmax = 0
    MatPmax = 0
    biggerMAxis = ""
    if MxMax > MyMax
        Mmax = MxMax
        PatMmax = PatMxMax
        Mb = Mxb
        MatPmax = MxAtPmax
        biggerMAxis = "x"
    else
        Mmax = MyMax
        PatMmax = PatMyMax
        Mb = Myb
        MatPmax = MyAtPmax
        biggerMAxis = "y"
    end
    horizScale = 200 / (1.05 * Mmax)

    fill "transparent"
    envelope = [baseX, baseY + vertScale * Pmax]
    for i in 2:length(Px)
        P = min(Pmax, Px[i])
        envelope = vcat(envelope, [baseX + horizScale * Mx[i], baseY + vertScale * P])
    end
    path envelope

    if needTwoCurves
        yEnvelope = [baseX, baseY + vertScale * Pmax]
        for i in 2:length(Py)
            P = min(Pmax, Py[i])
            yEnvelope = vcat(yEnvelope, [baseX + horizScale * My[i], baseY + vertScale * P])
        end
        path yEnvelope
    end

    # Write Pmax
    text [baseX, topOfGraph], string(Pmax, "f0") & " kips", "aboveright"
    # Mpure
    text [baseX + horizScale * Mb, baseY], string(Mb / 12, "f0"), "aboveleft"

    # PatMmax & Mmax
    xPos = baseX + horizScale * Mmax
    text [xPos, baseY + vertScale * PatMmax + 6], "ϕP~n~ = " & string(PatMmax, "f0") & " kips", "right"
    text [xPos, baseY + vertScale * PatMmax - 14], "ϕM~n~ = " & string(Mmax/12, "f0") & " kip·ft", "right"
    
    # M at Pmax
    text [baseX + (MatPmax * horizScale) - 8, topOfGraph + 8], "ϕM~n~ = " & string(MatPmax / 12, "f0") & " kip·ft", "right"

    # Draw the background lines
    iVert = backgroundLineSpacing(Pmax)
    iHoriz = backgroundLineSpacing(Mmax / 12)
    endX = baseX + 200
    horizScale = horizScale * 12
    strokewidth 0.5
    strokedasharray "5 10"

    for i in 1:20
        if (baseY + i * iVert * vertScale > topOfGraph)  break
        line [baseX, baseY + i * iVert * vertScale; endX, baseY + i * iVert * vertScale]
        text [baseX; baseY + i * iVert * vertScale], string(i * iVert, "f0"), "left"
    end

    for i in 1:20
        if (i * iHoriz * horizScale > 200) break
        line [baseX + i * iHoriz * horizScale, baseY; baseX + i * iHoriz * horizScale, topOfGraph]
        text [baseX + i * iHoriz * horizScale, baseY], string(i * iHoriz, "f0"), "below"
    end

    if Pu > 0 and Mu > 0
        fill "black"
        circle [baseX + Mu * horizScale, baseY + Pu * vertScale], 3
    end
end

function curveArray(fc′, fy, b, H, Ag, numBars, numTopBars, numLayers, areaPerBar, barRadius, barPattern, cover, tieDia, margin)
    # Calculate the interaction diagram for a concrete column.
    # Return the results in a pair of arrays.
    ϕb = 0.9    # ACI 318 strength reduction factor for bending
    ϕc = 0.65   # ditto for compression
    tieIR = 2 * tieDia  # inside radius of tie bending curve
    y = []          # y-coordinate of each layer of reinforcing steel
    A = []          # area of each layer of reinforcing steel
    P = []          # axial strength
    M = []          # bending strength
    β1 = {
        0.85 if fc′ <= 4;
        0.65 if fc′ >= 8;
        0.85 - (fc′ - 4) / 20 otherwise
    }
    As = numBars * areaPerBar
    ρ = As / Ag
  
    # find the locations of the layers of steel bars
    spacing = (H - 2 * margin) / (numLayers - 1)
    for i in 1:numLayers
        y = y & {margin if i == 1; y[i - 1] + spacing otherwise}
        numBarsInLayer = {
            numTopBars if i == 1 or i == numLayers;
            2 otherwise
        }
        A = A & (numBarsInLayer * areaPerBar)
    end

    # Find maximum compression, where moment = 0
    Po = 0.85 * fc′ * (Ag - As) + fy * As
    Pmax = 0.8 * ϕc * Po
  
    # Load the initial values of axial strength, P, and bending strength, M
    P = P & Pmax
    M = M & 0
  
    # get ready to start the while loop that will find axial and moment strengths
    phiHasChanged = false
    gotMatPmax = false
    ε_c = 0.003         # concreteStrain
    d = y[numLayers]  # depth from compression face to max steel layer
    # definition: ε_s = steelStrain = strain at "d"
    ε_s = 0.003         # initial value, starting at max compression
    
    # The upcoming loop will iterate with varying locations of the neutral axis.
    # That varies the steel strain, starting at max compression, and adding tension
    # to the steel in increments. At each iteration, we calculate column axial strength
    # and bending strenth. We'll iterate until we reach the pure bending condition
    # for the column.
    δ = -(fy / 29000 + ε_c) / 25   # difference in strain per step
    
    # Before we start the while loop, we're going to jump ahead to the point where a = h
    # a = h = β1*c = β1*d * (0.003 / (-ε_s + 0.003)), so from algebra:
    ε_s = 0.003 - 0.003 * β1 * d / H
    ε_s = ε_s - δ  # correct for the first increment
  
    while true
        ε_s = ε_s + δ
        c = d * (0.003 / (-ε_s + 0.003))  # distance to neutral axis
        s = (0.003 - ε_s) / d             # slope of strain line
        a = β1 * c                        # distance to edge of Whitney stress block
        Cc = 0                            # concrete compression chord force
        Mc = 0                            # moment about center of column due to Cc
  
        Cc = 0.85 * fc′ * b * min(a, H)      # Concrete compression chord force (first pass)
        Mc = {                               # Moment about center of column due to Cc
            Cc * (H / 2 - a / 2) if a < H;
            0 otherwise
        }
  
        Psteel = 0
        Msteel = 0
        for j in 1:numLayers
            # Find forces contributed by each layer of reinforcing bars.
            yLocal = c - y[j]
            localStrain = s * yLocal      # strain in this layer of bars
            fs = 29000 * localStrain      # steel stress
            fs = max(-fy, min(fs, fy))    # limited to yield stress
            PsLocal = A[j] * fs           # area of this layer times steel stress
            Psteel = Psteel + PsLocal
            localArm = H / 2 - y[j]
            Msteel = Msteel + PsLocal * localArm
            if a > y[j]
                CcLocal = A[j] * 0.85 * fc′
                Cc = Cc - CcLocal
                Mc = Mc - CcLocal * localArm
            end
        end
        
        P = P & ϕc * (Cc + Psteel)
        M = M & ϕc * (Mc + Msteel)
      
        if !gotMatPmax and P[length(P)] < Pmax
            # This is where the curve dips below Pmax
            MatPmax = M[length(M) - 1]
            gotMatPmax = true
        end
  
        if !phiHasChanged and ε_s <= -0.005
            iWhenPhiChanged = length(P)
            phiHasChanged = true
            δ = 2 * δ
        end
        if P[length(P)] < 0
            break
        end
    end
  
    # The last iteration of that loop just barely crossed into a point of negative axial load.
    # Do a linear interpolation to find M at the pure bending point (where axial load = 0).
    numPoints = length(P)
    n = numPoints - 1  # index of next-to-last point
    M = M[1:n-1] & M[n - 1] + (M[n] - M[n - 1]) / (P[n] - P[n - 1]) * (0 - P[n - 1])
    P = P[1:n-1] & 0
    numPoints = n
  
    # Rework the points after the change in phi
    numPointsAfterChange = numPoints - iWhenPhiChanged
    if numPointsAfterChange > 0
        f = (ϕb - ϕc) / ϕc
        δ = f / (numPoints - iWhenPhiChanged)
        f = [1:δ:(1 + f)]
        M = vcat(M[1:iWhenPhiChanged - 1], f ∘ M[iWhenPhiChanged:numPoints])
    end
    Mb = M[n]
  
    if ρ < 0.01
        # In low seismic areas, ACI 318 sections 10.6.1 and 10.3.1.2 allow ρ to less than 1%, but only
        # if the area of the column is assumed to be small enough to make an effective ρ of 1%.
        # Here, we do that reduction.
        r = ρ / 0.01  # reduction factor
        P = r × P
        M = r × M
        Pmax = r × Pmax
        Mb = r × Mb
    end
  
    Mmax = max(M)
    iMmax = findfirst(M == Mmax)
    PatMmax = P[iMmax]

    return P, M, Mb, Pmax, Mmax, PatMmax, MatPmax
end

function backgroundLineSpacing(num)
    exponent = floor(log10(num))
    numStr = string(num, "f0")
    n1 = floor(num / 10^exponent)  # first digit
    n2 = floor((num - n1 * 10^exponent) / 10^(exponent - 1))
    if n1 == 1
        arry = [0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.4, 0.4, 0.5, 0.5]
        n = arry[n2 + 1]
    else
        arry = [0, 0.5, 0.75, 1.0, 1.0, 1.5, 2.0, 2.0, 2.0]
        n = arry[n1]
    end
    return n * 10^exponent
end