# Runge-Kutta integration of the Taylor-Maccoll eq

13 visualizaciones (últimos 30 días)
Robert Wake el 26 de Jul. de 2021
Comentada: Robert Wake el 26 de Jul. de 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!
##### 0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

Alan Stevens el 26 de Jul. de 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 comentariosMostrar 1 comentario más antiguoOcultar 1 comentario más antiguo
Alan Stevens el 26 de Jul. de 2021
LIke so (you had too many h's!).
% initial values for velocity just after the shock
vr = 0.9078;
vrdash = -0.1262;
theta_s = 14.35;
%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
%%%%%% DEFINE x BEFORE USING IT TO CONSTRUCT y
%initial y matrix
y1 = zeros(1, length(x));
y2 = zeros(1, length(x));
%initial conditions for velocity
y1(1) = vr;
y2(1) = vrdash;
g = 1.4; %ratio of specific heats
t(1) = 0; %%%%% t IS NEVER USED
%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)-1
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 NO NEED FOR ANY MORE h's
y1(i+1) = y1(i) + (1/6)*(k1+2*k2+2*k3+k4);
y2(i+1) = y2(i) + (1/6)*(m1+2*m2+2*m3+m4);
end
plot(x,y1),grid
xlabel('x'),ylabel('y1')
Robert Wake el 26 de Jul. de 2021
Alan thank you so much! really appreciate it!

Iniciar sesión para comentar.