I have written a Roe solver with Harten entropy fix code in Matlab to numerically solve the one-dimensional Shallow Water Equations. :
\begin{eqnarray} \dfrac{\partial h(x,t)}{\partial t} + \dfrac{\partial}{\partial x}(h(x,t).v(x,t)) & = & 0\\ \dfrac{\partial h(x,t).v(x,t)}{\partial t} + \dfrac{\partial}{\partial x}(h(x,t).v(x,t)^2 + \dfrac{1}{2}g.h(x,t)) & = & 0 \\ x \in [0, L], t \in \mathbb{R} \end{eqnarray}
where $h$ is the fluid height, $v$ its speed and $g$ the gravitational acceleration. The boundary conditions on the speed are $v(0,t) = v(L,t) = 0, \forall t$. Because of the lack of boundary condition on the height, I impose that the flux applied on the height is zero at the boundaries.
I used the Roe solver as described in "Riemann Solvers and Numerical Methods for Fluid Dynamics" by Eleuterio F. Toro (p. 357) and added an Harten entropy fix. My variables are $h(x,t)$ and $h(x,t).v(x,t)$.
To verify my code I implemented two method :
- I compute the total energy of the fluid E as : $$E(t) = \dfrac{1}{2} \int_0^L \left( h.v^2 + g.h^2\right)dx $$ and see if it remains constant.
- I compute the integrals of the equation on the whole "volume": $$\dfrac{\partial}{\partial t} \int_0^L h dx = 0$$ $$\dfrac{\partial}{\partial t} \int_0^L (h.v) dx + \dfrac{1}{2}g[h(L,t)^2 - h(0,t)^2] = 0$$
After running my code, I see that the second method seems to works as both integrals are around $10^{-13}$. However by the plot of the total energy through time show a loss of energy that I can't explain. I try to increase the number of cells, to decrease the step-time and to change the Harten entropy fix parameter, without success.
To compare my results I also used different codes online, such as James Adams' one on Mathworks (https://fr.mathworks.com/matlabcentral/fileexchange/46475-1d-shallow-water-equations-dam-break). But they all shown this energy loss.
Intuitively I would explain this loss by a numerical viscosity, but I'm not sure as CFD is not my domain of competence. So I come here to have both an explanation and a solution please.
The code I wrote is:
clear all;
close all;
clc;
warning off;
%% ------------------------------------------------------------ PREPARATION
%% Parametres :
N = 250; % Nombre de cellules
rho = 1000; % Densite (kg/m^3);
g = 9.81; % acceleration gravitationnelle (kg/m/s^2)
L = 1; % Longueur du reservoir (m)
h_bar = 1; % Niveau de fluide moyen
Dx = L/N; % Pas spatial (m)
t_fin = 5; % Temps final de la simulation (s)
xx = linspace(0,L,N); % Vecteur pour plot_fluid.m
x_vec = ((1:N)-1/2)*Dx; % Vecteur des les centres des cellules
C = 20; % Fréquence d'affichage des plots pour la vidéo
%% Switchs:
VIDEO = 0; % Sauvegarde une vidéo
%% Initialisation/Allocation :
H = NaN(N,1); % Hauteur du fluide
V = NaN(N,1); % Moment du fluide
D = 0; % Position du réservoir
dot_D = 0; % Vitesse du réservoir
U_new = NaN(N,2);
u = NaN(N,1); % Vitesse du fluide dans le ref. du réservoir
X_f = NaN(1);
lambda = NaN(N,2);
time = 0;
n = 1;
telapsed = 0;
compteur = 1;
newSubFolder = sprintf('./IMAGES');
if ~exist(newSubFolder, 'dir')
mkdir(newSubFolder);
else
rmdir(newSubFolder,'s');
mkdir(newSubFolder);
end
%% -------------------------------------------------- SIMULATION TEMPORELLE
fprintf('Début de la simulation temporelle\n')
while (time(end) < t_fin)
tstart = tic;
if(mod(n,200) == 0)
if(telapsed > 60)
fprintf('It. n°: %d - T. Simu: %.3f/%.3f s - T. Exec: %0.3f min\n',n,time(n),t_fin, telapsed/60)
else
fprintf('It. n°: %d - T. Simu: %.3f/%.3f s - T. Exec: %0.3f s\n',n,time(n),t_fin, telapsed)
end
end
if(n==1)
%% Conditions initiales
%H(:,n) = h_bar*ones(N,1); % Fluide à plat
H(:,n) = h_bar*(1+1*exp(-(x_vec-L/2).^2/(2*(0.1)^2))); % Hump
%H(:,n) = h_bar*(1+0.25*sin(x_vec*2*pi/L)); % Sin wave
V(:,n) = zeros(N,1); % Fluide au repos
else
%% Mise à jour des variables
H(:,n) = U_new(:,1);
V(2:N-1,n) = U_new(2:N-1,2);
end
% Application des CL
V(1,n) = dot_D;
V(N,n) = dot_D;
% Calcul du vecteur vitesse
u(:,n) = V(:,n) - dot_D;
%% Calcul du pas de temps
% Valeurs propres de la matrice de Roe
Dt = 0.99*Dx/max( abs(u(:,n)) + sqrt(g*H(:,n)) );
%% Calcul du vecteur [α1; α2] au temps suivant par un schéma Godunov-Roe
U = [H(:,n), V(:,n)];
U_new = Roe_Solver(U,dot_D,g,Dt,Dx,N);
%% Calcul de la position du centre de gravité du fluide:
mf(n) = rho*Dx*sum(H(:,n));
X_f(n) = rho*Dx*sum(x_vec'.*H(:,n))/mf(n);
%% Calcul de l'énergie totale du fluide
E_tot(n) = Dx*sum(0.5*H(:,n).*V(:,n).^2) + Dx*sum(0.5*g*H(:,n).^2);
%% Verification from Saint-Venant equations
if (n>1)
EQ1_Verif(n) = Dx*(sum(H(:,n))-sum(H(:,n-1)))/Dt;
EQ2_Verif(n) = Dx*(sum(V(:,n))-sum(V(:,n-1)))/Dt + g*H(N,n) - g*H(1,n);
end
%% Plots pour la vidéo:
if(VIDEO==1 && ( mod(n,C) == 0 || n == 1) )
fig = figure('visible','off');
set(fig, 'Position', [0 0 600 500])
titlefig = sprintf('Système ballotant (t=%.4f s)',time(n));
plot_fluid(time,D,u(:,n),H,X_f,L,h_bar,titlefig,n,N,xx);
filename = sprintf('IMG-%d.png',compteur);
saveas(fig,strcat(newSubFolder,'/',filename));
clf(fig);
compteur = compteur +1;
end
% Incrementation
n = n+1;
% Vecteur temporel
time(n) = time(end)+Dt;
% Temps écoulé
telapsed = toc(tstart) + telapsed;
end
time = time(1:end-1);
%% Calculs d'erreurs
figure(101)
plot(time,mf)
title('Fluid total mass');
figure(102)
hold on
plot(time,EQ1_Verif,'r');
plot(time,EQ2_Verif,'b')
legend('H','V');
figure(103)
hold on
plot(time,E_tot,'r');
title('Total Energy');
fprintf('Fin d''exécution\n');
and
function U_new = Roe_Solver(U,dot_D,g,Dt,Dx,N)
U_new = NaN(N,2);
for i=1
% Calcul de F_{i-1/2}
F_m = [0 ; 0];
% Calcul de F_{i+1/2}
UL = U(1,:);
UR = U(2,:);
F_p = Flux(UL,UR,dot_D,g);
U_new(1,:) = U(1,:) - (Dt/Dx)*(F_p - F_m)';
end
for i=N
% Calcul de F_{i-1/2}
UL = U(N-1,:);
UR = U(N,:);
F_m = Flux(UL,UR,dot_D,g);
% Calcul de F_{i+1/2}
F_p = [0 ; 0];
U_new(N,:) = U(N,:) - (Dt/Dx)*(F_p - F_m)';
end
for i=2:N-1
% Calcul de F_{i-1/2}
UL = U(i-1,:);
UR = U(i,:);
F_m = Flux(UL,UR,dot_D,g);
% Calcul de F_{i+1/2}
UL = U(i,:);
UR = U(i+1,:);
F_p = Flux(UL,UR,dot_D,g);
U_new(i,:) = U(i,:) - (Dt/Dx)*(F_p - F_m)';
end
end
function F_p = Flux(UL,UR,dot_D,g)
% cf. E. Toro, "Riemann Solvers and Numerical Methods for Fluid Dynamics", P. 357
HL = UL(1);
VL = UL(2);
HR = UR(1);
VR = UR(2);
% Roe average values
H_tilde = (HL + HR)/2;
V_tilde = (VL + VR)/2;
A_tilde = [V_tilde - dot_D, H_tilde; g, V_tilde - dot_D];
[EVec, EVal] = eig(A_tilde);
% Averaged eigenvalues
lambda_tilde_1 = EVal(1,1);
lambda_tilde_2 = EVal(2,2);
% Averaged right eigenvectors
K_tilde_1 = EVec(:,1);
K_tilde_2 = EVec(:,2);
% Wave strenghts
alpha_tilde_vec = inv(EVec)*[HR - HL; VR - VL];
alpha_tilde_1 = alpha_tilde_vec(1);
alpha_tilde_2 = alpha_tilde_vec(2);
% Flux
FR = [HR*(VR - dot_D); (VR^2)/2 - VR*dot_D + g*HR];
FL = [HL*(VL - dot_D); (VL^2)/2 - VL*dot_D + g*HL];
F_p = 0.5*(FL + FR) ...
- 0.5*Harten_fix(lambda_tilde_1)*alpha_tilde_1*K_tilde_1 ...
- 0.5*Harten_fix(lambda_tilde_2)*alpha_tilde_2*K_tilde_2;
end
function abs_lambda_Harten = Harten_fix(lambda)
delta = 0.5; % Harten entropy fix parameter
for(i=1:2)
if(lambda < -delta || lambda > delta)
abs_lambda_Harten = abs(lambda);
else
abs_lambda_Harten = 0.5*( (lambda^2)/delta + delta);
end
end
end
Thank you for your help.
