Assigning values for a variable and solving symbolic equations numerically

24 ビュー (過去 30 日間)
Rodrigo Pena
Rodrigo Pena 2021 年 10 月 13 日
コメント済み: Rodrigo Pena 2021 年 10 月 14 日
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 !

採用された回答

Paul
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)
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 件のコメント
Paul
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)
sol = struct with fields:
x: -y parameters: [1×0 sym] conditions: y ~= z
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)
sol = struct with fields:
x: pi*k parameters: k conditions: in(k, 'integer')
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.
Rodrigo Pena
Rodrigo Pena 2021 年 10 月 14 日
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 !

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeConversion Between Symbolic and Numeric についてさらに検索

製品


リリース

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by