-
Notifications
You must be signed in to change notification settings - Fork 21
/
Copy pathyu.m
161 lines (126 loc) · 5.3 KB
/
yu.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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
for i=1:nx
for j=1:ny
if isfluid(i,j) == 0
%% b=2
if isfluid(i,j+1) == 1
b=2;
a=4;
y= jj + sqrt( R^2 - (i-ii)^2 ) ;
delta(i,j,a)= abs(j+1 - y);
%%%%%%%%%%%collision%%%%%%%%%%%%%
UW(i,j)= R * wall_rotation;
VW(i,j)=0;
fwall= f(i,j+2,a) + delta(i,j,a)* ( f(i,j+1,a) - f(i,j+2,a) ) + 2 * w(a) * Rho(i,j+1) * 3 * ( ex(a)* UW(i,j) + ey(a) * VW(i,j));
f(i,j,b)= fwall + delta(i,j,a)/( 1 + delta(i,j,a) ) * ( f(i,j+1,b) - fwall );
%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% 4
if isfluid(i,j-1) == 1
b=4;
a=2;
y= jj - sqrt( R^2 - (i-ii)^2 ) ;
delta(i,j,a)= abs(y - (j-1));
%%%%%%%%%%%collision%%%%%%%%%%%%%
UW(i,j)= - R * wall_rotation;
VW(i,j)=0;
fwall= f(i,j-2,a) + delta(i,j,a)* ( f(i,j-1,a) - f(i,j-2,a) )+ 2 * w(a) * Rho(i,j-1) * 3 * ( ex(a)* UW(i,j) + ey(a) * VW(i,j)) ;
f(i,j,b)= fwall + delta(i,j,a)/( 1 + delta(i,j,a) ) * ( f(i,j-1,b) - fwall );
%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% 1
if isfluid(i+1,j) == 1
a=3;
b=1;
x= ii + sqrt( R^2 - (j-jj)^2 ) ;
delta(i,j,a)= abs(i+1 - x);
%%%%%%%%%%%collision%%%%%%%%%%%%%
UW(i,j)= 0;
VW(i,j)= - R * wall_rotation;
fwall= f(i+2,j,a) + delta(i,j,a)* ( f(i+1,j,a) - f(i+2,j,a) ) + 2 * w(a) * Rho(i+1,j) * 3 * ( ex(a)* UW(i,j) + ey(a) * VW(i,j));
f(i,j,b)= fwall + delta(i,j,a)/( 1 + delta(i,j,a) ) * ( f(i+1,j,b) - fwall );
%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% 3
if isfluid(i-1,j) == 1
a=1;
b=3;
x= ii - sqrt( R^2 - (j-jj)^2 ) ;
delta(i,j,a)= abs(x - (i-1));
%%%%%%%%%%%collision%%%%%%%%%%%%%
UW(i,j)= 0;
VW(i,j)= R * wall_rotation;
fwall= f(i-2,j,a) + delta(i,j,a)* ( f(i-1,j,a) - f(i-2,j,a) )+ 2 * w(a) * Rho(i-1,j) * 3 * ( ex(a)* UW(i,j) + ey(a) * VW(i,j)) ;
f(i,j,b)= fwall + delta(i,j,a)/( 1 + delta(i,j,a) ) * ( f(i-1,j,b) - fwall );
%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% 5
if isfluid(i+1,j+1) == 1
b = j-ii-jj-i;
c = 1/2*( ii^2 + (j-jj-i)^2 - R^2 );
x = 1/2*( -b + sqrt( b^2 - 4*c ) );
y = j + (x-i);
a=7;
bb=5;
delta(i,j,a) = sqrt(( (i+1 - x)^2 + (j+1 - y)^2) )/sqrt(2);
%%%%%%%%%%%collision%%%%%%%%%%%%%
teta= atan( abs( (y - jj) / (x - ii ) )); % angel of the intersection of wall point
UW(i,j)= abs( R * wall_rotation * cos(teta));
VW(i,j)= -abs(R * wall_rotation * sin(teta));
fwall= f(i+2,j+2,a) + delta(i,j,a)* ( f(i+1,j+1,a) - f(i+2,j+2,a) ) + 2 * w(a) * Rho(i+1,j+1) * 3 * ( ex(a)* UW(i,j) + ey(a) * VW(i,j));
f(i,j,bb)= fwall + delta(i,j,a)/( 1 + delta(i,j,a) ) * ( f(i+1,j+1,bb) - fwall );
%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% 7
if isfluid(i-1,j-1) == 1
b = j-ii-jj-i;
c = 1/2*( ii^2 + (j-jj-i)^2 - R^2 );
x = 1/2*( -b - sqrt( b^2 - 4*c ) );
y = j + (x-i);
a=5;
bb=7;
delta(i,j,a) = sqrt(1/2*( (x - (i-1) )^2 + (y - (j-1) )^2) );
%%%%%%%%%%%collision%%%%%%%%%%%%%
teta= atan( abs( (y - jj) / (x - ii ) ));
UW(i,j)= -abs(R * wall_rotation * cos(teta));
VW(i,j)= abs(R * wall_rotation * sin(teta));
fwall= f(i-2,j-2,a) + delta(i,j,a)* ( f(i-1,j-1,a) - f(i-2,j-2,a) )+ 2 * w(a) * Rho(i-1,j-1) * 3 * ( ex(a)* UW(i,j) + ey(a) * VW(i,j)) ;
f(i,j,bb)= fwall + delta(i,j,a)/( 1 + delta(i,j,a) ) * ( f(i-1,j-1,bb) - fwall );
%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% 6
if isfluid(i-1,j+1) == 1
b = -j-ii+jj-i;
c = 1/2*(ii^2 + (j-jj+i)^2 - R^2 );
x = 1/2*( -b - sqrt( b^2 - 4*c ) );
y = j - (x-i);
a=8;
bb=6;
delta(i,j,a) = sqrt(1/2*( (x - (i-1) )^2 + (j+1 -y)^2) );
%%%%%%%%%%%collision%%%%%%%%%%%%%
teta= atan( abs( (y - jj) / (x - ii ) ));
UW(i,j)= abs(R * wall_rotation * cos(teta));
VW(i,j)= abs(R * wall_rotation * sin(teta));
fwall= f(i-2,j+2,a) + delta(i,j,a)* ( f(i-1,j+1,a) - f(i-2,j+2,a) ) + 2 * w(a) * Rho(i-1,j+1) * 3 * ( ex(a)* UW(i,j) + ey(a) * VW(i,j));
f(i,j,bb)= fwall + delta(i,j,a)/( 1 + delta(i,j,a) ) * ( f(i-1,j+1,bb) - fwall );
%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% 8
if isfluid(i+1,j-1) == 1
b = -j-ii+jj-i;
c = 1/2*( ii^2 + (j-jj+i)^2 - R^2 );
x = 1/2*( -b + sqrt( b^2 - 4*c ) );
y = j - (x-i);
a=6;
bb=8;
delta(i,j,a) = sqrt(1/2*( (i+1 - x)^2 + (y - (j-1))^2) );
%%%%%%%%%%%collision%%%%%%%%%%%%%
teta= atan( abs( (y - jj) / (x - ii ) ));
UW(i,j)= -abs(R * wall_rotation * cos(teta));
VW(i,j)= -abs(R * wall_rotation * sin(teta));
fwall= f(i+2,j-2,a) + delta(i,j,a)* ( f(i+1,j-1,a) - f(i+2,j-2,a) ) + 2 * w(a) * Rho(i+1,j-1) * 3 * ( ex(a)* UW(i,j) + ey(a) * VW(i,j));
f(i,j,bb)= fwall + delta(i,j,a)/( 1 + delta(i,j,a) ) * ( f(i+1,j-1,bb) - fwall );
%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
end
end