# Rectangle Vessels With Stay Plate - ASME VIII # Figure 13-2(a) - Sketch (7)

rectangle_vessel_one_stay_plate.R

# Rectangle Vessels With Stay Plates 
# Figure 13-2(a) - Sketch (7)

# 13-9 STAYED VESSELS OF RECTANGULAR CROSS SECTION 
# [FIGURE 13-2(A) SKETCHES (7) AND (8)]

#--- Begin single Script 

# Material SA-516 Gr 70 (S)
S = 20000
# Design Pressure @ 350F(P)
P = 900
# Inside Length (Short Side) (H)
H = 4.0 # U1-A Span ID
# Inside Length (Long Side) (h)
h = 4.0625
# Over all vessel length (Lv)
Lv = 126.25
# Plate Thickness (Short Side) (t1)
t1 = .750
# Plate Thickness (Long Side) (t2)
t2 = 1.0
# Plate Thickness (Stay Plate) (t3)
t3 = 5/8
# Corrosion Allowance (CA)
CA = .125
# Weld Joint Efficiency (Corner Joint) E = 1.0 100% UT and Spot RT
E = 1.0
# Tube Outside Diamter 
d = .750
# Tube Pitch
p = 1.625

# Ligament Efficincey 
em = (p - d) / p
em = .495
eb = .495

# (b) Vessel Stayed by a Single Plate. Figure 13-2(a)
# Figure 13-2(a) sketch (7) shows a vessel with a central stay plate.

# Adjust variables for corrosion.
h = h + 2*CA
# adjust h for stay plate 
#h = h / 2 + (t3/2) # note dim is to edge of stay plate
H = H +2*(CA)
t1 = t1 - CA
t2 = t2 - CA
t3 = t3 - CA

# The design equations in this paragraph
# are based on vessels in which the ratio of the length of the vessel to the long side or short side
# lengths (aspect ratio) is greater than four. These equations are conservatively applicable to vessels
# of aspect ratio less than four. Vessels with aspect ratios less than four may be designed in
# accordance with the provisions of U-2(g).

if (Lv/h > 4) {
  print("Aspect Ratio Satisfied")
} else if (Lv/h < 4) {
  print("Aspect Rati Not Satisifed")
}

# Paragraph 13-5 – Calculate the equation constants.
b = 1.0 
c_short_side = t1 / 2 # The sign of ci is positive. The sign of co is negative
c_long_side = t2 / 2 # The sign of ci is positive. The sign of co is negative
I1 = b*(t1^3) / 12
I2 = b*(t2^3) / 12
a = H/h
K = (I2/I1)*a

#-- Membrane Stress 
# Short Side Plates
Sm_eq1 = ((P * h) / (4 * t1)) * (4- ((2 + K*(5-a^2))/(1 + 2*K)))

# Long Side Plates
Sm_eq2 = (P * H) / (2* t2 * em)

# Stay Plate 
Sm_eq3 = ((P * h) / (2 * t3)) * ((2 + (K*(5 - (a^2)))) / (1 + 2 * K))

#-- Bending Stress 
# Short Side Plates
SbN_eq4 = ((P*c_short_side) / (24*I1)) * (-3*(H)^2 + 2*(h)^2 * ((1+2*(a^2)*K)/(1+2*K)))
SbN_eq4i = SbN_eq4 # Inside is negative
SbN_eq4o = -SbN_eq4 # Outside is positive
                               

SbQ_eq5 = ((P * h^2 * c_short_side) / (12 * I1)) * ((1 + 2*a^2*K) / (1 + 2 * K))
SbQ_eq5i = SbQ_eq5 # Inside is positive
SbQ_eq5o = -SbQ_eq5 # Outside is negative

# Long Side Plates
SbM_eq6 = ((P * h^2 * c_long_side) / (12 * I2 * eb)) * ((1 + K*(3-a^2))/(1 + 2*K))
SbM_eq6i = SbM_eq6 # Inside is positive
SbM_eq6o = -SbM_eq6 # Outside is negative

SbQ_eq7 = ((P * h^2 * c_long_side) / (12 * I2)) * ((1 + 2*a^2*K) / (1 + 2 * K))
SbQ_eq7i = SbQ_eq7 # Inside is positive
SbQ_eq7o = -SbQ_eq7 # Outside is negative

#-- Total Stress + Acceptance Criteria 
if (Sm_eq1 <= (S * E)) { 
  cat("True Short Side Membrane Stress",Sm_eq1,"psig is below Allowable Stress (S * E)", S*E,"psig")
} else if (Sm_eq1 > (S * E)) {
  cat("False - Short Side Membrane Stress",Sm_eq1,"psig is above Allowable Stress (S * E)", S*E,"psig")
}

# Short side plate at Location N, Membrane + Bending Stress:
StNi = Sm_eq1 + SbN_eq4i
StNo = Sm_eq1 + SbN_eq4o
if (StNi <= 1.5 * (S * E)) { 
  cat("True - Short Side Membrane + Bending Stress at location N (inside)",StNi,"psig is below Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
} else if (StNi > (S * E)) {
  cat("False - Short Side Membrane + Bending Stress at location N (inside)",StNi,"psig is above Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
}

if (StNo <= 1.5 * (S * E)) { 
  cat("True - Short Side Membrane + Bending Stress at location N (outside)",StNo,"psig is below Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
} else if (StNo > (S * E)) {
  cat("False - Short Side Membrane + Bending Stress at location N (outside)",StNo,"psig is above Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
}

# Short side plate at Location Q, Membrane + Bending Stress:
StQi = Sm_eq1 + SbQ_eq5i
StQo = Sm_eq1 + SbQ_eq5o
if (StQi <= 1.5 * (S * E)) { 
  cat("True - Short Side Membrane + Bending Stress at location Q (inside)",StQi,"psig is below Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
} else if (StQi > (S * E)) {
  cat("False - Short Side Membrane + Bending Stress at location Q (inside)",StQi,"psig is above Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
}

if (StQo <= 1.5 * (S * E)) { 
  cat("True - Short Side Membrane + Bending Stress at location Q (outside)",StQo,"psig is below Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
} else if (StQo > (S * E)) {
  cat("False - Short Side Membrane + Bending Stress at location Q (outside)",StQo,"psig is above Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
}

# Long side plate, Membrane Stress:
if (Sm_eq2 <= (S * E)) { 
  cat("True - Short Side Membrane Stress",Sm_eq2,"psig is below Allowable Stress (S * E)", S*E,"psig")
} else if (Sm_eq2 > (S * E)) {
  cat("False - Short Side Membrane Stress",Sm_eq2,"psig is above Allowable Stress (S * E)", S*E,"psig") 
}

# Long side plate at Location M, Membrane + Bending Stress:
StMi = Sm_eq2 + SbM_eq6i
StMo = Sm_eq2 + SbM_eq6o
if (StMi <= 1.5 * (S * E)) { 
  cat("True - Long Side Membrane + Bending Stress at location M (inside)",StMi,"psig is below Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
} else if (StMi > 1.5 * (S * E)) {
  cat("Flase - Long Side Membrane + Bending Stress at location M (inside)",StMi,"psig is above Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
}
if (StMo <= 1.5 * (S * E)) { 
  cat("True - Long Side Membrane + Bending Stress at location M (outside)",StMo,"psig is below Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
} else if (StMo > 1.5 * (S * E)) {
  cat("False - Long Side Membrane + Bending Stress at location M (outside)",StMo,"psig is above Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
}

# Long side plate at Location Q, Membrane + Bending Stress:
StQi = Sm_eq2 + SbQ_eq7i
StQo = Sm_eq2 + SbQ_eq7o
if (StQi <= 1.5 * (S * E)) { 
  cat("True - Long Side Membrane + Bending Stress at location Q (inside)",StQi,"psig is below Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
} else if (StQi > 1.5 * (S * E)) {
  cat("False - Long Side Membrane + Bending Stress at location Q (inside)",StQi,"psig is above Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
}
if (StQo <= 1.5 * (S * E)) { 
  cat("True - Long Side Membrane + Bending Stress at location Q (outside)",StQo,"psig is below Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
} else if (StQo > 1.5 * (S * E)) {
  cat("False - Long Side Membrane + Bending Stress at location Q (outside)",StQo,"psig is above Allowable Stress 1.5*(S * E)", 1.5*(S*E),"psig")
}

# Stay Plate Acceptance 
St = Sm_eq3
if (St <= (S * E)) { 
  cat("True Stay Plate Membrane Stress",St,"psig is below Allowable Stress (S * E)", S*E,"psig")
} else if (St > (S * E)) {
  cat("False - Stay Plate Membrane Stress",St,"psig is above Allowable Stress (S * E)", S*E,"psig")
}

#--- End single script 

#--- Begin iterative script 

# Iterate to find the failure thickness's 
# Material SA-516 Gr 70 (S)
S = 20000
# Design Pressure @ 350F(P)
P = 900
# Inside Length (Short Side) (H)
H = 4.0 # U1-A Span ID
# Inside Length (Long Side) (h)
h = 4.0625
# Over all vessel length (Lv)
Lv = 126.25
# Plate Thickness (Short Side) (t1)
t1 = .750
# Plate Thickness (Long Side) (t2)
t2 = 1.0
# Plate Thickness (Stay Plate) (t3)
t3 = 5/8
# Corrosion Allowance (CA)
CA = .125
# Weld Joint Efficiency (Corner Joint) E = 1.0 100% UT and Spot RT
E = 1.0
# Tube Outside Diamter 
do = .750
# Tube Pitch
p = 1.625

# Ligament Efficincey 
em = (p - d) / p
em = .495
eb = .495

# Adjust variables for corrosion.
h = h + 2*CA
# adjust h for stay plate 
H = H +2*(CA)

# Substitute t1 and t2 for a range of thickness
t1 = seq(from=0.01,to=1.750,by=0.001)
t2 = seq(from=0.01,to=2.0,by=0.001)
t3 = seq(from=0.01,to=1.625,by=0.001)

# Fix Lengths 
# find shortest vector length 
if (length(t1) < min(length(t2),length(t3))) {
  subset_vector_len = length(t1)
  subset_vector_len = length(t1)
  t2 = t2[(length(t2)-subset_vector_len):length(t2)]
  t3 = t3[(length(t3)-subset_vector_len):length(t3)]
  } else if (length(t2) < min(length(t1),length(t3))) {
    subset_vector_len = length(t2) 
    t1 = t1[(length(t1)-subset_vector_len):length(t1)]
    t3 = t3[(length(t3)-subset_vector_len):length(t3)]
} else if (length(t3) < min(length(t1),length(t2))) {
  subset_vector_len = length(t3)
  t1 = t1[(length(t1)-subset_vector_len):length(t1)]
  t2 = t2[(length(t2)-subset_vector_len):length(t2)]
}

# Adjust Variables For Corrosion
for(i in 1:length(t3)) {
t1[i] = t1[i] - CA
t2[i] = t2[i] - CA
t3[i] = t3[i] - CA
}

# Establish this appendix can be used
if (Lv/h > 4) {
  print("Satisfied")
} else if (Lv/h < 4) {
  print("Not Satisifed")
}

# Iterate through varying t1 and t2 thickness to find failure points 
# Initialize output 
output = data.frame(short_side_membrane_stress = numeric(0), short_side_membrane_plus_bending_stress_at_location_N_inside = numeric(0), short_side_membrane_plus_bending_stress_at_location_N_outside = numeric(0),
                    short_side_membrane_plus_bending_stress_at_location_Q_inside = numeric(0), short_side_membrane_plus_bending_stress_at_location_Q_outside = numeric(0),
                    long_side_membrane_stress = numeric(0),long_side_membrane_plus_bending_stress_at_location_M_inside = numeric(0),long_side_membrane_plus_bending_stress_at_location_M_outside = numeric(0),
                    long_side_membrane_plus_bending_stress_at_location_Q_inside = numeric(0), long_side_membrane_plus_bending_stress_at_location_Q_outside = numeric(0),stay_plate_membrane_stress = numeric(0), t1_thickness = numeric(0),t2_thickness = numeric(0), all_pass = numeric(0))

i=1
for (i in 1:length(t3)) { 
  # Paragraph 13-5 – Calculate the equation constants.
  b = 1.0 
  c_short_side = t1[i] / 2 # The sign of ci is positive. The sign of co is negative
  c_long_side = t2[i] / 2 # The sign of ci is positive. The sign of co is negative
  I1 = b*(t1[i]^3) / 12
  I2 = b*(t2[i]^3) / 12
  a = H/h
  K = (I2/I1)*a
  
  
  #-- Membrane Stress 
  # Short Side Plates
  Sm_eq1 = ((P * h) / (4 * t1[i])) * (4- ((2 + K*(5-a^2))/(1 + 2*K)))
  
  # Long Side Plates
  Sm_eq2 = (P * H) / (2* t2[i] * em)
  
  # Stay Plate 
  Sm_eq3 = ((P*h) / (2*t3[i])) * ((2 + (K*(5-a^2))) / (1 + 2*K))
  
  #-- Bending Stress 
  # Short Side Plates
  SbN_eq4 = ((P*c_short_side) / (24*I1)) * (-3*(H)^2 + 2*(h)^2 * ((1+2*(a^2)*K)/(1+2*K)))
  SbN_eq4i = SbN_eq4 # Inside is negative
  SbN_eq4o = -SbN_eq4 # Outside is positive
  
  
  SbQ_eq5 = ((P * h^2 * c_short_side) / (12 * I1)) * ((1 + 2*a^2*K) / (1 + 2 * K))
  SbQ_eq5i = SbQ_eq5 # Inside is positive
  SbQ_eq5o = -SbQ_eq5 # Outside is negative
  
  # Long Side Plates
  SbM_eq6 = ((P * h^2 * c_long_side) / (12 * I2 * eb)) * ((1 + K*(3-a^2))/(1 + 2*K))
  SbM_eq6i = SbM_eq6 # Inside is positive
  SbM_eq6o = -SbM_eq6 # Outside is negative
  
  SbQ_eq7 = ((P * h^2 * c_long_side) / (12 * I2)) * ((1 + 2*a^2*K) / (1 + 2 * K))
  SbQ_eq7i = SbQ_eq7 # Inside is positive
  SbQ_eq7o = -SbQ_eq7 # Outside is negative
  
  # Paragraphs 13-4(b), 13-4(c), and 13-7, Equations (7) through (10) – Acceptance Criteria:
  # 1 = true, 0 = false
  if (Sm_eq1 <= (S * E)) { 
    short_side_membrane_stress = 1 
  } else if (Sm_eq1 > (S * E)) {
    short_side_membrane_stress = 0
  }
  
  # Short side plate at Location N, Membrane + Bending Stress:
  StNi = Sm_eq1 + SbN_eq4i
  StNo = Sm_eq1 + SbN_eq4o
  if (StNi <= 1.5 * (S * E)) { 
    short_side_membrane_plus_bending_stress_at_location_N_inside = 1
  } else if (StNi > 1.5 * (S * E)) {
    short_side_membrane_plus_bending_stress_at_location_N_inside = 0
  }
  if (StNo <= 1.5 * (S * E)) { 
    short_side_membrane_plus_bending_stress_at_location_N_outside = 1
  } else if (StNo > 1.5 * (S * E)) {
    short_side_membrane_plus_bending_stress_at_location_N_outside = 0
  }
  
  # Short side plate at Location Q, Membrane + Bending Stress:
  StQi = Sm_eq1 + SbQ_eq5i
  StQo = Sm_eq1 + SbQ_eq5o
  if (StQi <= 1.5 * (S * E)) { 
    short_side_membrane_plus_bending_stress_at_location_Q_inside = 1
  } else if (StQi > 1.5 * (S * E)) {
    short_side_membrane_plus_bending_stress_at_location_Q_inside = 0
  }
  if (StQo <= 1.5 * (S * E)) { 
    short_side_membrane_plus_bending_stress_at_location_Q_outside = 1
  } else if (StQo > 1.5 * (S * E)) {
    short_side_membrane_plus_bending_stress_at_location_Q_outside = 0
  }

  # Long side plate, Membrane Stress:
  if (Sm_eq2 <= (S * E)) { 
    long_side_membrane_stress = 1 
  } else if (Sm_eq2 > (S * E)) {
    long_side_membrane_stress = 0
  }
  
  # Long side plate at Location M, Membrane + Bending Stress:
  StMi = Sm_eq2 + SbM_eq6i
  StMo = Sm_eq2 + SbM_eq6o
  if (StMi <= 1.5 * (S * E)) { 
    long_side_membrane_plus_bending_stress_at_location_M_inside = 1
  } else if (StMi > 1.5 * (S * E)) {
    long_side_membrane_plus_bending_stress_at_location_M_inside = 0
  }
  if (StMo <= 1.5 * (S * E)) { 
    long_side_membrane_plus_bending_stress_at_location_M_outside = 1
  } else if (StMo > 1.5 * (S * E)) {
    long_side_membrane_plus_bending_stress_at_location_M_outside = 0
  }
  
  # Long side plate at Location Q, Membrane + Bending Stress:
  StQi = Sm_eq2 + SbQ_eq7i
  StQo = Sm_eq2 + SbQ_eq7o
  if (StQi <= 1.5 * (S * E)) { 
    long_side_membrane_plus_bending_stress_at_location_Q_inside = 1
  } else if (StQi > 1.5 * (S * E)) {
    long_side_membrane_plus_bending_stress_at_location_Q_inside = 0
  }
  if (StQo <= 1.5 * (S * E)) { 
    long_side_membrane_plus_bending_stress_at_location_Q_outside = 1
  } else if (StQo > 1.5 * (S * E)) {
    long_side_membrane_plus_bending_stress_at_location_Q_outside = 0
  }
  
  # Stay Plate Membrane Stress
  if (Sm_eq3 <= (S * E)) { 
    stay_plate_membrane_stress = 1 
  } else if (Sm_eq3 > (S * E)) {
    stay_plate_membrane_stress = 0
  }
  
  t1_thickness = t1[i]
  t2_thickness = t2[i]
  
  # For varying t1,t2 check if Stresses at N,M and Q (inside and outside) Pass
  test_all_pass = sum(short_side_membrane_stress,short_side_membrane_plus_bending_stress_at_location_N_inside,short_side_membrane_plus_bending_stress_at_location_N_outside,
                      short_side_membrane_plus_bending_stress_at_location_Q_inside,short_side_membrane_plus_bending_stress_at_location_Q_outside,long_side_membrane_stress,
                      long_side_membrane_plus_bending_stress_at_location_M_inside,long_side_membrane_plus_bending_stress_at_location_M_outside,long_side_membrane_plus_bending_stress_at_location_Q_inside,
                      long_side_membrane_plus_bending_stress_at_location_Q_outside,stay_plate_membrane_stress)
  
  if(test_all_pass == 11) {
    all_pass = 1
  } else if (test_all_pass != 11) {
    all_pass = 0
  }
  
  
  temp_row = data.frame(short_side_membrane_stress,short_side_membrane_plus_bending_stress_at_location_N_inside,short_side_membrane_plus_bending_stress_at_location_N_outside,
                        short_side_membrane_plus_bending_stress_at_location_Q_inside,short_side_membrane_plus_bending_stress_at_location_Q_outside,long_side_membrane_stress,
                        long_side_membrane_plus_bending_stress_at_location_M_inside,long_side_membrane_plus_bending_stress_at_location_M_outside,long_side_membrane_plus_bending_stress_at_location_Q_inside,
                        long_side_membrane_plus_bending_stress_at_location_Q_outside,stay_plate_membrane_stress,stay_plate_membrane_stress,t1_thickness,t2_thickness,all_pass)
  output = rbind(output,temp_row)
  
  cat("This is Iteration",i,"\n")
  
} # end

# Find first passing thickness (tmin)
for (i in 2:nrow(output)) { 
  if(output$all_pass[i-1] == 0 & output$all_pass[i] == 1) {
    tmin_t1 = output$t1_thickness[i]
    tmin_t2 = output$t2_thickness[i]
  }
}

cat("Resulting Short Side Pressure Tmin = ",tmin_t1)
cat("Resulting Long Side Pressure Tmin = ",tmin_t2)
cat("Short Side Nominal Thickness = ",.750)
cat("Long Side Nominal Thickness = ",1.0)

# plot data
t1_out = output$t1_thickness
t2_out = output$t2_thickness
index = 1:nrow(output)

# Short side Pressure Tmin
#par(mfrow=c(1,1))
plot(x = index,y = t1_out, xlab = "No.", ylab = "t", main = "MC807A HAL-1030 A/B West Header Box\nThickness Needed For Stresses At Location N, M and Q To All Pass\nShort Side Pressure Tmin", cex.main=1.0,axes = FALSE,col = ifelse(t1_out < tmin_t1,'red','green'), pch = 19 )
axis(side = 1, at = 1:length(t1_out),labels = 1:length(t1_out), las =1,cex.axis=.7,tck = 0, lty = 2, col = "grey")
axis(side = 2, at = seq(from=0.01,to=3.05,by=0.05),labels = seq(from=0.01,to=3.05,by=0.05),  las =0,cex.axis=.7,tck = 0, lty = 2)
box()
abline(h=tmin_t1)
text(1000, .65, expression("Tmin = 0.573 / Tnominal = .750"), col = "black")
text(180, .5, expression("Thickness = Fail"), col = "red")
text(1380, 1.0, expression("Thickness = Pass"), col = "green")

# Long side Pressure Tmin
plot(x = index,y = t2_out, xlab = "No.", ylab = "t", main = "MC807A HAL-1030 A/B West Header Box\nThickness Needed For Stresses At Location N, M and Q To All Pass\nLong Side Pressure Tmin", cex.main=1.0,axes = FALSE,col = ifelse(t2_out < tmin_t2,'red','green'), pch = 19 )
axis(side = 1, at = 1:length(t2_out),labels = 1:length(t2_out), las =1,cex.axis=.7,tck = 0, lty = 2, col = "grey")
axis(side = 2, at = seq(from=0.01,to=3.05,by=0.05),labels = seq(from=0.01,to=3.05,by=0.05),  las =0,cex.axis=.7,tck = 0, lty = 2)
box()
abline(h=tmin_t2)
text(1000, .90, expression("Tmin = 0.823 / Tnominal = 1.0"), col = "black")
text(180, .72, expression("Thickness = Fail"), col = "red")
text(1450, 1.1, expression("Thickness = Pass"), col = "green")

#-- End iterative script 

#write.csv(output, file = "C:/Users/Andrew.Bannerman/Desktop/MARS/tmin_out_rec.csv")