MMPP | Ústav technické matematiky

Matematické modelování problémů proudění

kód předmětu: 2011083

Obsah

Plán předmětu
Seznam doporučené literatury
Materiály k předmětu

Přednášející:

Doc. Ing. Jiří Fürst, PhD.

Ústav technické matematiky, FS ČVUT v Praze

Karlovo nám. 13

121 35 Praha 2

místnost D301

Plán předmětu

Jednorozměrný zákon zachování
- Rovnice kontinuity v 1D
- Lineární advekce
  - Charakteristiky, analytické řešení počáteční úlohy
  - Formulace smíšené úlohy, okrajové podmínky
- Nelineární skalární hyperbolická rovnice (Burgersova rovnice)
  - Vznik nespojitosti a neexistence globálního klasického řešení
  - Definice slabého řešení a jeho nejednoznačnost
  - Koncept entropického řešení
  - Rankine-Hugoniotova podmínka (rychlost pohybu nespojitosti)
- Další příklady nelineárních rovnic
- Lineární hyperbolický systém rovnic prvního řádu
  - Diagonalizace matice A, charakteristické proměnné
  - Analytické řěšení počáteční úlohy
  - Formulace smíšené úlohy, okrajové podmínky
Numerické řešení jednorozměrných skalárních zákonů zachování
- Metoda konečných objemů pro 1D problém
- Chyba aproximace, konzistence, stabilita, globální chyba
- Laxova věta o ekvivalenci
- Spektrální kriterium stability
- Monotónní schéma a princip maxima pro skalární problém
- Princip umělé vazkosti
- TVD schéma pro skalární problém
- Příklady schémat:
  - Centrální schéma
  - Laxovo-Wendroffovo schéma
  - Schéma upwind
  - Laxovo-Friedrichsovo (Rusanovovo) schéma
  - Laxovo-Wendroffovo schéma s umělou vazkostí
  - TVD schéma
Numerické řešení hyperbolických systémů v 1D
- Schéma upwind
- Laxovo-Freidrichsovo schéma
- Laxovo-Wendroffovo (MacCormackovo) schéma
- AUSM schéma pro Eulerovy rovnice
- Příklady:
  - Řešení Riemannova problému pro systém Eulerových rovnic (rázová trubice)
  - Řešení rovnic proudění mělké vody (protržení hráze)
Princip metody konečných objemů pro vícerozměrné problémy

Plán bude doplněn v průběhu semestru

Seznam doporučené literatury

J. Fürst, K. Kozel: Numerické metody řešení problémů proudění I, skripta ČVUT, 2001
J. Fořt, K. Kozel: Numerické metody řešení problémů proudění II, skripta ČVUT, 2002
J. Fořt, K. Kozel, P. Louda, J. Fürst: Numerické metody řešení problémů proudění III, skripta ČVUT, 2004
J.H. Ferziger, M. Peric: Computational Methods for Fluid Dynamics, Springer 2002
J. Blazek: Computational Fluid Dynamics: Principles and Applications, Elsevier, 2001
E.F.Toro: Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Springer 2009

Materiály k předmětu

Program pro numerické řešení smíšené úlohy pro rovnici linární advekce v 1D,

schéma upwind

clear all;

% Reseni smisene ulohy pro rovnici u_t + a*u_x = 0 pomoci MKO
a = 1;

% Delka intervalu a pocet bunek
L = 1;
n = 100;
dx = L / n;
x = linspace(dx/2, L-dx/2, n);
% x = dx/2:dx:(L-dx/2);

u(1:n)  = 0;    % aktualni hodnota reseni
un(1:n) = 0;    % nova hodnota reseni

% Pocatecni podminka
for i = 1:n
  if (x(i)>0.25 && x(i)<0.5)
    u(i) = 1-cos(2*pi*(x(i)/0.25-1));
  end
end

plot(x, u); axis([0 1 -0.5 2.5]);
disp("Stiskni enter pro pokracovani"); pause;

dt = 0.8*dx;
t = 0;

for iter = 1:(n/2)
  % Okrajova podminka na leve casti intervalu
  f(1) = 0;

  % Numericky tok f
  for i = 2:n+1
    f(i) = a*u(i-1);
  end

  un(1:n) = u(1:n) - dt/dx * (f(2:n+1) - f(1:n));
  u = un;
  t = t + dt;

  if (mod(iter,10)==0)
    plot(x, u); axis([0 1 -0.5 2.5]);
    pause(1);
  end
end

%clear all;

% Reseni smisene ulohy pro rovnici u_t + a*u_x = 0 pomoci MKO
a = 1;

% Delka intervalu a pocet bunek
L = 1;
n = 100;

dx = L / n;
x = linspace(dx/2, L-dx/2, n);

u(1:n)  = 0;
un(1:n) = 0;

% Pocatecni podminka
for i = 1:n
  if (x(i)>0.25 && x(i)<0.5)
     u(i) = 1-cos(2*pi*(x(i)/0.25-1));
  end
end

plot(x, u); axis([0 1 -0.5 2.5]);
disp("Stiskni enter pro pokracovani"); pause;

dt = 0.8*dx;

for iter = 1:(n/2)

    % Okrajova podminka na leve casti intervalu
    f(1) = 0;
    f(n+1) = a*u(n);

    % Numericky tok f
    for i = 2:n
        f(i) = a/2*(u(i-1)+u(i)) - a^2/2*dt/dx*(u(i)-u(i-1));
    end

    un(1:n) = u(1:n) - dt/dx * (f(2:n+1) - f(1:n));

    u = un;

    if (mod(iter,10)==0)
      plot(x, u); axis([0 1 -0.5 2.5]);
      pause(1);
    end
end

Program pro numerické řešení smíšené úlohy pro Burgersovu rovnici v 1D,

schéma upwind

clear all;

% Reseni smisene ulohy pro rovnici u_t + (u^2/2)_x = 0 pomoci MKO

% Delka intervalu a pocet bunek
L = 1;
n = 100;

dx = L / n;
x = linspace(dx/2, L-dx/2, n);

u(1:n)  = 0;
un(1:n) = 0;

% Pocatecni podminka
for i = 1:n
  if (x(i)>0.25 && x(i)<0.5)
     u(i) = 1-cos(2*pi*(x(i)/0.25-1));
  end
end

plot(x, u); axis([0 1 -0.5 2.5]);
disp("Stiskni enter pro pokracovani"); pause;

dt = 0.4*dx;

t = 0;

for iter = 1:n

    % Okrajova podminka na leve casti intervalu
    f(1) = 0;

    % Numericky tok f
    for i = 2:n+1
        f(i) = u(i-1)^2/2;
    end

    un(1:n) = u(1:n) - dt/dx * (f(2:n+1) - f(1:n));
    t = t + dt;

    u = un;

    if (mod(iter,10)==0)
      plot(x, u); axis([0 1 -0.5 2.5]);
      pause(1);
    end
end

Laxovo-Wendroffovo schéma

clear all;

% Reseni smisene ulohy pro rovnici u_t + (u^2/2)_x = 0 pomoci MKO

% Delka intervalu a pocet bunek
L = 1;
n = 400;

dx = L / n;
x = linspace(dx/2, L-dx/2, n);

u(1:n)  = 0;
un(1:n) = 0;

% Pocatecni podminka
for i = 1:n
  if (x(i)>0.25 && x(i)<0.5)
     u(i) = 1-cos(2*pi*(x(i)/0.25-1));
  end
end

plot(x, u); axis([0 1 -0.5 2.5]);
disp("Stiskni enter pro pokracovani"); pause;

dt = 0.1*dx;

t = 0;

for iter = 1:n

    % Okrajova podminka na leve casti intervalu
    f(1) = 0;
    f(n+1) = u(n)^2/2;

    % Numericky tok f
    for i = 2:n
        a = (u(i-1)+u(i))/2;
        f(i) = a/2*(u(i-1)+u(i)) - a^2/2*dt/dx*(u(i)-u(i-1));
    end

    un(1:n) = u(1:n) - dt/dx * (f(2:n+1) - f(1:n));
    t = t + dt;

    u = un;

    if (mod(iter,10)==0)
      plot(x, u); axis([0 1 -0.5 2.5]);
      pause(1);
    end
end

TVD schéma

%clear all;

% Reseni smisene ulohy pro rovnici u_t + a*u_x = 0 pomoci MKO
a = 1;

% Delka intervalu a pocet bunek
L = 1;
n = 100;

dx = L / n;
x = linspace(dx/2, L-dx/2, n);

u(1:n)  = 0;
un(1:n) = 0;

% Pocatecni podminka
for i = 1:n
  if (x(i)>0.25 && x(i)<0.5)
     u(i) = 1-cos(2*pi*(x(i)/0.25-1));
  end
end

plot(x, u); axis([0 1 -0.5 2.5]);
disp("Stiskni enter pro pokracovani"); pause;

dt = 0.8*dx;

for iter = 1:(n/2)

    % Okrajova podminka na leve casti intervalu
    f(1) = 0;
    f(2) = 0;
    f(n+1) = a*u(n-1);
    f(n+1) = a*u(n);

    % Numericky tok f
    for i = 3:n-1
	fl = a*u(i-1);
	fr = a*u(i);
	flw = (fr+fl)/2 - a/2*dt/dx*(fr-fl);
	fup = (fr+fl)/2 - abs(a)/2*(u(i)-u(i-1));

	if ( (u(i)-u(i-1))*(u(i+1)-u(i)) <= 0)
	   phi = 0;
	else
	  r = (u(i-1)-u(i-2)) / (u(i)-u(i-1));
	  phi = min(1, 2*r);
	end
	f(i) = fup + phi*(flw-fup);
    end

    un(1:n) = u(1:n) - dt/dx * (f(2:n+1) - f(1:n));

    u = un;

    if (mod(iter,10)==0)
      plot(x, u); axis([0 1 -0.5 2.5]);
      pause(1);
    end
end


clear all;

% Reseni smisene ulohy pro rovnici u_t + (u^2/2)_x = 0 pomoci MKO

% Delka intervalu a pocet bunek
L = 1;
n = 100;

dx = L / n;
x = linspace(dx/2, L-dx/2, n);

u(1:n)  = 0;
un(1:n) = 0;

% Pocatecni podminka
for i = 1:n
if (x(i)>0.25 && x(i)<0.5)
u(i) = 1-cos(2*pi*(x(i)/0.25-1));
end
end

plot(x, u); axis([0 1 -0.5 2.5]);
disp("Stiskni enter pro pokracovani"); pause;

dt = 0.1*dx;

t = 0;

for iter = 1:n

% Okrajova podminka na leve casti intervalu
f(1) = 0;
f(2) = 0;
f(n) = 0;
f(n+1) = 0;

% Numericky tok f
for i = 2:n-1
a = (u(i-1)+u(i))/2;
fl = u(i-1)^2/2;
fr = u(i)^2/2;

flw = (fr+fl)/2 - a/2*dt/dx*(fr-fl);
fup = (fr+fl)/2 - abs(a)/2*(u(i)-u(i-1));

if ( (u(i)-u(i-1))*(u(i+1)-u(i)) <= 0)
phi = 0;
else
rl = (u(i-1)-u(i-2)) / (u(i)-u(i-1));
rr = (u(i+1)-u(i)) / (u(i)-u(i-1));
if (rr<0)
phi = 0;
else
phi = min(1, 2*min(rl, rr));
end
end
f(i) = fup + phi*(flw-fup);
end

un(1:n) = u(1:n) - dt/dx * (f(2:n+1) - f(1:n));
t = t + dt;

u = un;

if (mod(iter,10)==0)
plot(x, u); axis([0 1 -0.5 2.5]);
pause(1);
end
end

Programy pro řešení rovnice mělké vody

Rusanovovo schéma

%
%   rovnice resime ve tvaru W_t + F(W)_x = 0, kde
%   W = [h, q], F(W) = [q, q^2/h + gh^2/2] kde q = hu
%

g = 10;

% Delka intervalu a pocet bunek
L = 1;
n = 400;

dx = L / n;
x = dx/2:dx:(L-dx/2);

W(1:2,1:n)  = 0;

% Pocatecni podminka
for i = 1:n
  if (x(i)<0.5)
    W(1,i) = 1.0;
  else
    W(1,i) = 0.01;
  end
end

plot(x, W(1,:)); axis([0 1 0 1.5]);
disp("Stiskni enter pro pokracovani"); pause;
t = 0;

for iter = 1:n

    h = W(1,:);
    q = W(2,:);
    u = q ./ h;

    % Vypocet fyzikalniho toku
    FF(1,:) = q;
    FF(2,:) = q.^2./h + g*h.^2/2;

    % Vypocet spektralniho polomeru jakobianu
    sigma = abs(u) + sqrt(g*h);
    %

    dt = 0.4 * dx / max(sigma);

    % Numericky tok na leve a prave casti intervalu
    F(:,1) = FF(:,1);
    F(:,n+1) = FF(:,n);

    % Numericky tok f
    for i = 2:n
	s = max(sigma(i-1),sigma(i));
	F(:,i) = (FF(:,i-1)+FF(:,i))/2 - s/2*(W(:,i)-W(:,i-1));
    end

    W(:,1:n) = W(:,1:n) - dt/dx * (F(:,2:n+1) - F(:,1:n));
    t = t + dt;

    if (mod(iter,10)==0)
      plot(x, W(1,:)); axis([0 1 0 1.5]);
      pause(1);
    end
end

HLL schéma s rekonstrukcí


clear all;
% Reseni pocatecni ulohy pro rovnice melke vody
%   h_t + (uh)_x = 0
%   (hu)_t + (hu^2 + gh^2/2)_x = 0
%
%   rovnice resime ve tvaru W_t + F(W)_x = 0, kde
%   W = [h, q], F(W) = [q, q^2/h + gh^2/2] kde q = hu
%

g = 10;

% Delka intervalu a pocet bunek
L = 1;
n = 400;

dx = L / n;
x = dx/2:dx:(L-dx/2);

W(1:2,1:n)  = 0;

% Pocatecni podminka
for i = 1:n
if (x(i)<0.5)
W(1,i) = 1.0;
else
W(1,i) = 0.01;
end
end

plot(x, W(1,:)); axis([0 1 0 1.5]);
disp("Stiskni enter pro pokracovani"); pause;
t = 0;

for iter = 1:n

% Vypocet "derivace" W
DW(1:2,1) = 0;
DW(1:2,n) = 0;
for i = 2:n-1
for k = 1:2
if ( (W(k,i)-W(k,i-1))*(W(k,i+1)-W(k,i)) <= 0 )
DW(k,i) = 0;
else
if ( abs(W(k,i)-W(k,i-1)) < abs(W(k,i+1)-W(k,i)) )
DW(k,i) = W(k,i)-W(k,i-1);
else
DW(k,i) = W(k,i+1)-W(k,i);
end
end
end
end

h = W(1,:);
q = W(2,:);
u = q ./ h;

% Vypocet spektralniho polomeru jakobianu
sigma = abs(u) + sqrt(g*h);
%

dt = 0.4 * dx / max(sigma);

% Numericky tok na leve a prave casti intervalu
F(:,1) = [q(1); q(1)^2/h(1) + g*h(1)^2/2] ;
F(:,n+1) = [q(n); q(n)^2/h(n) + g*h(n)^2/2] ;

% Numericky tok f
for i = 2:n
sl = min( u(i-1)-sqrt(g*h(i-1)), u(i)-sqrt(g*h(i)) );
sr = min( u(i-1)+sqrt(g*h(i-1)), u(i)+sqrt(g*h(i)) );

Wl = W(:,i-1) + DW(:,i-1)/2;
Wr = W(:,i) - DW(:,i)/2;

Fl = [Wl(2); Wl(2)^2/Wl(1) + g*Wl(1)^2/2];
Fr = [Wr(2); Wr(2)^2/Wr(1) + g*Wr(1)^2/2];

if (0<=sl)
F(:,i) = Fl;
elseif (0<sr)
F(:,i) = (sr*Fl - sl*Fr + sl*sr*(Wr-Wl)) / (sr-sl);
else
F(:,i) = Fr;
end

end

W(:,1:n) = W(:,1:n) - dt/dx * (F(:,2:n+1) - F(:,1:n));
t = t + dt;

if (mod(iter,10)==0)
plot(x, W(1,:)); axis([0 1 0 1.5]);
pause(1);
end
end

</pre>
clear all;

% Reseni proudeni melke vody pres prekazku
%
%   h_t + (uh)_x = 0
%   (hu)_t + (hu^2 + gh^2/2)_x = -g h b'
%
%   rovnice resime ve tvaru W_t + F(W)_x = Q, kde
%   W = [h, q], F(W) = [q, q^2/h + gh^2/2] kde q = hu a
%   Q = [0, -g h b']
%
%
%  Jakobian dF/dW = [ 0, 1; gh^2-u, 2u]
%  vlastni cisla  u +- sqrt(g*h)
%  vlastni vektory [1; u+-sqrt(g*h)]
%
%  charakteristicke promenne V = [u+2*c,u-2*c]     c =sqrt(g*h)

g = 10;

% Rychlost vtoku m3/s
q_in = 1;
h_out = 0.8;

% Delka intervalu a pocet bunek
L = 1;
n = 200;

dx = L / n;
x = dx/2:dx:(L-dx/2);

% tvar dna
for i = 1:n
if (0.3 < x(i) && x(i)<0.7)
b(i) = (cos(pi*(x(i)-0.5)/0.2) + 1)/2 + (1-x(i))*0.2;
bx(i) = ( (cos(pi*(x(i)+dx/2-0.5)/0.2) + 1)/2 -
(cos(pi*(x(i)-dx/2-0.5)/0.2) + 1)/2 ) / dx - 0.2;
else
b(i) = (1-x(i))*0.2;
bx(i)= -0.2;
end
end

W(1,1:n)  = max(h_out,max(b)+0.1) - b;
W(2,1:n)  = 0;

plot(x, b, x, W(1,:)+b); axis([0 1 0 3]);
disp("Stiskni enter pro pokracovani"); pause;
t = 0;

tend = 1;

Q(1:2,1:n) = 0;

for iter = 1:(40*n)

h = W(1,:);
q = W(2,:);
u = q ./ h;

% Vypocet spektralniho polomeru jakobianu
sigma = abs(u) + sqrt(g*h);
%

dt = 0.4 * dx / max(sigma);

% Numericky tok na vstupu (zadany prutok)
F(:,1) = [q_in; q_in^2/h(1) + g*h(1)^2/2] ;

% Numericky tok na vystupu (zadana vyska hladiny)
F(:,n+1) = [q(n); q(n)^2/h_out + g*h_out^2/2] ;

% Numericky tok f
for i = 2:n
sl = min( u(i-1)-sqrt(g*h(i-1)), u(i)-sqrt(g*h(i)) );
sr = min( u(i-1)+sqrt(g*h(i-1)), u(i)+sqrt(g*h(i)) );

Wl = W(:,i-1);
Wr = W(:,i);

Fl = [Wl(2); Wl(2)^2/Wl(1) + g*Wl(1)^2/2];
Fr = [Wr(2); Wr(2)^2/Wr(1) + g*Wr(1)^2/2];

if (0<=sl)
F(:,i) = Fl;
elseif (0<sr)
F(:,i) = (sr*Fl - sl*Fr + sl*sr*(Wr-Wl)) / (sr-sl);
else
F(:,i) = Fr;
end

end

% Zdrojovy clen
Q(2,1:n) = - g * h .* bx;

W(:,1:n) = W(:,1:n) - dt/dx * (F(:,2:n+1) - F(:,1:n)) + dt*Q;
t = t + dt;

if (mod(iter,10)==0)
plot(x, b, x, W(1,:)+b); axis([0 1 0 3]);
pause(0.1);
end
end

Řešení Eulerových rovnic v 1D


clear all;

% Reseni pocatecni ulohy pro Eulerovy rovnice
%   rho_t + (rho u)_x = 0
%   (rho u)_t + (rho u^2 + p)_x = 0
%   (rho E)_t + [u (rho E + p)]_x = 0
%
%   p = (kappa-1) * rho * (E - u^2/2)
%
%   rovnice resime ve tvaru W_t + F(W)_x = 0, kde
%   W = [rho, m, e], F(W) = [rho u, rho u^2 + p, u (rho E + p)]
%   kde m = rho u, e = rho E
%

kappa = 1.4;

% Delka intervalu a pocet bunek
L = 1;
n = 400;

dx = L / n;
x = dx/2:dx:(L-dx/2);

W(1:3,1:n)  = 0;

% Pocatecni podminka
for i = 1:n
if (x(i)<0.5)
W(1:3,i) = [ 10; 0; 1e6/(kappa-1) ];
else
W(1:3,i) = [ 1; 0; 1e5/(kappa-1) ];
end
end

plot(x, W(1,:)); axis([0 1 0 20]);
disp("Stiskni enter pro pokracovani"); pause;
t = 0;

for iter = 1:n

rho = W(1,:);
u   = W(2,:) ./ rho;
E   = W(3,:) ./ rho;

p = (kappa-1) * rho .* (E - 0.5*u.*u);
a = sqrt(kappa * p ./ rho);

FF(1,:) = rho .* u;
FF(2,:) = rho .* u.^2 + p;
FF(3,:) = u .* (rho .* E + p);

% Vypocet spektralniho polomeru jakobianu
sigma = abs(u) + a;
%

dt = 0.4 * dx / max(sigma);

% Numericky tok na leve a prave casti intervalu
F(:,1) = FF(:,1);
F(:,n+1) = FF(:,n);

% Numericky tok f
for i = 2:n
sl = min( u(i-1)-a(i-1), u(i)-a(i) );
sr = max( u(i-1)+a(i-1), u(i)+a(i) );

Wl = W(:,i-1);
Wr = W(:,i);

Fl = FF(:,i-1);
Fr = FF(:,i);

if (0<=sl)
F(:,i) = Fl;
elseif (0<sr)
F(:,i) = (sr*Fl - sl*Fr + sl*sr*(Wr-Wl)) / (sr-sl);
else
F(:,i) = Fr;
end

end

W(:,1:n) = W(:,1:n) - dt/dx * (F(:,2:n+1) - F(:,1:n));
t = t + dt;

if (mod(iter,10)==0)
plot(x, W(1,:)); axis([0 1 0 20]);
pause(1);
end
end

clear all;

% Reseni proudeni 1D dyzou
 %   (A rho)_t + (A rho u)_x = 0
 %   (A rho u)_t + ( A rho u^2 + A p)_x = p A_x
 %   (A rho E)_t + [A u (rho E + p)]_x = 0
 %
 %   p = (kappa-1) * rho * (E - u^2/2)
 %
 %   rovnice resime ve tvaru A W_t + (A F(W))_x = Q, kde
 %   W = [rho, m, e], F(W) = [rho u, rho u^2 + p, u (rho E + p)]
 %   kde m = rho u, e = rho E
 %
 %   na vstupu zadvame celkovy tlak a celkovou teplotu
 %   na vystupu zadavame staticky tlak
 %
 kappa = 1.4;

% Delka intervalu a pocet bunek
 L = 1;
 n = 200;

dx = L / n;
 x = dx/2:dx:(L-dx/2);

% Prurezy dyzy na hranicich bunek  A(x) = 1 + (x-0.5)^2
 A(1:n+1) = ((0:dx:L)-L/2) .^2 + 1;
 Ai(1:n) = (A(1:n) + A(2:n+1)) / 2;

W(1:3,1:n) = 0;
 Q(1:3,1:n) = 0;

% Pocatecni podminka
 for i = 1:n
 W(1:3,i) = [ 1.1; 0; 1e5/(kappa-1) ];
 end

plot(x, W(1,:)); axis([0 1 0 1.2]);
 disp("Stiskni enter pro pokracovani"); pause;
 t = 0;

Rezidua = [];

for iter = 1:20000

rho = W(1,:);
 u   = W(2,:) ./ rho;
 E   = W(3,:) ./ rho;

p = (kappa-1) * rho .* (E - 0.5*u.*u);
 a = sqrt(kappa * p ./ rho);

FF(1,:) = rho .* u;
 FF(2,:) = rho .* u.^2 + p;
 FF(3,:) = u .* (rho .* E + p);

% Vypocet spektralniho polomeru jakobianu
 sigma = abs(u) + a;
 %

dt = 0.4 * dx / max(sigma);

% Numericky tok na leve a prave casti intervalu
 pTot = 1.e5;
 rhoTot = 1.2;

pb =  min(p(1),pTot);
 M2 = ( (pb/pTot) ^ ((1-kappa)/kappa) - 1 ) * 2/(kappa-1);
 rhob = rhoTot * (1+(kappa-1)/2*M2) ^ (1/(1-kappa));
 cb = sqrt(kappa*pb/rhob);
 ub = cb*sqrt(M2);
 rEb = pb/(kappa-1) + 0.5*rhob*ub^2;
 F(:,1) = [rhob*ub; rhob*ub^2+pb; ub*(rEb + pb)];

p2 = 0.75 * pTot;

pb = p2;
 rhob = rho(n);
 ub = u(n);
 rEb = pb/(kappa-1) + 0.5*rhob*ub^2;
 F(:,n+1) = [rhob*ub; rhob*ub^2+pb; ub*(rEb + pb)];

% Numericky tok f
 for i = 2:n
 sl = min( u(i-1)-a(i-1), u(i)-a(i) );
 sr = max( u(i-1)+a(i-1), u(i)+a(i) );

Wl = W(:,i-1);
 Wr = W(:,i);

Fl = FF(:,i-1);
 Fr = FF(:,i);

if (0<=sl)
 F(:,i) = Fl;
 elseif (0<sr)
 F(:,i) = (sr*Fl - sl*Fr + sl*sr*(Wr-Wl)) / (sr-sl);
 else
 F(:,i) = Fr;
 end

end
 Q(2,:) = p .* (A(2:n+1)-A(1:n))/dx;

for k=1:3
 F(k,:) = A .* F(k,:);
 Wn(k,1:n) = W(k,1:n) - dt/dx./Ai .* (F(k,2:n+1) - F(k,1:n)) + dt./Ai.*Q(k,:);
 end

if (mod(iter,100)==0)
 Rezidua = [ Rezidua,  sqrt(sum((Wn-W).**2,dim=2)) / dt];
 plot(x, p/1e5); axis([0 1 0 1.2]);
 pause(0.01);
 end

W = Wn;
 t = t + dt;

end

Obsah

Plán předmětu

Seznam doporučené literatury

Materiály k předmětu

Program pro numerické řešení smíšené úlohy pro rovnici linární advekce v 1D,

Program pro numerické řešení smíšené úlohy pro Burgersovu rovnici v 1D,

Programy pro řešení rovnice mělké vody

Řešení Eulerových rovnic v 1D

Základní předměty

Další předměty

Konzultace

Přihlásit

Uživatelské jméno
Heslo

	Pamatuj si mě Zapomněli jste heslo?