MATLAB Answers

Runge-Kutta integration of the Taylor-Maccoll eq

조회 수: 6(최근 30일)
Robert Wake
Robert Wake 2021년 7월 26일
댓글: Robert Wake 2021년 7월 26일
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
Alan Stevens 2021년 7월 26일
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
Robert Wake
Robert Wake 2021년 7월 26일
Alan thank you so much! really appreciate it!

댓글을 달려면 로그인하십시오.

추가 답변(0개)

태그

제품

Community Treasure Hunt

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

Start Hunting!

Translated by