Assigning values for a variable and solving symbolic equations numerically

6 views (last 30 days)
Rodrigo Pena
Rodrigo Pena on 13 Oct 2021
Commented: Rodrigo Pena on 14 Oct 2021
Hello, how are you?
I am trying to assign values for a variable 't' to solve a symbolic equation. I used the command solve to solve the equations algebrically to the desired variable 'd0' and I would like to assign variables for t and obtain 'd0' numerically. How can I do this ?
Please find my code below:
syms d0 t positive
h = 10; % Height of the tank [m]
H = 30; % Height of the base [m]
D = 10; % Diameter of the tank [m]
w = 700; % Wind pressure of the tank [N/m²]
e = 10/100; % Eccentricity [m]
delta_a = 20/100; % Allowable deflection [m]
gamma_w = 10*10^3; % Unit weight of water [N/m³]
gamma_s = 80*10^3; % Unit weight of steel [N/m³]
E = 210*10^9; % Modulus of elasticity [Pa]
t_t = 1.5/100; % Average thickness of the tank wall
sigma_b = 165 * 10^6; % Allowable bending stress [Pa]
R = d0/2 - t;
I = pi/64*(d0^4 - (d0 - 2*t)^4);
A = pi*t*(d0 - t);
r = sqrt(I/A); % radius of gyration [m]
sigma_a = (12*pi^2*E)/(92*(H/r)^2); % Allowable axial stress [Pa]
V = 1.2*pi*D^2 * h; % Volume of the tank [m³]
A_s = 1.25*pi*D^2; % Surface area of tank [m²]
A_p = (2*D*h)/3; % Projected area for wind loading [m²]
P = V*gamma_w + A_s*t_t*gamma_s; % Load on the column due to weight of water and steel tank
W = w*A_p; % Lateral load at the tank C.G [m]
delta1 = (W*H^2)/(12*E*I) * (4*H + 3*h);
delta2 = H/(2*E*I) * (0.5*W*h + P*e)*(H + h);
delta = delta1 + delta2; % Deflection at C.G of the tank
M = W*(H + 0.5*h) + (delta + e)*P; % Moment at base [N.m]
f_b = M/(2*I) *d0; % Bending stress
f_a = V*gamma_w + A_s*t_t*gamma_s; % Axial stress
t = [1/100:0.01:20/100]; % t range #ok<NBRAK>
g1 = sigma_a - f_a == 0;
g1_d0 = solve(g1,d0)
g2 = sigma_b - f_b == 0;
g2_d0 = solve(g2,d0)
g3 = delta_a - delta ==0;
g3_d0 = solve(g3,d0)
g4 = R - 20 == 0;
g4_d0 = solve(g4,d0)
g5 = -R + 0.35 == 0;
g5_d0 = solve(g5,d0)
Thanks in advance !

Accepted Answer

Paul
Paul on 13 Oct 2021
You can use subs() and double(), at least for the "easy" versions of d0:
syms d0 t positive
h = 10; % Height of the tank [m]
H = 30; % Height of the base [m]
D = 10; % Diameter of the tank [m]
w = 700; % Wind pressure of the tank [N/m²]
e = 10/100; % Eccentricity [m]
delta_a = 20/100; % Allowable deflection [m]
gamma_w = 10*10^3; % Unit weight of water [N/m³]
gamma_s = 80*10^3; % Unit weight of steel [N/m³]
E = 210*10^9; % Modulus of elasticity [Pa]
t_t = 1.5/100; % Average thickness of the tank wall
sigma_b = 165 * 10^6; % Allowable bending stress [Pa]
R = d0/2 - t;
I = pi/64*(d0^4 - (d0 - 2*t)^4);
A = pi*t*(d0 - t);
r = sqrt(I/A); % radius of gyration [m]
sigma_a = (12*pi^2*E)/(92*(H/r)^2); % Allowable axial stress [Pa]
V = 1.2*pi*D^2 * h; % Volume of the tank [m³]
A_s = 1.25*pi*D^2; % Surface area of tank [m²]
A_p = (2*D*h)/3; % Projected area for wind loading [m²]
P = V*gamma_w + A_s*t_t*gamma_s; % Load on the column due to weight of water and steel tank
W = w*A_p; % Lateral load at the tank C.G [m]
delta1 = (W*H^2)/(12*E*I) * (4*H + 3*h);
delta2 = H/(2*E*I) * (0.5*W*h + P*e)*(H + h);
delta = delta1 + delta2; % Deflection at C.G of the tank
M = W*(H + 0.5*h) + (delta + e)*P; % Moment at base [N.m]
f_b = M/(2*I) *d0; % Bending stress
f_a = V*gamma_w + A_s*t_t*gamma_s; % Axial stress
tvals = [1/100:0.01:20/100]; % t range #ok<NBRAK> % so as not to conflict with the sym variable t
g1 = sigma_a - f_a == 0;
g1_d0 = solve(g1,d0,'ReturnConditions',true)
g1_d0 = struct with fields:
d0: [2×1 sym] parameters: [1×0 sym] conditions: [2×1 sym]
g2 = sigma_b - f_b == 0;
g2_d0 = solve(g2,d0,'ReturnConditions',true)
g2_d0 = struct with fields:
d0: z1 parameters: z1 conditions: 260082878480646144000000000000000*pi^2*t^8 + 1820580149364523008000000000000000*pi^2*t^6*z1^2 + 1105352233542746112000000000000000*pi^2*t^4*z1^4 + 65020719620161536000000000000000*pi^2*t^2*z1^6 + 17182393990160644319477760000000*pi*t^4…
g3 = delta_a - delta ==0;
g3_d0 = solve(g3,d0,'ReturnConditions',true)
g3_d0 = struct with fields:
d0: x parameters: x conditions: - 112742891520000000*pi*t^4 + 225485783040000000*pi*t^3*x - 169114337280000000*pi*t^2*x^2 + 56371445760000000*pi*t*x^3 - 31731444346091383 == 0 & 0 < x
g4 = R - 20 == 0;
g4_d0 = solve(g4,d0,'ReturnConditions',true)
g4_d0 = struct with fields:
d0: 2*t + 40 parameters: [1×0 sym] conditions: symtrue
g5 = -R + 0.35 == 0;
g5_d0 = solve(g5,d0,'ReturnConditions',true)
g5_d0 = struct with fields:
d0: 2*t + 7/10 parameters: [1×0 sym] conditions: symtrue
double(subs(g5_d0.d0,t,sym(tvals)))
ans = 1×20
0.7200 0.7400 0.7600 0.7800 0.8000 0.8200 0.8400 0.8600 0.8800 0.9000 0.9200 0.9400 0.9600 0.9800 1.0000 1.0200 1.0400 1.0600 1.0800 1.1000
You'll have to decide what you want to do with g1_d0, which has two solutions, and g2_d0 and g3_d0 which have to satisfy conditions.
  3 Comments
Rodrigo Pena
Rodrigo Pena on 14 Oct 2021
aah ok. I managed to implement what you wrote above. I will search more on how to get rid of these z1 and x variables that appears in the calculations. thanks !

Sign in to comment.

More Answers (0)

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!