1
$\begingroup$

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 :

  1. 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.
  2. 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.

|cite|improve this question
$\endgroup$
1
$\begingroup$

Energy will not be conserved. You guessed it right that it is due to numerical viscosity. For smooth solutions, energy should be conserved to a better accuracy. You can use a high order method to get better energy conservation but it will never be perfectly conserved. The energy actually plays the role of entropy for this model.

You can use a central difference scheme with say RK4 time integration which would be stable for smooth solutions and will have better energy conservation.

You will find more literature in numerical weather prediction papers regarding energy conservation. There are staggered grid schemes that can achieve this, see [1].

There are some FE methods that also conserve energy, see [2].

[1] Arakawa and Lamb, Monthly Weather Review, Vol. 109, 1980. https://doi.org/10.1175/1520-0493(1981)109%3C0018:APEAEC%3E2.0.CO;2

[2] Mcrae and Cotter, QJRMS, 2013 https://doi.org/10.1002/qj.2291

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thanks for your answer and the references Praveen Chandrashekar ! So it's simply that even with an entropy fix the energy is still not conserved. I will try to change the scheme and see if the energy is "less lost". $\endgroup$ – Anthony Jul 26 '18 at 12:02

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.