Problems in solving Taylor-Maccoll Equation via ode45

Hello everyone,

For a study project, I have to solve the Taylor-Maccoll Equation that describes the hypersonic filed around a cone.

My textbook provides the system of ODEs in this fashion

The γ and the are constants. Now I have tried to implement the ode45 via writing a function that codes the second members of the system. In particular, I have substituted

I don't know if this make sense, but then I have solved for x and y in this way

This is how I have implemented

function f = TaylorMaccoll(omega, V)

global gamma V_lim

Vr     = V(1);
Vomega = V(2);


% a building

tris          = ((V_lim*V_lim) - (Vomega*Vomega) - (Vr*Vr));
den           = (gamma - 1)*(tris);

rap  = (2*Vomega*Vomega)/(den);

a = 1 - rap;

% b building

b = (2*Vomega*Vr)/den;

% Fs building

F1 = Vomega;
F2 = -((Vomega*cotd(omega))+(2*Vr));
Fq = (b/a)*F1 + (1/a)*F2;

f = [F1; Fq];
end

Then, in the main, I have called the ode45 in this way

V0  = [V_inf*cosd(beta_c), epsilon*V_inf*sind(beta_c)];

[omega,V] = ode45(@TaylorMaccol,[beta_c:-0.001:eps],V0);

The stepsize is negative because I have to integrate from the shock angle to the body wall. The problem is that the solution is not stable

I have seen on the internet several other attempts that actually work, but I would like to follow my class and my textbook. I am pretty sure that the error is in the way I implement the ODEs but I have never seen this kind of problem so far. If some one can provide me some help I would be very grateful, also with some hints I don't want the entire code.

[Alan Stevens]

You will need to upload the whole code, or specify all the data values, so we can test possible soutions.

[Leonardo Molino]

Sure, I will provide you the main

clc; clear all; close all;

global V_lim gamma

R      = 287;
gamma  = 1.4;
cp     = (gamma*R)/(gamma-1);

Ma      = 10;
T_inf   = 237;
p_inf   = 0.1;
h_inf   = cp*T_inf;
a0      = sqrt(gamma*R*T_inf);
V_inf   = Ma*a0;
H       = (V_inf*V_inf)/2 + h_inf;
T0      = H/cp;
p0_inf  = (T0/T_inf)^(gamma/(gamma-1))*p_inf;
rho_inf = 0.129;

n      = 2/(gamma - 1);
V_lim   = sqrt(2*H);


delta_c = 10;
%beta_c  = deg2rad(ObliqueShockSolver(n, Ma, delta_c, 0.001)); % Can use 14.42 deg 
beta_c = deg2rad(14.42);

primo     = 1/(n+1);
rapporto1 = n/(Ma*Ma*sin(beta_c)*sin(beta_c));
secondo   = 1 + rapporto1;

epsilon = primo*secondo;

V0  = [V_inf*cos(beta_c), epsilon*V_inf*sin(beta_c)];

[sol] = ode45(@TaylorMaccoll, [beta_c:-0.01:eps], V0);

omega = convang(sol.x,'rad','deg');
plot(omega,sol.y(2,:))

function f = TaylorMaccoll(omega, V)
global gamma V_lim
Vr     = V(1);
Vomega = V(2);
% a building
tris          = ((V_lim*V_lim) - (Vomega*Vomega) - (Vr*Vr));
den           = (gamma - 1)*(tris);
rap  = (2*Vomega*Vomega)/(den);
a = 1 - rap;
% b building
b = (2*Vomega*Vr)/den;
% Fs building
F1 = Vomega;
F2 = -((Vomega*cotd(omega))+(2*Vr));
Fq = (b/a)*F1 + (1/a)*F2;
f = [F1; Fq];
end

The problem is that I have to stop the integration when the . I have read something about the 'Events' option, but I don't really know how to implement it.

[Torsten]

There is still code missing (see above).

[Leonardo Molino]

Oh yes I know, I have added the value you can replace :D Is 14.42 deg

[Torsten]

Why do you use cotd and not cot since omega seems to be in radians ?

And why cot(omega) and not cot(g*omega) ?

[Leonardo Molino]

Regarding the first question, it is beacuse some times I forget the "d" so I prefer to convert everything to rad and using the classic "cos" or "sin".

Regarding the second one, it is a typo of my book, in italian literature we write cotg as cotangent.

[Alan Stevens]

In F1 you have cotd(omega), though you pass omega as radians to the function.

Also, your initial value of Vomega is positive. Should it be?

[Torsten]

@Leonardo Molino

Regarding the first question, it is beacuse some times I forget the "d" so I prefer to convert everything to rad and using the classic "cos" or "sin".

But you use cotd instead of cot although omega is in radians ...

[Leonardo Molino]

@Torsten You are right, I have made a little bit of confusion... I will fix this.

@Alan Stevens I came up with this idea. If, as you said, we start to integrate from a negative value of the Vomega and If we keep the "positive" step size, probably the system may converge. Also, I don't need to integrate up to the zero value of Vomega, I can stop integrating when the Vomega changes its sign, because if this happens, this mean that I have reached the "wall". The problem is, how can I instruct the ode45 to stop when the Vomega turns positive?

[Torsten]

clc; clear all; close all;

global V_lim gamma

R      = 287;
gamma  = 1.4;
cp     = (gamma*R)/(gamma-1);

Ma      = 10;
T_inf   = 237;
p_inf   = 0.1;
h_inf   = cp*T_inf;
a0      = sqrt(gamma*R*T_inf);
V_inf   = Ma*a0;
H       = (V_inf*V_inf)/2 + h_inf;
T0      = H/cp;
p0_inf  = (T0/T_inf)^(gamma/(gamma-1))*p_inf;
rho_inf = 0.129;

n      = 2/(gamma - 1);
V_lim   = sqrt(2*H);


delta_c = 10;
%beta_c  = deg2rad(ObliqueShockSolver(n, Ma, delta_c, 0.001)); % Can use 14.42 deg 
beta_c = deg2rad(14.42);

primo     = 1/(n+1);
rapporto1 = n/(Ma*Ma*sin(beta_c)*sin(beta_c));
secondo   = 1 + rapporto1;

epsilon = primo*secondo;

V0  = [V_inf*cos(beta_c), epsilon*V_inf*sin(beta_c)];

options = odeset('Events',@myEvent);
[sol] = ode45(@TaylorMaccoll, [beta_c:-0.01:eps], V0, options);

omega = convang(sol.x,'rad','deg');
plot(omega,sol.y(2,:))

function f = TaylorMaccoll(omega, V)
  global gamma V_lim
  Vr     = V(1);
  Vomega = V(2);
  % a building
  tris          = ((V_lim*V_lim) - (Vomega*Vomega) - (Vr*Vr));
  den           = (gamma - 1)*(tris);
  rap  = (2*Vomega*Vomega)/(den);
  a = 1 - rap;
  % b building
  b = (2*Vomega*Vr)/den;
  % Fs building
  F1 = Vomega;
  F2 = -((Vomega*cotd(omega))+(2*Vr));
  Fq = (b/a)*F1 + (1/a)*F2;
  f = [F1; Fq];
end

function [value,isterminal,direction] = myEvent(t,y)
  value = y(2);
  isterminal = 1;
  direction = 0;
end

[Leonardo Molino]

@Torsten Thank you very much, my idea seems to work! That straight line in the beginning of the plot is very promising

[Torsten]

The event function makes the solver stop when y(2) = Vomega = 0.

[Leonardo Molino]

I have decided to completely rewrite the ODEs.

The idea I came up with is basically to wrote everthing in terms of , so in this way I can easly explicit the . Starting from this, I simply have wrote the system of ODEs via substituting and

function dy = TyMc (omega, y)

global Vlim gamma

dy = zeros(2,1);

A = 2/(gamma - 1);

a = 1 - ((A)*(y(2)^2/(Vlim^2-y(2)^2-y(1)^2)));

b = (A)*((y(2)*y(1))/(Vlim^2-y(2)^2-y(1)^2));

dy(1) = y(2);
dy(2) = (b/a)*y(2) - (1/a)*((y(2)*cot(omega))+(2*y(1)));

end

Thanks to your help, I know how to stop the iteration. The BC can be negative (the ) with a positive step size or positive with a negative one (I have chosen the first one).

With huge surprise, the whole thing works and I have solved the hypersonic field around the cone!

Thank you very much for your help!