446
11 Probl` emes aux limites en dimension un
Les Programmes 83 et 84 calculent respectivement la viscosit´ e artificielle pour
le sch´ ema d´ ecentr´ e et le sch´ ema de fitting exponentiel. Ces m´ ethodes sont
choisies par l’utilisateur en prenant un param` etre stabfun (Programme 82)
´ egal ` a @upwind ou @sgvisc. La fonction sgvisc utilise la fonction bern (Programme 85) pour ´ evaluer l’inverse de la fonction de Bernoulli.
Programme 83 - upwind : Viscosit´ e artificielle du sch´ ema d´ ecentr´ e
function [Kupw,rhsbc] = upwind(Nx,h,nu,beta)
%UPWIND viscosit´ e artificielle du sch´ ema d´ ecentr´ e
% Matrice de raideur et second membre.
Peclet=0.5*h*abs(beta);
dd=(Peclet(1:Nx-1)+Peclet(2:Nx))/h;
ld=-Peclet(2:Nx-1)/h;
ud=-Peclet(1:Nx-2)/h;
Kupw=spdiags([[ld 0]’,dd’,[0 ud]’],-1:1,Nx-1,Nx-1);
rhsbc = - [Peclet(1)/h, Peclet(Nx)/h];
return
Programme 84 - sgvisc : Viscosit´ e artificielle du sch´ ema de fitting exponentiel
function [Ksg,rhsbc] = sgvisc(Nx, h, nu, beta)
%SGVISC viscosit´ e artificielle du sch´ ema de fitting exponentiel
% Matrice de raideur et second membre.
Peclet=0.5*h*abs(beta)./nu;
[bp, bn]=bern(2*Peclet);
Peclet=Peclet-1+bp;
dd=(nu(1:Nx-1).*Peclet(1:Nx-1)+nu(2:Nx).*Peclet(2:Nx))/h;
ld=-nu(2:Nx-1).*Peclet(2:Nx-1)/h;
ud=-nu(1:Nx-2).*Peclet(1:Nx-2)/h;
Ksg=spdiags([[ld 0]’,dd’,[0 ud]’],-1:1,Nx-1,Nx-1);
rhsbc = - [nu(1)*Peclet(1)/h, nu(Nx)*Peclet(Nx)/h];
return
Programme 85 - bern : Evaluation de la fonction de Bernoulli
function [bp,bn]=bern(x)
%BERN Evaluation de la fonction de Bernoulli
xlim=1e-2; ax=abs(x);
if ax==0,
bp=1; bn=1; return
end
if ax>80
if x>0, bp=0.; bn=x; return
else, bp=-x; bn=0; return, end
end