-
Notifications
You must be signed in to change notification settings - Fork 0
/
fluidDriver.m
139 lines (112 loc) · 2.59 KB
/
fluidDriver.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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
% A driver for just the fluid solver. See ibDriver.m for details...
%
% License: This code is free to use for any purposes, provided
% any publications resulting from the use of this code
% reference the original code/author.
%
% Author: Samuel Isaacson ([email protected])
% Date: 11/2007
%
% Please notify the author of any bugs, and contribute any
% modifications or bug fixes back to the original author.
%
% Disclaimer:
% This code is provided as is. The author takes no responsibility
% for its results or effects.
% params:
N = 64;
rho = 1;
mu = .01; %.0005;
h = 1 / N;
dt = .75 * h;
Nt = 500;
% Euclidean Mesh
X = (0:(N-1))' * h;
[X,Y] = meshgrid(X,X);
X = X';
Y = Y';
x = reshape(X, N*N, 1);
y = reshape(Y, N*N, 1);
% initial fluid flow and force:
u0 = zeros(N*N,2);
u0 = [(sin(2*pi*x) .* sin(2*pi*y)) (sin(2*pi*x) .* cos(2*pi*y)) ];
f = zeros(N*N,2);
% u0(30:N:(N*N-N+30),2) = .1;
% u0(31:N:(N*N-N+31),2) = .1;
% u0(32:N:(N*N-N+32),2) = .1;
% u0(33:N:(N*N-N+33),2) = .1;
% u0(257:512,1) = .1;
%----------END OF INITIALIZATION CODE-------------%
% declare fluid variables
uX = zeros(N*N,Nt+1);
uY = uX;
u = u0;
uX(:,1) = u0(:,1);
uY(:,1) = u0(:,2);
p = zeros(N*N,Nt+1);
% time step
for(i = 1:Nt)
i
[u pT] = fluidSolver(u, f, rho, mu, dt, h, N);
uX(:,i+1) = u(:,1);
uY(:,i+1) = u(:,2);
p(:,i+1) = pT;
end
% plotting stuff:
xmax = max(max(uX));
xmin = min(min(uX));
ymax = max(max(uY));
ymin = min(min(uY));
mmax = max([xmax ymax]);
mmin = min([xmin ymin]);
pmax = max(max(p));
pmin = min(min(p));
hh = figure('DoubleBuffer', 'on');
subplot(2,2,1);
pcolor(X - h/2, Y - h/2, reshape(uX(:,1),N,N));
caxis([mmin mmax]);
colorbar;
title(['u at t = 0']);
xlabel('x');
ylabel('y');
subplot(2,2,2);
pcolor(X - h/2, Y - h/2, reshape(uY(:,1),N,N));
caxis([mmin mmax]);
colorbar;
title(['v at t = 0']);
xlabel('x');
ylabel('y');
subplot(2,2,3);
pcolor(X - h/2, Y - h/2, reshape(p(:,1),N,N));
caxis([pmin pmax]);
colorbar;
title(['p at t = 0']);
xlabel('x');
ylabel('y');
drawnow
pause
for(i = 2:(Nt+1))
subplot(2,2,1);
pcolor(X - h/2, Y - h/2, reshape(uX(:,i),N,N));
caxis([mmin mmax]);
colorbar;
title(['u at t = ' num2str(dt*(i-1))]);
xlabel('x');
ylabel('y');
subplot(2,2,2);
pcolor(X - h/2, Y - h/2, reshape(uY(:,i),N,N));
caxis([mmin mmax]);
colorbar;
title(['v at t = ' num2str(dt*(i-1))]);
xlabel('x');
ylabel('y');
subplot(2,2,3);
pcolor(X - h/2, Y - h/2, reshape(p(:,i),N,N));
caxis([pmin pmax]);
colorbar;
title(['p at t = ' num2str(dt*(i-1))]);
xlabel('x');
ylabel('y');
drawnow;
pause(.1);
end