MATLAB Answers

Runge-Kutta integration of the Taylor-Maccoll eq

1 view (last 30 days)
Robert Wake
Robert Wake on 26 Jul 2021
Commented: Robert Wake on 26 Jul 2021
Hi, I'm trying to integrate the Taylor-Maccoll equation in order to determine the velocity flowfield in a conical shockwave, moving at hypersonic speed. As explained in the photo ive attached, the equation can be represented as a standard ODE with solution vector. I have initial conditions for both vr and vr', which I am then supposed to place into the second vector, y'. I then integrate y' using the Runge-Kutta method and at some point in the integration vr' should equal (roughly) 0. I have a step size of h, which defines the increment in theta (polar coordinates) I take from the shockwave, towards the cone angle.
I am aware of two codes which have been uploaded to the MATLAB servers which I can download, however these codes do not include the discrete values of velocity at every angle of theta, and so I must code these a different way.
In my code I have manually typed up the Runge-Kutta method but I think there is some problem in my code as my vrdash values are not equating to zero at the correct theta value (noweher near in fact!).
% initial values for velocity just after the shock
vr = 0.9078;
vrdash = -0.1262;
theta_s = 14.35;
%initial y matrix
y1 = zeros(1, length(x));
y2 = zeros(1, length(x));
%initial conditions for velocity
y1(1) = vr;
y2(1) = vrdash;
%defining the interval and step size (thetas_s = half angle of shock wave)
h = -0.01; %negative step size as moving towards zero in theta direction
x = theta_s:h:0; %integrate between shock angle and zero
g = 1.4; %ratio of specific heats
t(1) = 0;
%defining function from y' equation in photo
A = (g-1)/2;
F1 = @(t, y1, y2) y2;
F2 = @(t, y1, y2) (y1*y2^2-A*(1-y1^2-y2^2)*(2*y1+y2*cotd(theta_s)))...
/(A*(1-y1^2-y2^2) - y2^2);
%initial loop for Runge-Kutta method
for i=1:length(x)
k1 = h*F1(t, y1(i), y2(i));
m1 = h*F2(t, y1(i), y2(i));
k2 = h*F1(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*k1);
m2 = h*F2(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*m1);
k3 = h*F1(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*k2));
m3 = h*F2(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*m2));
k4 = h*F1(t+h, (y1(i)+h),(y2(i)+k3*h));
m4 = h*F2(t+h, (y1(i)+h),(y2(i)+m3*h));
% main equation
y1(i+1) = y1(i) + (1/6)*(k1+2*k2+2*k3+k4)*h;
y2(i+1) = y2(i) + (1/6)*(m1+2*m2+2*m3+m4)*h;
end
end
I'm sorry if I havent explained this properly, if there's anything i'm missing I'd be happy to explain. Any help would be greatly appreciated as I've been stuck on this code for nearly a month now. Thanks!

Accepted Answer

Alan Stevens
Alan Stevens on 26 Jul 2021
Shouldn't
k1 = h*F1(t, y1(i), y2(i));
m1 = h*F1(t, y1(i), y2(i));
k2 = h*F2(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*k1);
m2 = h*F2(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*k1);
k3 = h*F1(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*k2));
m3 = h*F2(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*k2));
k4 = h*F1(t+h, (y1(i)+h),(y2(i)+k3*h));
m4 = h*F2(t+h, (y1(i)+h),(y2(i)+k3*h));
be
k1 = h*F1(t, y1(i), y2(i));
m1 = h*F2(t, y1(i), y2(i));
k2 = h*F1(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*k1);
m2 = h*F2(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*m1);
k3 = h*F1(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*k2));
m3 = h*F2(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*m2));
k4 = h*F1(t+h, (y1(i)+h),(y2(i)+k3*h));
m4 = h*F2(t+h, (y1(i)+h),(y2(i)+m3*h));
  3 Comments
Robert Wake
Robert Wake on 26 Jul 2021
Alan thank you so much! really appreciate it!

Sign in to comment.

More Answers (0)

Tags

Products

Community Treasure Hunt

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

Start Hunting!

Translated by