Skip to content

Instantly share code, notes, and snippets.

/TubeProblemLW.m Secret
Created Aug 23, 2017

Embed
What would you like to do?
function TubeProblemLW
%Resolution by explicit LaxWendroff. Boundary conditions imposed by
%implicit discretization of compatibility equations.
%Not fully optimized.
close all
clc
%Points at the mesh and length of the tube:
N = 100;
L = 1;
xv = linspace(0,1,N);
Deltax = L/(N-1);
%Final t and relation Deltat/Deltatmax:
rDeltat = 0.95;
tfin = 0.05;
%Artificial viscosity:
D = 0;
%Solver options:
options = optimset('TolFun',1e-8,'TolX',1e-8,'display','off');
%Gas properties:
Rg = 287;
gamma = 1.4;
cp = gamma/(gamma-1)*Rg;
cv = 1/(gamma-1)*Rg;
%Boundary conditions:
pa = 1e5;
Ta = 300;
% pr = @(t) pa+0.5*pa*atan(t/0.001)/(pi/2);
pr = @(t) 1.2*pa;
Tr = 300;
%Initial conditions:
pv = pa*ones(1,N);
Tv = Ta*ones(1,N);
rhov = pv./Rg./Tv;
vv = zeros(1,N);
av = sqrt(gamma*pv./rhov);
%Initialize matrices containing the values of p,rho and v at each time:
pm = pv;
rhom = rhov;
vm = vv;
%Loop until tfin
tn = 0;
[Wn,Fn] = CompWF(reshape([pv;rhov;vv],[],1),gamma);
Wn = reshape(Wn,3,[]);
Fn = reshape(Fn,3,[]);
iter = 0;
while tn<=tfin
iter = iter+1;
%Compute Deltat:
Deltat = rDeltat*Deltax/max(abs(vv)+av);
tn1 = tn+Deltat;
%Add artificial viscosity to Fn:
Fn(:,2:end) = Fn(:,2:end) - D*Deltax*(Wn(:,2:end)-Wn(:,1:end-1));
%Prediction of W_(j+0.5)^(n+0.5) and F_(j+0.5)^(n+0.5):
Wpred = 0.5*(Wn(:,1:end-1)+Wn(:,2:end)) - ...
Deltat/2/Deltax*(Fn(:,2:end)-Fn(:,1:end-1));
Fpred = compF(Wpred,gamma);
%Add artificial viscosity to Fpred:
Fpred(:,1:end-1) = Fpred(:,1:end-1) - D*Deltax*(Wpred(:,2:end)-Wpred(:,1:end-1));
%Next W in the center points:
Wn1 = Wn(:,2:end-1) - Deltat/Deltax * (Fpred(:,2:end)-Fpred(:,1:end-1));
%Get p,rho,T at j=2 and t=t^(n+1):
[p2,rho2,v2] = compU(Wn1(:,1),gamma);
%Get third Riemann invariant at j=2 ant t=t^(n+1):
r3_2 = v2 - 2/(gamma-1)*sqrt(gamma*p2/rho2);
%Get p,rho,T at j=1 and t=t^(n):
[p1,rho1,v1] = compU(Wn(:,1),gamma);
%Get third Riemann invariant at j=1 ant t=t^(n):
r3_1_0 = v1 - 2/(gamma-1)*sqrt(gamma*p1/rho1);
%Get p,rho,u at the j=1 and t=t_n^1 imposing boundary conditions:
% u1sol = Broyden(@(u)fLeft(u,r3_1_0,r3_2,pr(tn1),cp*Tr,Deltat,Deltax,gamma),...
% [p2;rho2;v2],[],options,[0,0,-Inf],Inf(1,3));
% u1sol = NewtonRaphson(@(u)fLeft(u,r3_1_0,r3_2,pr(tn1),cp*Tr,Deltat,Deltax,gamma),...
% [p2;rho2;v2],options,[0,0,-Inf],Inf(1,3));
u1sol = fsolve(@(u)fLeft(u,r3_1_0,r3_2,pr(tn1),cp*Tr,Deltat,Deltax,gamma),...
[p2;rho2;v2],options);
if u1sol(3)<0
u1sol(3) = v2;
end
%Get p,rho,T at j=2 and t=t^(n+1):
[pN_1,rhoN_1,vN_1] = compU(Wn1(:,end),gamma);
%Get first and second Riemann invariants at j=N-1 ant t=t^(n+1):
r1_Nm1 = pN_1/rhoN_1^gamma;
r2_Nm1 = vN_1 + 2/(gamma-1)*sqrt(gamma*pN_1/rhoN_1);
%Get p,rho,T at j=N and t=t^(n):
[pN,rhoN,vN] = compU(Wn(:,end),gamma);
%Get first and second Riemann invariants at j=N ant t=t^(n):
r1_0 = pN/rhoN^gamma;
r2_0 = vN + 2/(gamma-1)*sqrt(gamma*pN/rhoN);
%Get p,rho,u at j=N and t=t_n^1 imposing boundary conditions:
% uNsol = Broyden(@(u)fRight(u,r1_Nm1,r1_0,r2_Nm1,r2_0,pa,Deltat,Deltax,gamma),...
% [pN_1;rhoN_1;vN_1],[],options,0,Inf);
% uNsol = NewtonRaphson(@(u)fRight(u,r1_Nm1,r1_0,r2_Nm1,r2_0,pa,Deltat,Deltax,gamma),...
% [pN_1;rhoN_1;vN_1],options,0,Inf);
uNsol = fsolve(@(u)fRight(u,r1_Nm1,r1_0,r2_Nm1,r2_0,pa,Deltat,Deltax,gamma),...
[pN_1;rhoN_1;vN_1],options);
if uNsol(3)<0
uNsol(3) = vN_1;
end
%Complete with W at boundaries:
Wn1 = cat(2,CompWF(u1sol,gamma),Wn1,CompWF(uNsol,gamma));
Fn1 = compF(Wn1,gamma);
%Extract new values of p, T and v:
[pv,rhov,vv] = compU(Wn1,gamma);
av = sqrt(gamma*pv/rhov);
pm = cat(2,pm,pv);
rhom = cat(2,rhom,rhov);
vm = cat(2,vm,vv);
%Actualize variables:
Wn = Wn1;
Fn = Fn1;
tn = tn1;
%plot variables:
figure(1)
subplot(2,2,1)
plot(xv,pv/pa,'-*')
ylabel('p')
% title(['p_1 = ',num2str(pr(tn1)/pa)])
%
subplot(2,2,2)
plot(xv,rhov,'-*')
ylabel('\rho')
%
subplot(2,2,3)
plot(xv,vv,'-*')
ylabel('v')
%
subplot(2,2,4)
plot(xv,abs(vv)./av)
ylabel('M')
title(['t = ',num2str(tn1)])
drawnow
% disp([p1/pa,p2/pa,abs(vv(end)/av(end))])
if ~isreal(vv) | isnan(vv)
keyboard
end
end
end
function F = compF(W,gamma)
%Receives matrix W and computes F:
%Extract p,rho,v from W:
rho = W(1,:);
v = W(2,:)./rho;
p = (gamma-1)*(W(3,:)-0.5*rho.*v.^2);
%Compute F:
F = [rho.*v; p+rho.*v.^2; gamma/(gamma-1)*p.*v+0.5*rho.*v.^3];
end
function [p,rho,v] = compU(W,gamma)
%Receives matrix W and computes U:
%Extract p,rho,v from W:
rho = W(1,:);
v = W(2,:)./rho;
p = (gamma-1)*(W(3,:)-0.5*rho.*v.^2);
end
function f = fLeft(u,r3_1_0,r3_2,pLeft,h0Left,Deltat,Deltax,gamma)
%Receives u=[p1;rho1;v1] and r3_(j=2) and imposes p1=pLeft, h01=hLeft and
%D_(u-a) w_3 = 0:
%Extract p,rho,v:
p1 = u(1);
rho1 = u(2);
v1 = u(3);
c1 = sqrt(gamma*p1/rho1);
%Compute r3:
r3 = v1 - 2/(gamma-1)*c1;
%Residual:
f = [ p1-pLeft;
gamma/(gamma-1)*p1/rho1 + v1^2/2 - h0Left;
(r3 - r3_1_0)/Deltat + (v1-c1)*(r3_2-r3)/Deltax ];
end
function f = fRight(u,r1_Nm1,r1_N,r2_Nm1,r2_N,pRight,Deltat,Deltax,gamma)
%Receives u=[pN;rhoN;vN] and imposes pN=pRight, D_u(w1)=0 and
%D_(u+a)(w2) = 0:
%Extract p,rho,v:
pN = u(1);
rhoN = u(2);
vN = u(3);
cN = sqrt(gamma*pN/rhoN);
%Compute r1 and r2:
r1 = pN/rhoN^gamma;
r2 = vN + 2/(gamma-1)*cN;
%Residual:
f = [ pN-pRight;
(r1 - r1_N)/Deltat + vN*(r1 - r1_Nm1)/Deltax;
(r2 - r2_N)/Deltat + (vN+cN)*(r2 - r2_Nm1)/Deltax ];
end
function TubeProblemLWexplicit
%Resolution by explicit LaxWendroff. Boundary conditions imposed by
%implicit discretization of compatibility equations.
%Not fully optimized.
close all
clc
%Points at the mesh and length of the tube:
N = 100;
L = 1;
xv = linspace(0,1,N);
Deltax = L/(N-1);
%Final t and relation Deltat/Deltatmax:
rDeltat = 0.95;
tfin = 0.05;
%Artificial viscosity:
D = 8;
%Solver options:
options = optimset('TolFun',1e-8,'TolX',1e-8,'display','off');
%Gas properties:
Rg = 287;
gamma = 1.4;
cp = gamma/(gamma-1)*Rg;
cv = 1/(gamma-1)*Rg;
%Boundary conditions:
pa = 1e5;
Ta = 300;
% pr = @(t) pa+0.1*pa*atan(t/0.001)/(pi/2);
pr = @(t) 1.2*pa;
Tr = 300;
%Initial conditions:
pv = pa*ones(1,N);
Tv = Ta*ones(1,N);
rhov = pv./Rg./Tv;
vv = zeros(1,N);
av = sqrt(gamma*pv./rhov);
%Initialize matrices containing the values of p,rho and v at each time:
pm = pv;
rhom = rhov;
vm = vv;
%Loop until tfin
tn = 0;
[Wn,Fn] = CompWF(reshape([pv;rhov;vv],[],1),gamma);
Wn = reshape(Wn,3,[]);
Fn = reshape(Fn,3,[]);
iter = 0;
while tn<=tfin
iter = iter+1;
%Compute Deltat:
Deltat = rDeltat*Deltax/max(abs(vv)+av);
tn1 = tn+Deltat;
%Add artificial viscosity to Fn:
Fn(:,2:end) = Fn(:,2:end) - D*Deltax*(Wn(:,2:end)-Wn(:,1:end-1));
%Prediction of W_(j+0.5)^(n+0.5) and F_(j+0.5)^(n+0.5):
Wpred = 0.5*(Wn(:,1:end-1)+Wn(:,2:end)) - ...
Deltat/2/Deltax*(Fn(:,2:end)-Fn(:,1:end-1));
Fpred = compF(Wpred,gamma);
%Add artificial viscosity to Fpred:
Fpred(:,1:end-1) = Fpred(:,1:end-1) - D*Deltax*(Wpred(:,2:end)-Wpred(:,1:end-1));
%Next W in the center points:
Wn1 = Wn(:,2:end-1) - Deltat/Deltax * (Fpred(:,2:end)-Fpred(:,1:end-1));
%Get p,rho,T at j=1,2 and t=t^(n):
[p12,rho12,v12] = compU(Wn(:,[1 2]),gamma);
c12 = sqrt(gamma*p12./rho12);
%Get third Riemann invariant at j=1,2 and t=t^(n):
[r3_12] = v12 - 2/(gamma-1)*sqrt(gamma*p12./rho12);
%Compute third Riemann invariant at j=1 and t=t^(n+1):
r3_1 = r3_12(1)-Deltat*(v12(1)-c12(1))*(r3_12(2)-r3_12(1))/Deltax;
%Get p,rho,u at the j=1 and t=t_n^1 imposing boundary conditions:
% u1sol = Broyden(@(u)fLeft(u,pr(tn1),cp*Tr,r3_1,gamma),...
% [p23(1);rho23(1);v23(1)],[],options,[0,0,-Inf],Inf(1,3));
% u1sol = NewtonRaphson(@(u)fLeft(u,pr(tn1),cp*Tr,r3_1,gamma),...
% [p12(2);rho12(2);v12(2)],options,[0,0,-Inf],Inf(1,3));
u1sol = fsolve(@(u)fLeft(u,pr(tn1),cp*Tr,r3_1,gamma),...
[p12(2);rho12(2);v12(2)],options);
if u1sol(3)<0
u1sol(3) = v12(2);
end
%Get p,rho,T at j=N-1,N and t=t^(n):
[pN_10,rhoN_10,vN_10] = compU(Wn(:,end-1:end),gamma);
%Get first and second Riemann invariants at j=N-2,N-1 ant t=t^(n+1):
r1_Nm10 = pN_10./rhoN_10.^gamma;
r2_Nm10 = vN_10 + 2/(gamma-1)*sqrt(gamma*pN_10./rhoN_10);
%Compute Riemann invariants at j=N and t=t^(n+1):
r1_N = r1_Nm10(end) - Deltat*vN_10(end)*(r1_Nm10(end)-r1_Nm10(end-1))/Deltax;
r2_N = r2_Nm10(end) - Deltat*(vN_10(end)+vN_10(end))*(r2_Nm10(end)-r2_Nm10(end-1))/Deltax;
%Get p,rho,u at j=N and t=t_n^1 imposing boundary conditions:
% uNsol = Broyden(@(u)fRight(u,pa,r1_N,r2_N,gamma),...
% [pN_21(2);rhoN_21(2);vN_21(2)],[],options,[0,0,-Inf],Inf(1,3));
% uNsol = NewtonRaphson(@(u)fRight(u,pa,r1_N,r2_N,gamma),...
% [pN_10(1);rhoN_10(1);vN_10(1)],options,[0,0,-Inf],Inf(1,3));
uNsol = fsolve(@(u)fRight(u,pa,r1_N,r2_N,gamma),...
[pN_10(1);rhoN_10(1);vN_10(1)],options);
if uNsol(3)<0
uNsol(3) = vN_10(end-1);
end
%Complete with W at boundaries:
Wn1 = cat(2,CompWF(u1sol,gamma),Wn1,CompWF(uNsol,gamma));
Fn1 = compF(Wn1,gamma);
%Extract new values of p, T and v:
[pv,rhov,vv] = compU(Wn1,gamma);
av = sqrt(gamma*pv/rhov);
pm = cat(2,pm,pv);
rhom = cat(2,rhom,rhov);
vm = cat(2,vm,vv);
%Actualize variables:
Wn = Wn1;
Fn = Fn1;
tn = tn1;
%plot variables:
figure(1)
subplot(2,2,1)
plot(xv,pv/pa,'-*')
ylabel('p')
% title(['p_1 = ',num2str(pr(tn1)/pa)])
%
subplot(2,2,2)
plot(xv,rhov,'-*')
ylabel('\rho')
%
subplot(2,2,3)
plot(xv,vv,'-*')
ylabel('v')
%
subplot(2,2,4)
plot(xv,abs(vv)./av)
ylabel('M')
title(['t = ',num2str(tn1)])
drawnow
% disp([p1/pa,p2/pa,abs(vv(end)/av(end))])
if ~isreal(vv) | isnan(vv)
keyboard
end
end
end
function F = compF(W,gamma)
%Receives matrix W and computes F:
%Extract p,rho,v from W:
rho = W(1,:);
v = W(2,:)./rho;
p = (gamma-1)*(W(3,:)-0.5*rho.*v.^2);
%Compute F:
F = [rho.*v; p+rho.*v.^2; gamma/(gamma-1)*p.*v+0.5*rho.*v.^3];
end
function [p,rho,v] = compU(W,gamma)
%Receives matrix W and computes U:
%Extract p,rho,v from W:
rho = W(1,:);
v = W(2,:)./rho;
p = (gamma-1)*(W(3,:)-0.5*rho.*v.^2);
end
function f = fLeft(u,pLeft,h0Left,r3_1,gamma)
%Receives u=[p1;rho1;v1] and r3_(j=2) and imposes p1=pLeft, h01=hLeft and
%D_(u-a) w_3 = 0:
%Extract p,rho,v:
p1 = u(1);
rho1 = u(2);
v1 = u(3);
c1 = sqrt(gamma*p1/rho1);
%Compute r3:
r3 = v1 - 2/(gamma-1)*c1;
%Residual:
f = [ p1-pLeft;
gamma/(gamma-1)*p1/rho1 + v1^2/2 - h0Left;
r3-r3_1 ];
end
function f = fRight(u,pRight,r1_N,r2_N,gamma)
%Receives u=[pN;rhoN;vN] and imposes pN=pRight, D_u(w1)=0 and
%D_(u+a)(w2) = 0:
%Extract p,rho,v:
pN = u(1);
rhoN = u(2);
vN = u(3);
cN = sqrt(gamma*pN/rhoN);
%Compute r1 and r2:
r1 = pN/rhoN^gamma;
r2 = vN + 2/(gamma-1)*cN;
%Residual:
f = [ pN-pRight;
r1-r1_N;
r2-r2_N ];
end
function TubeProblemLWextrap
%Resolution by explicit LaxWendroff. Boundary conditions imposed by
%extrapolation of Riemann invariants.
%Not fully optimized.
close all
clc
%Points at the mesh and length of the tube:
N = 100;
L = 1;
xv = linspace(0,1,N);
Deltax = L/(N-1);
%Final t and relation Deltat/Deltatmax:
rDeltat = 0.95;
tfin = 0.05;
%Artificial viscosity:
D = 0;
%Solver options:
options = optimset('TolFun',1e-8,'TolX',1e-8,'display','off');
%Gas properties:
Rg = 287;
gamma = 1.4;
cp = gamma/(gamma-1)*Rg;
cv = 1/(gamma-1)*Rg;
%Boundary conditions:
pa = 1e5;
Ta = 300;
% pr = @(t) pa+0.1*pa*atan(t/0.001)/(pi/2);
pr = @(t) 1.2*pa;
Tr = 300;
%Initial conditions:
pv = pa*ones(1,N);
Tv = Ta*ones(1,N);
rhov = pv./Rg./Tv;
vv = zeros(1,N);
av = sqrt(gamma*pv./rhov);
%Initialize matrices containing the values of p,rho and v at each time:
pm = pv;
rhom = rhov;
vm = vv;
%Loop until tfin
tn = 0;
[Wn,Fn] = CompWF(reshape([pv;rhov;vv],[],1),gamma);
Wn = reshape(Wn,3,[]);
Fn = reshape(Fn,3,[]);
iter = 0;
while tn<=tfin
iter = iter+1;
%Compute Deltat:
Deltat = rDeltat*Deltax/max(abs(vv)+av);
tn1 = tn+Deltat;
%Add artificial viscosity to Fn:
Fn(:,2:end) = Fn(:,2:end) - D*Deltax*(Wn(:,2:end)-Wn(:,1:end-1));
%Prediction of W_(j+0.5)^(n+0.5) and F_(j+0.5)^(n+0.5):
Wpred = 0.5*(Wn(:,1:end-1)+Wn(:,2:end)) - ...
Deltat/2/Deltax*(Fn(:,2:end)-Fn(:,1:end-1));
Fpred = compF(Wpred,gamma);
%Add artificial viscosity to Fpred:
Fpred(:,1:end-1) = Fpred(:,1:end-1) - D*Deltax*(Wpred(:,2:end)-Wpred(:,1:end-1));
%Next W in the center points:
Wn1 = Wn(:,2:end-1) - Deltat/Deltax * (Fpred(:,2:end)-Fpred(:,1:end-1));
%Get p,rho,T at j=2,3 and t=t^(n+1):
[p23,rho23,v23] = compU(Wn1(:,[1 2]),gamma);
%Get third Riemann invariant at j=2,3 and t=t^(n+1):
[r3_23] = v23 - 2/(gamma-1)*sqrt(gamma*p23./rho23);
%Extrapolate third Riemann invariant at j=1 and t=t^(n+1):
r3_1 = r3_23*[2;-1];
%Get p,rho,u at the j=1 and t=t_n^1 imposing boundary conditions:
% u1sol = Broyden(@(u)fLeft(u,pr(tn1),cp*Tr,r3_1,gamma),...
% [p23(1);rho23(1);v23(1)],[],options,[0,0,-Inf],Inf(1,3));
% u1sol = NewtonRaphson(@(u)fLeft(u,pr(tn1),cp*Tr,r3_1,gamma),...
% [p23(1);rho23(1);v23(1)],options,[0,0,-Inf],Inf(1,3));
u1sol = fsolve(@(u)fLeft(u,pr(tn1),cp*Tr,r3_1,gamma),...
[p23(1);rho23(1);v23(1)],options);
if u1sol(3)<0
u1sol(3) = v23(1);
end
%Get p,rho,T at j=N-2,N-1 and t=t^(n+1):
[pN_21,rhoN_21,vN_21] = compU(Wn1(:,end-1:end),gamma);
%Get first and second Riemann invariants at j=N-2,N-1 ant t=t^(n+1):
r1_Nm21 = pN_21./rhoN_21.^gamma;
r2_Nm21 = vN_21 + 2/(gamma-1)*sqrt(gamma*pN_21./rhoN_21);
%Extrapolate Riemann invariants to j=N and t=t^(n+1):
r1_N = r1_Nm21*[-1;2];
r2_N = r2_Nm21*[-1;2];
%Get p,rho,u at j=N and t=t_n^1 imposing boundary conditions:
% uNsol = Broyden(@(u)fRight(u,pa,r1_N,r2_N,gamma),...
% [pN_21(2);rhoN_21(2);vN_21(2)],[],options,[0,0,-Inf],Inf(1,3));
% uNsol = NewtonRaphson(@(u)fRight(u,pa,r1_N,r2_N,gamma),...
% [pN_21(2);rhoN_21(2);vN_21(2)],options,[0,0,-Inf],Inf(1,3));
uNsol = fsolve(@(u)fRight(u,pa,r1_N,r2_N,gamma),...
[pN_21(2);rhoN_21(2);vN_21(2)],options);
if uNsol(3)<0
uNsol(3) = vN_21(end);
end
%Complete with W at boundaries:
Wn1 = cat(2,CompWF(u1sol,gamma),Wn1,CompWF(uNsol,gamma));
Fn1 = compF(Wn1,gamma);
%Extract new values of p, T and v:
[pv,rhov,vv] = compU(Wn1,gamma);
av = sqrt(gamma*pv/rhov);
pm = cat(2,pm,pv);
rhom = cat(2,rhom,rhov);
vm = cat(2,vm,vv);
%Actualize variables:
Wn = Wn1;
Fn = Fn1;
tn = tn1;
%plot variables:
figure(1)
subplot(2,2,1)
plot(xv,pv/pa,'-*')
ylabel('p')
% title(['p_1 = ',num2str(pr(tn1)/pa)])
%
subplot(2,2,2)
plot(xv,rhov,'-*')
ylabel('\rho')
%
subplot(2,2,3)
plot(xv,vv,'-*')
ylabel('v')
%
subplot(2,2,4)
plot(xv,abs(vv)./av)
ylabel('M')
title(['t = ',num2str(tn1)])
drawnow
% disp([p1/pa,p2/pa,abs(vv(end)/av(end))])
if ~isreal(vv) | isnan(vv)
keyboard
end
end
end
function F = compF(W,gamma)
%Receives matrix W and computes F:
%Extract p,rho,v from W:
rho = W(1,:);
v = W(2,:)./rho;
p = (gamma-1)*(W(3,:)-0.5*rho.*v.^2);
%Compute F:
F = [rho.*v; p+rho.*v.^2; gamma/(gamma-1)*p.*v+0.5*rho.*v.^3];
end
function [p,rho,v] = compU(W,gamma)
%Receives matrix W and computes U:
%Extract p,rho,v from W:
rho = W(1,:);
v = W(2,:)./rho;
p = (gamma-1)*(W(3,:)-0.5*rho.*v.^2);
end
function f = fLeft(u,pLeft,h0Left,r3_1,gamma)
%Receives u=[p1;rho1;v1] and r3_(j=2) and imposes p1=pLeft, h01=hLeft and
%D_(u-a) w_3 = 0:
%Extract p,rho,v:
p1 = u(1);
rho1 = u(2);
v1 = u(3);
c1 = sqrt(gamma*p1/rho1);
%Compute r3:
r3 = v1 - 2/(gamma-1)*c1;
%Residual:
f = [ p1-pLeft;
gamma/(gamma-1)*p1/rho1 + v1^2/2 - h0Left;
r3-r3_1 ];
end
function f = fRight(u,pRight,r1_N,r2_N,gamma)
%Receives u=[pN;rhoN;vN] and imposes pN=pRight, D_u(w1)=0 and
%D_(u+a)(w2) = 0:
%Extract p,rho,v:
pN = u(1);
rhoN = u(2);
vN = u(3);
cN = sqrt(gamma*pN/rhoN);
%Compute r1 and r2:
r1 = pN/rhoN^gamma;
r2 = vN + 2/(gamma-1)*cN;
%Residual:
f = [ pN-pRight;
r1-r1_N;
r2-r2_N ];
end
function TubeProblemRoeextrap
%Resolution by explicit Roe. Boundary conditions imposed by
%extrapolation of Riemann invariants.
%Not fully optimized.
close all
clc
%Points at the mesh and length of the tube:
N = 100;
L = 1;
xv = linspace(0,1,N);
Deltax = L/(N-1);
%Final t and relation Deltat/Deltatmax:
rDeltat = 0.95;
tfin = 0.05;
%Solver options:
options = optimset('TolFun',1e-8,'TolX',1e-8,'display','off');
%Gas properties:
Rg = 287;
gamma = 1.4;
cp = gamma/(gamma-1)*Rg;
cv = 1/(gamma-1)*Rg;
%Boundary conditions:
pa = 1e5;
Ta = 300;
% pr = @(t) pa+0.1*pa*atan(t/0.001)/(pi/2);
pr = @(t) 1.2*pa;
Tr = 300;
%Initial conditions:
pv = pa*ones(1,N);
Tv = Ta*ones(1,N);
rhov = pv./Rg./Tv;
vv = zeros(1,N);
av = sqrt(gamma*pv./rhov);
%Initialize matrices containing the values of p,rho and v at each time:
pm = pv;
rhom = rhov;
vm = vv;
tv = [];
%Loop until tfin
tn = 0;
[Wn,Fn] = CompWF(reshape([pv;rhov;vv],[],1),gamma);
Wn = reshape(Wn,3,[]);
Fn = reshape(Fn,3,[]);
iter = 0;
while tn<=tfin
iter = iter+1;
%Compute Deltat:
Deltat = rDeltat*Deltax/max(abs(vv)+av);
tn1 = tn+Deltat;
%Compute rhob, ub, Hb, cb:
Hv = gamma/(gamma-1)*pv./rhov + 0.5*vv.^2;
Rv = sqrt(rhov(2:end)./rhov(1:end-1));
rhobv = Rv.*rhov(1:end-1);
ubv = (Rv.*vv(2:end) + vv(1:end-1))./(Rv + 1);
Hbv = (Rv.*Hv(2:end) + Hv(1:end-1))./(Rv + 1);
cbv = sqrt((gamma-1)*(Hbv-ubv.^2/2));
%Valores de alfa_i (en un vector 1*Nx):
alfa1v = rhov(2:end) - rhov(1:end-1) - (pv(2:end)-pv(1:end-1))./cbv.^2;
alfa2v = vv(2:end) - vv(1:end-1) + (pv(2:end)-pv(1:end-1))./rhobv./cbv;
alfa3v = vv(2:end) - vv(1:end-1) - (pv(2:end)-pv(1:end-1))./rhobv./cbv;
%Autovalores (en un vector 1*Nx):
lambda1v = ubv;
lambda2v = ubv + cbv;
lambda3v = ubv - cbv;
%Autovectores (en una matriz 3*Nx; cada columna es el autovector de la
%celda correspondiente):
V1m = [ones(1,N-1); ubv; ubv.^2/2];
aux = rhobv/2./cbv;
V2m = [aux; aux.*(ubv+cbv); aux.*(Hbv+ubv.*cbv)];
V3m = -[aux; aux.*(ubv-cbv); aux.*(Hbv-ubv.*cbv)];
%Valores de F_(i+1/2) (en matriz 3*(N-1) llamada F05m):
F05m = 0.5*(Fn(:,1:end-1) + Fn(:,2:end)) - 0.5*V1m*spdiags(alfa1v.'.*abs(lambda1v.'),0,N-1,N-1) - ...
0.5*V2m*spdiags(alfa2v.'.*abs(lambda2v.'),0,N-1,N-1) - ...
0.5*V3m*spdiags(alfa3v.'.*abs(lambda3v.'),0,N-1,N-1);
%Next W in the center points:
Wn1 = Wn(:,2:end-1) - Deltat/Deltax * (F05m(:,2:end)-F05m(:,1:end-1));
%Get p,rho,T at j=2,3 and t=t^(n+1):
[p23,rho23,v23] = compU(Wn1(:,[1 2]),gamma);
%Get third Riemann invariant at j=2,3 and t=t^(n+1):
[r3_23] = v23 - 2/(gamma-1)*sqrt(gamma*p23./rho23);
%Extrapolate third Riemann invariant at j=1 and t=t^(n+1):
r3_1 = r3_23*[2;-1];
%Get p,rho,u at the j=1 and t=t_n^1 imposing boundary conditions:
% u1sol = Broyden(@(u)fLeft(u,pr(tn1),cp*Tr,r3_1,gamma),...
% [p23(1);rho23(1);v23(1)],[],options,[0,0,-Inf],Inf(1,3));
% u1sol = NewtonRaphson(@(u)fLeft(u,pr(tn1),cp*Tr,r3_1,gamma),...
% [p23(1);rho23(1);v23(1)],options,[0,0,-Inf],Inf(1,3));
u1sol = fsolve(@(u)fLeft(u,pr(tn1),cp*Tr,r3_1,gamma),...
[p23(1);rho23(1);v23(1)],options);
if u1sol(3)<0
u1sol(3) = v23(1);
end
%Get p,rho,T at j=N-2,N-1 and t=t^(n+1):
[pN_21,rhoN_21,vN_21] = compU(Wn1(:,end-1:end),gamma);
%Get first and second Riemann invariants at j=N-2,N-1 ant t=t^(n+1):
r1_Nm21 = pN_21./rhoN_21.^gamma;
r2_Nm21 = vN_21 + 2/(gamma-1)*sqrt(gamma*pN_21./rhoN_21);
%Extrapolate Riemann invariants to j=N and t=t^(n+1):
r1_N = r1_Nm21*[-1;2];
r2_N = r2_Nm21*[-1;2];
%Get p,rho,u at j=N and t=t_n^1 imposing boundary conditions:
% uNsol = Broyden(@(u)fRight(u,pa,r1_N,r2_N,gamma),...
% [pN_21(2);rhoN_21(2);vN_21(2)],[],options,[0,0,-Inf],Inf(1,3));
% uNsol = NewtonRaphson(@(u)fRight(u,pa,r1_N,r2_N,gamma),...
% [pN_21(2);rhoN_21(2);vN_21(2)],options,[0,0,-Inf],Inf(1,3));
uNsol = fsolve(@(u)fRight(u,pa,r1_N,r2_N,gamma),...
[pN_21(2);rhoN_21(2);vN_21(2)],options);
if uNsol(3)<0
uNsol(3) = vN_21(end);
end
%Complete with W at boundaries:
Wn1 = cat(2,CompWF(u1sol,gamma),Wn1,CompWF(uNsol,gamma));
Fn1 = compF(Wn1,gamma);
%Extract new values of p, T and v:
[pv,rhov,vv] = compU(Wn1,gamma);
av = sqrt(gamma*pv/rhov);
pm = cat(1,pm,pv);
rhom = cat(1,rhom,rhov);
vm = cat(1,vm,vv);
tv = cat(1,tv,tn1);
%Actualize variables:
Wn = Wn1;
Fn = Fn1;
tn = tn1;
%plot variables:
figure(1)
subplot(2,2,1)
plot(xv,pv/pa,'-*')
ylabel('p')
% title(['p_1 = ',num2str(pr(tn1)/pa)])
%
subplot(2,2,2)
plot(xv,rhov,'-*')
ylabel('\rho')
%
subplot(2,2,3)
plot(xv,vv,'-*')
ylabel('v')
%
subplot(2,2,4)
plot(xv,abs(vv)./av)
ylabel('M')
title(['t = ',num2str(tn1)])
drawnow
% disp([p1/pa,p2/pa,abs(vv(end)/av(end))])
if ~isreal(vv) | isnan(vv)
keyboard
end
end
end
function F = compF(W,gamma)
%Receives matrix W and computes F:
%Extract p,rho,v from W:
rho = W(1,:);
v = W(2,:)./rho;
p = (gamma-1)*(W(3,:)-0.5*rho.*v.^2);
%Compute F:
F = [rho.*v; p+rho.*v.^2; gamma/(gamma-1)*p.*v+0.5*rho.*v.^3];
end
function [p,rho,v] = compU(W,gamma)
%Receives matrix W and computes U:
%Extract p,rho,v from W:
rho = W(1,:);
v = W(2,:)./rho;
p = (gamma-1)*(W(3,:)-0.5*rho.*v.^2);
end
function f = fLeft(u,pLeft,h0Left,r3_1,gamma)
%Receives u=[p1;rho1;v1] and r3_(j=2) and imposes p1=pLeft, h01=hLeft and
%D_(u-a) w_3 = 0:
%Extract p,rho,v:
p1 = u(1);
rho1 = u(2);
v1 = u(3);
c1 = sqrt(gamma*p1/rho1);
%Compute r3:
r3 = v1 - 2/(gamma-1)*c1;
%Residual:
f = [ p1-pLeft;
gamma/(gamma-1)*p1/rho1 + v1^2/2 - h0Left;
r3-r3_1 ];
end
function f = fRight(u,pRight,r1_N,r2_N,gamma)
%Receives u=[pN;rhoN;vN] and imposes pN=pRight, D_u(w1)=0 and
%D_(u+a)(w2) = 0:
%Extract p,rho,v:
pN = u(1);
rhoN = u(2);
vN = u(3);
cN = sqrt(gamma*pN/rhoN);
%Compute r1 and r2:
r1 = pN/rhoN^gamma;
r2 = vN + 2/(gamma-1)*cN;
%Residual:
f = [ pN-pRight;
r1-r1_N;
r2-r2_N ];
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.