-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathall_ops_v2_INT.m
142 lines (128 loc) · 5.25 KB
/
all_ops_v2_INT.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
function [f_1,f_2,BUp_1, BUp_2, BLW_1, BLW_2]=all_ops_v2_INT(partition,G,omega,force,Cur_Iter)
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% Function for performing all local MCMP SC LBM operations before the
% streaming portion of the algorithm - version for multiscale modeling
% ---------------------------------------------------------------------
% Last modified: December 8th 2014
% ---------------------------------------------------------------------
% Input: partion (structure), omega (inverse relaxation time),
% force (bulk force, gravity or pressure drop),
% Cur_Iter for velocity computation
% ---------------------------------------------------------------------
% Output: density distributions for each of the two components, upper
% and lower density buffers (extra rows) to be send to corresponding
% neighboring partitions
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% --- RETRIEVE/PROCESS DATA FROM PARTITIONS
[Nr Mc N_c]=size(partition.f_1);
Nc=N_c;
f_1=partition.f_1;
f_2=partition.f_2;
Fxs_1=partition.Fxs_1;
Fxs_2=partition.Fxs_2;
Fys_1=partition.Fys_1;
Fys_2=partition.Fys_2;
forceINT=partition.forceINT;
omega_1=omega; omega_2=omega;
% Indexing
[iabw1 jabw1]=find(f_1(:,:,9)~=0);
lena=length(iabw1);
% Linear index
ija=(jabw1-1)*Nr+iabw1;
% Density
rho_1=sum(f_1,3);
rho_2=sum(f_2,3);
% Buffers
BUp_1=rho_1(Nr-1,:);
BUp_2=rho_2(Nr-1,:);
BLW_1=rho_1(2,:);
BLW_2=rho_2(2,:);
% --- LBM CONSTANTS
% Direction vectors
C_x=[1 0 -1 0 1 -1 -1 1 0];
C_y=[0 1 0 -1 1 1 -1 -1 0];
% Density weights
w0=16/36.;w1=4/36.;w2=1/36.;
W=[w1 w1 w1 w1 w2 w2 w2 w2 w0];
cs2=1/3;cs2x2=2*cs2;cs4x2=2*cs2.^2;
f1=3.; f2=4.5; f3=1.5;
NxM=Nr*Mc;
% --- IMMISCIBLE FLUIDS INTERACTIOS
psi_1=rho_1;
psi_2=rho_2;
% Fsum - interaction forces
[Fxtemp_1 Fytemp_1]=Fsum(Nr,Mc,iabw1,jabw1,lena,psi_1,partition.Up_1,partition.Low_1);
[Fxtemp_2 Fytemp_2]=Fsum(Nr,Mc,iabw1,jabw1,lena,psi_2,partition.Up_2,partition.Low_2);
% Interaction forces, final form
Fx_1=-G*psi_1.*Fxtemp_2; Fx_2=-G*psi_2.*Fxtemp_1;
Fy_1=-G*psi_1.*Fytemp_2; Fy_2=-G*psi_2.*Fytemp_1;
% --- MACROSCOPIC, COMPOSITE, AND EQUILIBRIUM VELOCITIES
% Velocity
ux_1=zeros(Nr,Mc);uy_1=zeros(Nr,Mc); ux_2=zeros(Nr,Mc);uy_2=zeros(Nr,Mc);
if Cur_Iter>0
ux_1=zeros(Nr,Mc);uy_1=zeros(Nr,Mc); ux_2=zeros(Nr,Mc);uy_2=zeros(Nr,Mc);
for ic=1:N_c-1
ux_1=ux_1+C_x(ic).*f_1(:,:,ic);
uy_1=uy_1+C_y(ic).*f_1(:,:,ic);
ux_2=ux_2+C_x(ic).*f_2(:,:,ic);
uy_2=uy_2+C_y(ic).*f_2(:,:,ic);
end
end
% Composite velocity
ux=(ux_1*omega_1+ux_2*omega_2)./(rho_1*omega_1+rho_2*omega_2);
uy=(uy_1*omega_1+uy_2*omega_2)./(rho_1*omega_1+rho_2*omega_2);
% Equilibrium velocities
for nzr=1:lena
if rho_1(ija(nzr))>0
ux_1(ija(nzr))=ux(ija(nzr))+Fx_1(ija(nzr))/rho_1(ija(nzr))/omega_1+Fxs_1(ija(nzr))/omega_1;
uy_1(ija(nzr))=uy(ija(nzr))+Fy_1(ija(nzr))/rho_1(ija(nzr))/omega_1+Fys_1(ija(nzr))/omega_1;
end
if rho_2(ija(nzr))>0
ux_2(ija(nzr))=ux(ija(nzr))+Fx_2(ija(nzr))/rho_2(ija(nzr))/omega_2+Fxs_2(ija(nzr))/omega_2;
uy_2(ija(nzr))=uy(ija(nzr))+Fy_2(ija(nzr))/rho_2(ija(nzr))/omega_2+Fys_2(ija(nzr))/omega_2;
end
end
% --- EQUILIBRIUM DENSITY DISTRIBUTIONS
feq_1=f_eqilibrium(ux_1,uy_1,rho_1,ija,NxM,Nr,Mc,N_c);
feq_2=f_eqilibrium(ux_2,uy_2,rho_2,ija,NxM,Nr,Mc,N_c);
% --- COLLISIONS
f_1=(1.-omega_1).*f_1+omega_1.*feq_1;
f_2=(1.-omega_2).*f_2+omega_2.*feq_2;
% --- INTERFACE TRACKING FOR MULTISCALE MODELING
% Volume and interface force addition
% Interface will be 3-4 nodes thick - simple tracking
% Change the index naming in "Include" loop at the end
% [iIn4 jIn4]=find((((rho_2-rho_1)<=0)&rho_1<=2&rho_1>6e-2));
% Complex tracking with filtering for diluted systems
% Here both fluids have nominal bulk densities of 2.0 and dissolved
% densities of 6e-2
rho_0=1.95;
% Considers bulk of the droplet
F1=(rho_1>=rho_0);
% Considers interface for non-diluted systems
% (in diluted it consumes part of the bulk)
F2=((((rho_2-rho_1)<=0)&rho_1<=2&rho_1>6e-2));
% 2-3 nodes thick interface, adjust thickness by adjusting
% rho_0; -1, 1 and 0 entries
FN=F2-F1;
% For visualization
% spy(FN>0)
% Get the nodes - interface entries are 1
[iInt3,jInt3]=find(FN>0);
% --- INCLUDE INTERFATIAL AND VOLUME/BODY FORCE
for ic=1:N_c;
for ia=1:lena
i=iabw1(ia); j=jabw1(ia);
% Interface
for jin=1:length(iInt3)
if (i==iInt3(jin)&j==jInt3(jin))
f_1(i,j,ic)=f_1(i,j,ic)+forceINT(i,j,ic);
f_2(i,j,ic)=f_2(i,j,ic)+forceINT(i,j,ic);
end
end
% Volume
f_1(i,j,ic)=f_1(i,j,ic)+force(ic);
f_2(i,j,ic)=f_2(i,j,ic)+force(ic);
end
end
end