19/07/2024, 02:44 TensionStiffening

Concrete Assignment 4 Question 2b

In [1]: import numpy as np
import scipy
from scipy.optimize import root
from scipy.optimize import brentq
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
sns.set(rc={"figure.dpi":100, "savefig.dpi":300})
sns.set_context("notebook")
sns.set_style("ticks")

In [2]: num = 1000

e_c = -0.0035
e_cc = np.linspace(-0.00002,e_c,num)
B = (4-e_cc/e_c)/(6-2*(e_cc-e_c))
a = 1/B*(e_cc/e_c-(e_cc/e_c)**2/3)

fig, ax = plt.subplots()
plt.plot(e_cc,B,label='B')
plt.xlabel('$\epsilon_{cc}$')
plt.title('Rectangular Stress Block Constants vs $\epsilon_{cc}$')
plt.plot(e_cc,a,label='a')
plt.legend()
plt.grid()

file:///Users/benjaminklassen/Downloads/TensionStiffening.html 1/8
19/07/2024, 02:44 TensionStiffening

In [3]: # Initial Matrix Definitions

ep = np.zeros(num)
es_top = np.zeros(num)
es_bot = np.zeros(num)
ec_bot = np.zeros(num)
fp = np.zeros(num)
fs_bot = np.zeros(num)
fs_top = np.zeros(num)
fc_bot = np.zeros(num)
c = np.zeros(num)
c2 = np.zeros(num)
Cc = np.zeros(num)
Cs = np.zeros(num)
C = np.zeros(num)
Tp = np.zeros(num)
Ts = np.zeros(num)
Tc = np.zeros(num)
T = np.zeros(num)

In [4]: #Constant Givens

## Section Properties
d_bot = 720
d_top = 50
b = 400
h = 800
Ag = b*h

file:///Users/benjaminklassen/Downloads/TensionStiffening.html 2/8
19/07/2024, 02:44 TensionStiffening

yg = h/2
Ig = 1/12*b*h**3

## Prestressed Tendons
A_RO = 0.015
B_RO = 119
C_RO = 18
Ep = 200000 #MPa
fpu = 1860 #MPa
fpy = 0.9*fpu
db = 16 #mm

## Concrete
f_c = -40
Ec = 4500*np.sqrt(-f_c)
f_t=0.6*np.sqrt(-f_c)
e_t = f_t/Ec # Approximate for strain at rupture
print(f'Ec: {Ec:1.2f} MPa')
print(f'f\'t: {f_t:1.2f} MPa')
print(f'e_t: {e_t:1.3e}')

## Transformation

n_p = Ep/Ec

Ec: 28460.50 MPa

f't: 3.79 MPa
e_t: 1.333e-04

i) Case One
In [5]: # Variables that change
Ap = 1000 #mm^2
fpo = 1350 #MPa
del_ep = fpo/Ep
e_cc = np.linspace(-0.00043,-0.0035,num)
B = (4-e_cc/e_c)/(6-2*(e_cc-e_c))
a = 1/B*(e_cc/e_c-(e_cc/e_c)**2/3)

In [6]: ## Transformations
Ac = Ag - Ap
Ap_tr = n_p*Ap

A_tr = Ac + Ap_tr
y_tr = (Ag*yg +(n_p-1)*(Ap*d_bot))/A_tr
I_tr = (Ig + Ag*(yg-y_tr)**2 + (n_p-1)*Ap*(d_bot-y_tr)**2)

print(f'Transformed Area: {A_tr:1.2f}mm^2')

print(f'Transformed y_bar (from top): {y_tr:1.2f}mm')
print(f'Transformed I: {I_tr:1.2f}mm^2')

Transformed Area: 326027.28mm^2

Transformed y_bar (from top): 405.92mm
Transformed I: 17672450421.00mm^2

file:///Users/benjaminklassen/Downloads/TensionStiffening.html 3/8
19/07/2024, 02:44 TensionStiffening

In [7]: ## Decompression Forces

N = 0
N_0 = Ap * fpo # No offsets to subtract
M_0 = -N_0 * (y_tr - d_bot) / 1000**2

## To find Mcr
e_cen = (N - N_0) / (Ec * A_tr)
f_cen = e_cen * Ec
M_cr = (f_t - f_cen) * I_tr / (h - y_tr) / 1000**2 + M_0
Phi_zero_mom = (0 - M_0 * 1000**2) / (Ec * I_tr) # rad/mm
Phi_M_cr = (M_cr - M_0) * 1000**2 / (Ec * I_tr)

## Moment-curvature
phi_1 = np.linspace(Phi_zero_mom, Phi_M_cr, 10)
M_1 = phi_1 * Ec * I_tr / 1000**2 + M_0

print(f'Cracking Moment: {M_cr:1.2f} kNm')

print(f'Decompression Moment: {M_0:1.2f} kNm')
print(f'Phi at Cracking Moment: {Phi_M_cr: 1.4e}')
print(f'Phi at Zero Moment: {Phi_zero_mom: 1.4e}')

Cracking Moment: 779.88 kNm

Decompression Moment: 424.01 kNm
Phi at Cracking Moment: 7.0753e-07
Phi at Zero Moment: -8.4303e-07

In [8]: for i in range(num):

e_conc = e_cc[i]
alpha = a[i]
beta = B[i]

def f(c):

# Strains
e_p = e_conc*(c-d_bot)/c+del_ep

e_c_bot = e_conc*(c-d_bot)/c

#Stresses
f_p = Ep*e_p*(A_RO+(1-A_RO)/(1+(B_RO*e_p)**C_RO)**(1/C_RO))
f_c_bot=e_c_bot*Ec

# Stress Limits
if f_c_bot>f_t:
Ac_eff = ((h-d_bot)+7.5*db)*b
pp = Ap/Ac_eff
alpha_1 = 0.7 #Tendons Only
alpha_2 = 1 # Short-term load assumed
fc1 = alpha_1*alpha_2*f_t/(1+np.sqrt(500*e_c_bot))
fc2 = pp*(fpy-f_p)
f_c_bot = max(min(fc1,fc2),0)
if f_p > fpy:
f_p = fpy

# Forces
Cc = (alpha*f_c)*(beta*c)*b

file:///Users/benjaminklassen/Downloads/TensionStiffening.html 4/8
19/07/2024, 02:44 TensionStiffening

C = Cc

A_ct = (h-c)*b #Assumed total area in tension since db not provided

Tp = Ap*f_p
Tc = A_ct*f_c_bot

T=Tp+Tc

return C+T

c[i] = brentq(f,1,800)

# Strains
ep[i] = e_conc*(c[i]-d_bot)/c[i]+del_ep
ec_bot[i] = e_conc*(c[i]-d_bot)/c[i]

if Ap==0:
ep[i]=0

#Stresses
fp[i] = Ep*ep[i]*(A_RO+(1-A_RO)/(1+(B_RO*ep[i])**C_RO)**(1/C_RO))
fc_bot[i]=ec_bot[i]*Ec
if fp[i] > fpy: #fpy used with tension stiffening instead of fpu
fp[i] = fpy

# Stress Limits
if fc_bot[i]>f_t:
Ac_eff = ((h-d_bot)+7.5*db)*b
pp = Ap/Ac_eff
alpha_1 = 0.7 #Tendons Only
alpha_2 = 1 # Short-term load assumed
fc1 = alpha_1*alpha_2*f_t/(1+np.sqrt(500*ec_bot[i]))
fc2 = pp*(fpy-fp[i])
fc_bot[i] = max(min(fc1,fc2),0)

# Forces
Cc[i] = (alpha*f_c)*(beta*c[i])*b
C[i] = (Cc[i])/1000

A_ct = (h-c[i])*b #Assumed total area in tension since db not provided

Tp[i] = Ap*fp[i]
Tc[i] = A_ct*fc_bot[i]

T[i]=(Tp[i]+Tc[i])/1000

# Tabulated Values
col_names = ['e_cc', 'c (mm)','a (mm)', 'B (mm)', 'ep','ec_bot', 'fc_bot','f
'T (kN)', 'Tp', 'Tc', 'C (kN)']
table = pd.DataFrame(list(zip(e_cc, c, a, B, ep, ec_bot,fc_bot, fp,
T,Tp/1000,Tc/1000,C)),columns=col_names)
print(table)

file:///Users/benjaminklassen/Downloads/TensionStiffening.html 5/8
19/07/2024, 02:44 TensionStiffening

e_cc c (mm) a (mm) B (mm) ep ec_bot fc_bot

\
0 -0.000430 716.590695 0.182153 0.646852 0.006752 0.000002 0.058224
1 -0.000433 713.103604 0.183440 0.646705 0.006754 0.000004 0.119199
2 -0.000436 709.736152 0.184727 0.646558 0.006756 0.000006 0.179510
3 -0.000439 706.481317 0.186015 0.646411 0.006758 0.000008 0.239198
4 -0.000442 703.332684 0.187301 0.646264 0.006760 0.000010 0.298302
.. ... ... ... ... ... ... ...
995 -0.003488 157.214545 1.329422 0.500587 0.019235 0.012485 0.000000
996 -0.003491 157.145010 1.330400 0.500441 0.019253 0.012503 0.000000
997 -0.003494 157.075658 1.331379 0.500294 0.019271 0.012521 0.000000
998 -0.003497 157.006488 1.332356 0.500147 0.019289 0.012539 0.000000
999 -0.003500 156.937500 1.333333 0.500000 0.019307 0.012557 0.000000

fp (MPa) T (kN) Tp Tc C (kN)

0 1348.983941 1350.926527 1348.983941 1.942586 -1350.926527
1 1349.403894 1353.547086 1349.403894 4.143192 -1353.547086
2 1349.819225 1356.300518 1349.819225 6.481293 -1356.300518
3 1350.230225 1359.178018 1350.230225 8.947793 -1359.178018
4 1350.637159 1362.171596 1350.637159 11.534437 -1362.171596
.. ... ... ... ... ...
995 1674.000000 1674.000000 1674.000000 0.000000 -1674.000000
996 1674.000000 1674.000000 1674.000000 0.000000 -1674.000000
997 1674.000000 1674.000000 1674.000000 0.000000 -1674.000000
998 1674.000000 1674.000000 1674.000000 0.000000 -1674.000000
999 1674.000000 1674.000000 1674.000000 0.000000 -1674.000000

[1000 rows x 12 columns]

In [9]: M = ((Tp)*d_bot + Tc*(h-c)/2 + Cc*B*c/2)/(1000**2)

phi = -e_cc/c
index = np.where(phi>Phi_M_cr)[0][0]
phi = phi[index:]
M = M[index:]

phi_1 = np.hstack((phi_1,phi))
M_1 = np.hstack((M_1,M))
Mn = M_1.max()
print(f'Nominal Moment Capacity: {Mn:1.2f} kNm')
fig, ax = plt.subplots()
plt.plot(phi_1,M_1,label='i)')
plt.xlabel('Curvature ($\phi$)')
plt.ylabel('Moment (kNm)')
plt.title('Moment vs Curvature')
plt.grid()

Nominal Moment Capacity: 1139.60 kNm

file:///Users/benjaminklassen/Downloads/TensionStiffening.html 6/8
19/07/2024, 02:44 TensionStiffening

In [10]: # To find phi(M=1000 kNm)

differences = np.abs(M_1 - 1000)
closest_index = np.argmin(differences)
phi_index = phi_1[closest_index]
print(phi_index)

3.693359432366014e-06

In [11]: Ac_eff = ((h-d_bot)+7.5*db)*b

pp = Ap/Ac_eff
alpha_1 = 0.7 #Tendons Only
alpha_2 = 1 # Short-term load assumed
fc1 = alpha_1*alpha_2*f_t/(1+np.sqrt(500*ec_bot))
fc2 = pp*(fpy-fp)

In [12]: plt.plot(ec_bot,fc1,label='fc1')
plt.plot(ec_bot,fc2,label='fc2')
plt.legend()

Out[12]: <matplotlib.legend.Legend at 0x143dd0f50>

file:///Users/benjaminklassen/Downloads/TensionStiffening.html 7/8
19/07/2024, 02:44 TensionStiffening

In [ ]:

file:///Users/benjaminklassen/Downloads/TensionStiffening.html 8/8