Adding boundary condition to dsolve argument
12 views (last 30 days)
Show older comments
Caleb Rudick
on 13 Oct 2021
Commented: Caleb Rudick
on 15 Oct 2021
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
% 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
% Let F_lift be constant along the length of the blade
liftDSTRBTN = diff(F_lift,r)
alpha = dsolve(liftDSTRBTN == 0) % NEED HELP HERE!
0 Comments
Accepted Answer
Cris LaPierre
on 14 Oct 2021
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
% 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
% Let F_lift be constant along the length of the blade
liftDSTRBTN = diff(F_lift,r)
% Set condition
cond = alpha(0) == pi/2;
% vvv -----------> used alpha(r) to preserver symfun
alpha(r) = dsolve(liftDSTRBTN == 0,cond)
You can test the solution at r=0 using subs. It does equal
subs(alpha,r,0)
More Answers (0)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!