Runge-Kutta integration of the Taylor-Maccoll eq

44 visualizzazioni (ultimi 30 giorni)
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!

Risposta accettata

Alan Stevens
Alan Stevens il 26 Lug 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 Commenti
Alan Stevens
Alan Stevens il 26 Lug 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
Robert Wake il 26 Lug 2021
Alan thank you so much! really appreciate it!

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Downloads in Help Center e File Exchange

Tag

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by