20 views (last 30 days)

Show older comments

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!

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));

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

Start Hunting!