Adding boundary condition to dsolve argument

8 次查看(过去 30 天)
I've used the dsolve function to get a solution for , but I want matlab to solve for such that . That is, I'm trying to add the boundary condition that when r = 0, α = . How can I input this boundary condition into the dsolve argument? I looked through the documentation but couldn't figure it out.
My code is below...
syms F_lift(r) C_lift(r) rho v_rel(r) A_blade r omega alpha(r) v_air P_in T_in
% r..........effective radius of turbine blade's rotation at a given distance along its length
% alpha......turbine blade's angle of attack at a given distance along its length
% F_lift.....lift force generated by turbine blade at a given distance along its length
% C_lift.....coefficient of lift for turbine airfoil at a given angle of attack
% v_rel......relative velocity between incoming air and turbine blade at a given distance along its length
% rho........mass density of the air after passing through the stator
% A_blade....surface area of turbine blade coming into contact with the airflow
% omega......angular velocity of the turbine and shaft
%%%%%%% major assumptions made
% Give the generic equation for lift of an airfoil
F_lift = 1/2*C_lift*rho*v_rel^2*A_blade
F_lift(r) = 
% Give C_lift as a function of alpha
C_lift = .12*alpha + .2;
% Give v_rel as a function of v_air and omega
omega = 10471.975511; % rad/s
v_rel = sqrt(v_air^2 + omega^2*r^2);
% Give area of each blade
A_blade = 5.8064e-5; % m^2
% Determine the air density at the turbine inlet using the ideal gas law
T_in = 1023.15; %%%%%%% 750 deg C
P_in = 1100*10^3; %%%%%%% ~10 atm = ~1100 kPa
R = 287.058; %%%%%%% R_air
rho = P_in/(R*T_in);
% Give v_air
v_rel = subs(v_rel,v_air,150); %%%%%%% assumed to be 150 m/s
% Give the new equation
F_lift = 1/2*C_lift*rho*v_rel^2*A_blade
F_lift(r) = 
% Let F_lift be constant along the length of the blade
liftDSTRBTN = diff(F_lift,r)
liftDSTRBTN(r) = 
alpha = dsolve(liftDSTRBTN == 0) % NEED HELP HERE!
alpha = 

采纳的回答

Cris LaPierre
Cris LaPierre 2021-10-14
I'm no expert here, so take this for what it is worth.
You mention needing to set a condition. Does the Solve Differential Equation with Conditions example help? You can then call dsolve with the conditions you have set.
If so, your code might look like this. One change I did make was to use alpha(r) everywhere to keep alpha a symfun. It may not matter, as I get the same result either way.
syms F_lift(r) C_lift(r) rho v_rel(r) A_blade r omega alpha(r) v_air P_in T_in
% Give the generic equation for lift of an airfoil
F_lift = 1/2*C_lift*rho*v_rel^2*A_blade
F_lift(r) = 
% Give C_lift as a function of alpha
C_lift = .12*alpha(r) + .2;
% ^^^^^^^^ -----------> used alpha(r) to preserver symfun
% Give v_rel as a function of v_air and omega
omega = 10471.975511; % rad/s
v_rel = sqrt(v_air^2 + omega^2*r^2);
% Give area of each blade
A_blade = 5.8064e-5; % m^2
% Determine the air density at the turbine inlet using the ideal gas law
T_in = 1023.15; %%%%%%% 750 deg C
P_in = 1100*10^3; %%%%%%% ~10 atm = ~1100 kPa
R = 287.058; %%%%%%% R_air
rho = P_in/(R*T_in);
% Give v_air
v_rel = subs(v_rel,v_air,150); %%%%%%% assumed to be 150 m/s
% Give the new equation
F_lift = 1/2*C_lift*rho*v_rel^2*A_blade
F_lift = 
% Let F_lift be constant along the length of the blade
liftDSTRBTN = diff(F_lift,r)
liftDSTRBTN = 
% Set condition
cond = alpha(0) == pi/2;
% vvv -----------> used alpha(r) to preserver symfun
alpha(r) = dsolve(liftDSTRBTN == 0,cond)
alpha(r) = 
You can test the solution at r=0 using subs. It does equal
subs(alpha,r,0)
ans(r) = 

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Calculus 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by