Last active
May 1, 2020 06:17
-
-
Save supersonic71/968b25f3e3e1c48b4721917b5450280c to your computer and use it in GitHub Desktop.
code
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
%NAMED CONSTANTS | |
RHO = 1; | |
T = 2; | |
V = 3; | |
var = 3; | |
n_points = 31; | |
n_time = 14000; | |
sol = zeros(var, n_points, n_time); | |
dx = .1; | |
%initiliaze at t = 0 | |
% A = 1 + 2.2(x-1.5)^2 | |
C = .02; | |
X = (0:n_points -1)./10; | |
A = 1+2.2.*(X-1.5).^2; | |
for i = 1:n_points | |
sol(RHO,i,1) = 1 - .3146.*X(i); | |
sol(T,i,1) = 1 - .2314.*X(i); | |
sol(V,i,1) = (.1 + 1.09.*X(i)).*(sol(T,i,1).^.5); | |
end | |
dts = zeros(1,n_points); | |
drho1 = zeros(1,n_points); | |
dV1 = zeros(1,n_points); | |
dT1 = zeros(1,n_points); | |
drho2 = zeros(1,n_points); | |
dV2 = zeros(1,n_points); | |
dT2 = zeros(1,n_points); | |
gamma = 1.4; | |
for t = 1:n_time-1 | |
%predictor step | |
for i = 2:n_points-1 | |
drho1(i) = -sol(RHO,i,t).*(sol(V,i+1,t)-sol(V,i,t))./dx - sol(RHO,i,t).*sol(V,i,t).*(log(A(i+1)) - log(A(i)))./dx - sol(V,i,t).*(sol(RHO,i+1,t)-sol(RHO,i,t))./dx; | |
dV1(i) = -sol(V,i,t).*(sol(V,i+1,t)-sol(V,i,t))./dx - (1./gamma).*( (sol(T,i+1,t) - sol(T,i,t))./dx + (sol(T,i,t)./sol(RHO,i,t)).*(sol(RHO,i+1,t)-sol(RHO,i,t))./dx ); | |
dT1(i) = -sol(V,i,t).*(sol(T,i+1,t)-sol(T,i,t))./dx - (gamma -1).*(sol(T,i,t)).*( (sol(V,i+1,t) - sol(V,i,t))./dx +sol(V,i,t).*(log(A(i+1)) - log(A(i)))./dx ); | |
dts(1,i) = C.*(dx./(sol(T,i,t).^.5 + sol(V,i,t))); | |
end | |
dt = min(dts(1,2:end-1)); | |
for i = 2:n_points-1 | |
sol(RHO,i,t+1) = sol(RHO,i,t) + drho1(i).*dt; | |
sol(V,i,t+1) = sol(V,i,t) + dV1(i).*dt; | |
sol(T,i,t+1) = sol(T,i,t) + dT1(i).*dt; | |
end | |
sol(V,1,t+1) = sol(V,1,t); | |
sol(RHO,1,t+1) = sol(RHO,1,t); | |
sol(T,1,t+1) = sol(T,1,t); | |
sol(V,n_points,t+1) = sol(V,n_points,t); | |
sol(RHO,n_points,t+1) = sol(RHO,n_points,t); | |
sol(T,n_points,t+1) = sol(T,n_points,t); | |
for i = 2:n_points-1 | |
drho2(i) = -sol(RHO,i,t+1).*(sol(V,i,t+1)-sol(V,i-1,t+1))./dx - sol(RHO,i,t+1).*sol(V,i,t+1).*(log(A(i)) - log(A(i-1)))./dx - sol(V,i,t+1).*(sol(RHO,i,t+1)-sol(RHO,i-1,t+1))./dx; | |
dV2(i) = -sol(V,i,t+1).*(sol(V,i,t+1)-sol(V,i-1,t+1))./dx - (1./gamma).*( (sol(T,i,t+1) - sol(T,i-1,t+1))./dx + (sol(T,i,t+1)./sol(RHO,i,t+1)).*(sol(RHO,i,t+1)-sol(RHO,i-1,t+1))./dx ); | |
dT2(i) = -sol(V,i,t+1).*(sol(T,i,t+1)-sol(T,i-1,t+1))./dx - (gamma -1).*(sol(T,i,t+1)).*( (sol(V,i,t+1) - sol(V,i-1,t+1))./dx +sol(V,i,t+1).*(log(A(i)) - log(A(i-1)))./dx ); | |
drho = mean([drho1(i) drho2(i)]); | |
dV = mean([dV1(i) dV2(i)]); | |
dT = mean([dT1(i) dT2(i)]); | |
sol(RHO,i,t+1) = sol(RHO,i,t) + drho.*dt; | |
sol(V,i,t+1) = sol(V,i,t) + dV.*dt; | |
sol(T,i,t+1) = sol(T,i,t) + dT.*dt; | |
end | |
sol(V,1,t+1) = 2.*sol(V,2,t+1) - sol(V,3,t+1); | |
sol(RHO,1,t+1) = sol(RHO,1,t); | |
sol(T,1,t+1) = sol(T,1,t); | |
sol(V,n_points,t+1) = 2.*sol(V,n_points-1,t+1) - sol(V,n_points-2,t+1); | |
sol(RHO,n_points,t+1) = 2.*sol(RHO,n_points-1,t+1) - sol(RHO,n_points-2,t+1); | |
sol(T,n_points,t+1) = 2.*sol(T,n_points-1,t+1) - sol(T,n_points-2,t+1); | |
fprintf(num2str(t) + "\n"); | |
end |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
change the time step(line 46) to dt = min( c*dx ./ (sqrt( sol(T, : ,t) ) + sol(V,:,t) ) )
replace the line 46 to line 39