Assigning values for a variable and solving symbolic equations numerically
12 次查看(过去 30 天)
显示 更早的评论
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 !
0 个评论
采纳的回答
Paul
2021-10-13
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)
g2 = sigma_b - f_b == 0;
g2_d0 = solve(g2,d0,'ReturnConditions',true)
g3 = delta_a - delta ==0;
g3_d0 = solve(g3,d0,'ReturnConditions',true)
g4 = R - 20 == 0;
g4_d0 = solve(g4,d0,'ReturnConditions',true)
g5 = -R + 0.35 == 0;
g5_d0 = solve(g5,d0,'ReturnConditions',true)
double(subs(g5_d0.d0,t,sym(tvals)))
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 个评论
Paul
2021-10-14
编辑:Paul
2021-10-14
You can access the second solution of g1_d0 by
g1_d0.d0(2)
Some equations only have solutions when certain conditions are true. Consdier the equation
(x+y)/(x+z) = 0.
This equation has a trivial solution, but only as long as y is not the same as z, which is what solve() show
syms x y z
sol = solve((x+y)/(x+z)==0,'ReturnConditions',true)
Other equations have more than one (or infinite) solutions that can still be expressed in terms of free parameter. Consider the equation
sin(x) = 0
Any value of x that is an integer multiple of pi is solution, which is what solve() shows
sol = solve(sin(x)==0,'ReturnConditions',true)
Note here that sol.conditions shows the restriction that the parameter, k, must be an integer.
I didn't look too closely at the equations to see what's going on with g2_d0 and g3_d0.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!