Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Replies: 1   Last Post: Nov 23, 2012 4:20 PM

 Messages: [ Previous | Next ]
 ARTEMIS CIVIL Posts: 3 Registered: 11/23/12
Posted: Nov 23, 2012 4:01 PM

Hello, I really have to deal with the problem below. I am a beginner in matlab and my deadline is quite pressing. Thank you in advance!
Well, I have to guess an initial value "n" in order to find the ultimate value "Pu". When the extracting Pu deviates from a specific value (Pe) more than 0.1% then I need the code to make use of the last value of "n" until the parameter "Pu"reach the Pe. Otherwise the parameter n is kept.
herein is the code:
function [Pu] = KN (c,B,d,rho,fck,fsy,Es,eta,fcube,h)
% Starting Assumption
c = 569;
r1 = 130;
r2 = 90;
B = (130 * 90)^0.5;
Ap = 0.1397;
h = 43;
rho = Ap / h;
d = 0.5 * h;
fck = 48.6;
sigmap = 1.84;
Fp = sigmap * h;
fpk = 2576;
fsy = fpk - (Fp / Ap);
Es = 1.83 * 10^5;
fcube = 0.8 * fck;
eta=0.45;
ddelta = 1;
j = 0;
delta = h;
while abs(ddelta/delta) > 0.01
j = j + 1;
if j > 300
break
end
Fc = 0.8 * 2/3 * fck * (h/2 - delta/4);
Ft = d * rho * fsy;
Fbmax = Fc - Ft
Mbmax = (Ft * (2 * d - h) - Fc * (d - 13*h/16 - 3 * delta/32));
Fb = eta * Fbmax
Mb = eta * Mbmax
X = calX(c,B,d,Es,fsy,fcube,rho,Fb,Mb);
y = caly(c,B,d,Es,fsy,fcube,rho,Fb,X);
if B / d < 2
tasi = 0.0035 * (1 - 0.22*(B/d))*(1 + B/(2*y));
else
tasi = 0.0019 * ( 1 + B/(2*y));
end
delta1 = 1/2*tasi*(c-B);
ddelta = delta - delta1;
delta = delta1;
end
TA = calTA (c,B,d,y,X);
p1 = P1 (B,y,d,fcube,TA);
p2 = P2 (c,B,y,d,Es,fsy,X,rho,Fb);
function [P1] = P1 (B,y,d,fcube,TA)
% B = diameter of the load area, given in [mm]
% y = depth of the compression zone in KN model, given in [mm]
% d = average effective depth of the slab, given in [mm]
% fcube = compressive strength of concrete, measured on standard cubes,
% given in [N/mm2]
% TA = tan(alpha)
%
% calculate ft
%
if B/d < 2
ft = 825 * (0.35 + 0.3* (fcube/150))*(1 - 0.22 * (B / d));
else
ft = 460 * (0.35 + 0.3 * (fcube/150));
end
%
% calculate falpha
%
falpha = TA * (1 - TA)/(1 + TA * TA);
%
% calculate P1
P1 = pi * (B / d) * (y / d) * (B + 2 * y) / (B + y) * ft * falpha * d * d /1000;
function [P2] = P2 (c,B,y,d,Es,fsy,X,rho,Fb)
% c = diameter of the slab, given in [mm]
% B = diameter of the load area, given in [mm]
% y = depth of the compression zone in KN model, given in [mm]
% d = average effective depth of the slab, given in [mm]
% fsy = yield stress of reinforcing steel, given in [N/mm2]
% Es = modulus of elastisity of reinforcing steel, given in [N/mm2]
% X = boundary moment and normal force ratio
% rho = reinforcement percentage
% Fb = horizontal lateral normal force, given in [kN]
%
% Calculate kz
ky = 3*(c - B)/(2*(3*d - y));
kz = ky - 3 * X * c / (4 * (3 * d - y));
% Calculate R1 R2overBeta
if B / d < 2
tasi = 0.0035 * (1 - 0.22*(B/d))*(1 + B/(2*y));
else
tasi = 0.0019 * ( 1 + B/(2*y));
end
rs = Es / fsy * tasi * (d - y);
C0 = B/2 + 1.8 * d;
if rs > C0
R1 = rho * fsy * d * ((rs - C0) + rs * log(c/(2 * rs)))/1000;
R2overBeta = rho * fsy * d * C0/1000;
else
R1 = rho * fsy * d * rs * log(c/(2*C0))/1000;
R2overBeta = rho * fsy * d * rs/1000;
end
%
% Calculate P2
%
P2 = 2 * pi / kz * (R1 + R2overBeta + Fb * (c/2/1000));
function [ta] = calTA (c,B,d,y,X)
ky = 3*(c - B)/(2*(3*d - y));
kz = ky - 3 * X * c / (4 * (3 * d - y));
A = 1/4.7*(1+y/B)*log(c/(B+2*y));
ta = ((kz + 1)- sqrt((kz + 1)*(kz + 1) - 4 * (kz + A)*(A + 1)))/(2*(kz+A));
function [y] = caly(c,B,d,Es,fsy,fcube,rho,Fb,X)
dp = 100;
j = 0;
y = d;
while abs(dp)>0.1
j = j+1;
if j > 300
break
end
TA = calTA (c,B,d,y,X);
p1 = P1(B,y,d,fcube,TA);
p2 = P2(c,B,y,d,Es,fsy,X,rho,Fb);
dp = p1 - p2;
y = y - dp/1000;
end
function [X] = calX(c,B,d,Es,fsy,fcube,rho,Fb,Mb)
dX = 1;
k = 0;
X = 0;
X1 = 1;
while abs (dX) > 0.0000001
k = k + 1;
if k > 300;
break
end
y = caly(c,B,d,Es,fsy,fcube,rho,Fb,X);
TA = calTA (c,B,d,y,X);
p1 = P1 (B,y,d,fcube,TA);
p2 = P2 (c,B,y,d,Es,fsy,X,rho,Fb);
X1 = 4*pi*(Mb/((p1+p2)*1000/2));
dX = X1 - X;
X = X + dX/20;
end

the code is a bit long but tahnk you again for your time and the advices!!!

Date Subject Author
11/23/12 ARTEMIS CIVIL
11/23/12 dpb