Deprecated: Assigning the return value of new by reference is deprecated in /customers/1/f/c/dahoiv.net/httpd.www/programmering/wp-includes/cache.php on line 36 Deprecated: Assigning the return value of new by reference is deprecated in /customers/1/f/c/dahoiv.net/httpd.www/programmering/wp-includes/query.php on line 21 Deprecated: Assigning the return value of new by reference is deprecated in /customers/1/f/c/dahoiv.net/httpd.www/programmering/wp-includes/theme.php on line 507 Numerical solution of the stationary advection diffusion equation. | Programming
« Bash script for SFV and uncompressing of .r00 files
My .emacs file »

Numerical solution of the stationary advection diffusion equation.

June 23rd, 2008 by Daniel Høyer Iversen

The stationary advection diffusion equation solved with the finite difference method.


newSuperDriver.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

function U = newSuperDriver2(M, N)
p=params(2*M,N);
[A,B]=matrise();
v=A(B); %Løser ligningsystemet med 2M skritt i x-retning og N skritt i y-retning
%v=solver(A,B);
 
p=params(M,N);
[A,B]=matrise();
u=A(B); %Løser ligningsystemet med M skritt i x-retning og N skritt i y-retning
%u=solver(A,B);
 
U=feval(p.u_star,zeros(N+1,1) ,[0;p.y]);
for (i =0:M-1)
     U=[U; feval(p.u_star,(i+1)*p.h,0)];
     for (j =1:N)
          U=[U; (4*v(j+i*2*N+N)-u(j+i*N))/3]; %Ekstrapolerer mellom punktene fra de to løsningene
      end
end

params.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

function p = params(m,n)
 
persistent S;
if (m ~= 0),
	S.N = n; %Antall skritt i y-retning
	S.M = m; %Antall skritt i x-retning
	S.ny =-0.009; %Verdien til ny
 
	S.h = 1/S.M; %Beregner skrittlengden i x-retning
	S.k = 1/S.N; %Beregner skrittlengden i y-retning
	S.NM = S.M*S.N;
 
	S.vektorfelt=@vektorfelt_praktisk;  %Vektorfeltet
	S.HS_func=@HS_func_praktisk2;        %f(x,y)
        S.u_star=@u_star_praktisk2;          %randbetingelsene
 
        S.x=ones(S.N,1);
        S.y=[S.k:S.k:1]';
 
        x=[S.h:S.h:1];
        t=x(ones(1,S.N),:);
        S.noder=[t(:) repmat(S.y,S.M,1)]; %Inneholder koordinatene til alle nodene
end
 
p = S;

matrise.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

function [A,B] = matrise()
 
p = params(0,0);
 
V=feval(p.vektorfelt, p.noder(:,1), p.noder(:,2)); %Beregner vektorfeltet
 
B=spdiags([1],[1],2,1);
P = ones(p.NM,1)*(-2*(1/p.h^2+1/p.k^2)); %Verdiene til diagonalen
 
konst1=1/p.h^2;
konst2=1/(2*p.h*p.ny)*V(:,1);
W=konst1-konst2;                         %Verdiene til vest-konstantene
E=konst1+konst2;                         %Verdiene til øst-konstantene
 
konst1=1/p.k^2;
konst2=1/(2*p.k*p.ny)*V(:,2);
N=konst1+konst2;                         %Verdiene til nord-konstantene
S=konst1-konst2;                         %Verdiene til sør-konstantete
 
B(1:p.N)=W(1:p.N).*feval(p.u_star,zeros(p.N,1),p.y); %første del av B-vektoren
B(1)=B(1)+S(1)*feval(p.u_star,p.h,0);    
 
%Retter opp sør/nord/B-konstantene pga randkravene
for i=1:p.M-1
	S(i*p.N)=S(i*p.N)+N(i*p.N);
	N(i*p.N)=0;
	B(i*p.N+1)=S(i*p.N+1)*feval(p.u_star,(i+1)*p.h,0);
	S(i*p.N+1)=0;
end
S(end)=S(end)+N(end);
 
%Retter opp W-konstantene pga randkravene
W(end-p.N+1:end)= W(end-p.N+1:end)+E(end-p.N+1:end);
 
W(1:p.N)=[];
W=[W;zeros(p.N,1)];
 
E(end-p.N+1:end)=[];
E=[zeros(p.N,1);E];
 
N(end)=[];
N=[0;N];
 
S(1)=[];
S=[S;0];
 
A = spdiags([W, S, P, N, E], [-p.N -1 0 1 p.N], p.NM,p.NM); %Setter de fem diagonalene inn i matrise A
B(p.NM)=0;
B= feval(p.HS_func, p.noder(:,1), p.noder(:,2) )/p.ny-B; %Beregner F-vektor og trekker fra B

HS_func_praktisk.m

1
2
3
4
5
6
7
8
9
10
11
12
13

function f = HS_func_praktisk(x,y)
 
p=params(0,0);
 
f = 0*x;
 
%f(round(p.NM/2+p.N/4))=10000;
 
for i=round(p.NM/2+p.N/2):round(p.NM/2+p.N/2+p.N*0.05)
     for j=0:round(p.M*0.05)
     f(i+p.N*j)=2000;
     end
end

u_star_praktisk.m

1
2
3
4
5
6
7
8
9
10
11
12
13

function u = u_star_praktisk(x,y)
 
p=params(0,0);
 
%u=5+0.2*log(50*y+exp(1));
%u=1+0.4*sin(pi*y-pi*0.5);
%u=2;
 
%u=y*4;
%u=cos(x*s)+4;
 
u=exp(-((y-0.65)/0.1).^2)*2.6;
%u=4;

vektorfelt_praktisk.m

1
2
3
4
5
6
7
8
9
10

 
function v = vektorfelt_praktisk(x,y)
p=params(0,0);
 
%v = 5*[(4*y-4).*sin(2*atan((y-1)./(x-1)))+2*(2*(y-1).^2 + 2*(2*(x-1).^2.*cos(2*atan((y-1)./(x-1)))./((x-1).*(1+((y-1).^2./(x-1).^2)))) ), -(4*x-4).*sin(2*atan((y-1)./(x-1)))+2*(2*(y-1).^2 + 2*(2*(x-1).^2.*cos(2*atan((y-1)./(x-1)).*(y-1))./((x-1).^2.*(1+((y-1).^2./(x-1).^2))))  )];
 
v= 5*[(2*y-2).*sin(2*atan((y-1)./(x-1))) + 2*(y.^2-2*y+2+x.^2-2*x).*cos(2*atan((y-1)./(x-1)))./((x-1).*(1+(y-1).^2./(x-1).^2)) ,   -(2*x-2).*sin(2*atan((y-1)./(x-1))) +  2*(y.^2-2*y+2+x.^2-2*x).*cos(2*atan((y-1)./(x-1)))./((x-1).^2.*(1+(y-1).^2./(x-1).^2)).*(y-1)];
 
v(p.NM-p.N:p.NM,1)=0;
v(p.NM-p.N:p.NM,2)=0;
1

 
1

 

Result of the simulation:
Result of the stationary advection diffusion equation.

Do you want to use this code?

3 Responses to “Numerical solution of the stationary advection diffusion equation.”

  1. bob&alice

    koden er til matlab, er den ikke?

    Nov 24th, 2008    4:54 pm   
  2. bob&alice

    duh, glem det, sorry. glemte å se i høyre marg.. :)

    Nov 24th, 2008    4:55 pm   
  3. fathi

    Iwould like to use this code, please send it to me .
    Best regard
    Fathi

    Feb 5th, 2013    11:51 am   

Leave a Response

« Bash script for SFV and uncompressing of .r00 files
My .emacs file »