« back to all changes in this revision
Viewing changes to ntable/src/ntable.hpp
- browse files at revision 1184
- compare with another revision
- download diff
- download tarball
- view history from revision 1184
- Committer: Qi Shao
- Date: 2017-08-22 01:18:11 UTC
- mfrom: (1179.1.4 ESyS-Darcy_clean)
- Revision ID: q.shao@uq.edu.au-20170822011811-aq63qj2swi99l7w4
Added fluid phase and its forces on solid particles. Flows of fluid are calculated according to the Darcy's Law. The coupling of fluid and DEM is based on the poro-elastic theory.
- files added:
- Doc/ESyS-Darcy
- Doc/ESyS-Darcy/examples
- Doc/ESyS-Darcy/examples/seepage
- files modified:
- Fields/FieldMaster.cpp
added
removed
ntable/src/ntable.hpp
21
21
m_array(),
22
22
m_idParticleMap(),
23
23
m_p0_global(),
24
m_pmax_global(), //fluid contents
24
25
m_dim(0),
25
26
m_alpha(0),
26
27
m_global_idx(0),
29
30
m_xsize(0),
30
31
m_ysize(0),
31
32
m_zsize(0),
32
m_valid(false)
33
m_valid(false),
34
m_fluidexist(false) //fluid contents
33
35
{
34
36
}
35
37
42
44
\param range grid spacing
43
45
\param alpha pair search cutoff
44
46
\param p0_global minimal corner (origin) of the global search space
45
\param ix x-index of the local origin
47
\param ix x-index of the local origin
46
48
\param iy y-index of the local origin
47
49
\param iz z-index of the local origin
48
50
*/
54
56
double range,
55
57
double alpha,
56
58
const Vec3& p0_global,
59
const Vec3& pmax_global, //fluid contents
57
60
int ix,
58
61
int iy,
59
62
int iz
60
63
)
61
64
: m_p0_global(p0_global),
65
m_pmax_global(pmax_global), //fluid contents
62
66
m_dim(range),
63
67
m_alpha(alpha),
64
68
m_xsize(x),
65
69
m_ysize(y),
66
70
m_zsize(z),
67
m_valid(true)
71
m_valid(true),
72
m_fluidexist(false) //fluid contents
68
73
{
69
74
m_global_idx=ix;
70
75
m_global_idy=iy;
71
m_global_idz=iz;
76
m_global_idz=iz;
72
77
m_min_corner=Vec3(m_p0_global.X()+m_dim*double(m_global_idx),
73
78
m_p0_global.Y()+m_dim*double(m_global_idy),
74
79
m_p0_global.Z()+m_dim*double(m_global_idz));
98
103
}
99
104
}
100
105
106
107
//fluid contents
108
/*!
109
clean up the cell array
110
*/
111
template<typename T>
112
void NeighborTable<T>::clear_cell_array()
113
{
114
for(unsigned int iter=0;iter<m_cellarray.size();iter++){
115
m_cellarray[iter].erase(m_cellarray[iter].begin(),m_cellarray[iter].end());
116
}
117
}
118
119
101
120
/*!
102
121
helper function which returns the index in the search array
103
into which a given position is mapped. Returns -1 if pos is
122
into which a given position is mapped. Returns -1 if pos is
104
123
outside the search space
105
124
106
125
\param the position
120
139
return idx;
121
140
}
122
141
142
143
//fluid contents
144
/*!
145
helper function which returns the index in the cell array
146
into which a given position is mapped. Returns -1 if pos is
147
outside the search space
148
149
\param the position
150
*/
151
template<typename T>
152
int NeighborTable<T>::cellindex(const Vec3& pos)
153
{
154
int idx=-1;
155
156
int ix=int(floor((pos.X()-m_min_corner.X()-m_dim)/m_xside))+1;
157
int iy=int(floor((pos.Y()-m_min_corner.Y()-m_dim)/m_yside))+1;
158
int iz=int(floor((pos.Z()-m_min_corner.Z()-m_dim)/m_zside))+1;
159
if((ix>=1)&&(ix<m_xcell-1)&&
160
(iy>=1)&&(iy<m_ycell-1)&&
161
(iz>=1)&&(iz<m_zcell-1)){
162
idx=m_ycell*m_zcell*ix+m_zcell*iy+iz;
163
}
164
return idx;
165
}
166
167
123
168
/*!
124
169
helper function which returns the index in the search array
125
170
from given x-,y- and z-indices.
135
180
return m_ysize*m_zsize*x+m_zsize*y+z;
136
181
}
137
182
138
/*!
139
check if a position is in the inner part
183
184
//fluid contents
185
/*!
186
helper function which returns the index in the cell array
187
from given x-,y- and z-indices.
188
189
\param x the x-index
190
\param y the y-index
191
\param z the z-index
192
\warning no checks
193
*/
194
template<typename T>
195
int NeighborTable<T>::cellindex(int x,int y,int z) const
196
{
197
return m_ycell*m_zcell*x+m_zcell*y+z;
198
}
199
200
201
/*!
202
check if a position is in the inner part
140
203
141
204
\param pos the position
142
205
*/
155
218
res=((ix>0) && (ix<m_xsize-1) && (iy>0) && (iy<m_ysize-1)&& (iz>0) && (iz<m_zsize-1));
156
219
}
157
220
158
return res;
221
return res;
159
222
}
160
223
161
224
175
238
// m_idParticleMap.insert(IdParticleMap::value_type(iter->getID(),&(*iter)));
176
239
m_idParticleMap.insert(make_pair(iter->getID(),&(*iter)));
177
240
m_array[idx].push_back(iter);
241
242
//fluid contents:
243
if(m_fluidexist){
244
int cellidx=cellindex(iter->getPos());
245
if(cellidx!=-1){
246
m_cellarray[cellidx].push_back(iter);
247
};
248
};
249
178
250
} else { // outside -> delete
179
251
typename list<T>::iterator h=iter;
180
252
iter--;
183
255
}
184
256
185
257
/*!
186
Build or rebuild the search array, calls clean_search_array.
258
Build or rebuild the search array, calls clean_search_array.
187
259
Particles outside the search space are removed from the list
188
*/
260
*/
189
261
template<typename T>
190
262
void NeighborTable<T>::build()
191
263
{
192
264
clear_search_array();
265
if(m_fluidexist){clear_cell_array();}; //fluid contents
193
266
// clean out id map
194
267
m_idParticleMap.clear();
195
268
for(typename list<T>::iterator iter=m_list.begin();
200
273
if(idx!=-1){
201
274
m_array[idx].push_back(iter);
202
275
m_idParticleMap[id]=&(*iter);
276
//fluid contents
277
if(m_fluidexist){
278
int cellidx=cellindex(iter->getPos());
279
if(cellidx!=-1){m_cellarray[cellidx].push_back(iter);};
280
}
203
281
} else { // outside -> delete
204
282
typename list<T>::iterator h=iter;
205
283
iter--;
210
288
}
211
289
212
290
/*!
213
Return pointer to particle at index.
214
*/
291
Return pointer to particle at index.
292
*/
215
293
template<typename T>
216
294
T* NeighborTable<T>::ptr(NeighborTable<T>::indextype idx)
217
295
{
229
307
230
308
/*!
231
309
Return pointer to particle with given id. Return NULL if the
232
table doesn't contain the particle.
310
table doesn't contain the particle.
233
311
234
312
\param id the id of the particle
235
*/
313
*/
236
314
template<typename T>
237
315
T* NeighborTable<T>::ptr_by_id(int id)
238
316
{
258
336
\param pos position
259
337
*/
260
338
template<typename T>
261
T *NeighborTable<T>::getNearestPtr(const Vec3& pos)
339
T *NeighborTable<T>::getNearestPtr(const Vec3& pos)
262
340
{
263
341
T* res=NULL;
264
342
265
343
// get grid index for pos
266
344
int idx=index(pos);
267
345
268
double dist=3.0*(m_dim*m_dim); // squared max. dist (grid diagonal)
346
double dist=3.0*(m_dim*m_dim); // squared max. dist (grid diagonal)
269
347
if(idx!=-1){
270
348
for(typename pointtype::iterator iter=m_array[idx].begin();
271
349
iter!=m_array[idx].end();
352
430
int xmin=int(floor((vmin.X()-m_p0_global.X())/m_dim))-m_global_idx;
353
431
int ymin=int(floor((vmin.Y()-m_p0_global.Y())/m_dim))-m_global_idy;
354
432
int zmin=int(floor((vmin.Z()-m_p0_global.Z())/m_dim))-m_global_idz;
355
// check for minimum outside
433
// check for minimum outside
356
434
xmin=(xmin < 0) ? 0 : xmin;
357
435
ymin=(ymin < 0) ? 0 : ymin;
358
436
zmin=(zmin < 0) ? 0 : zmin;
359
437
360
438
// maximum corner index
361
439
int xmax=int(floor((vmax.X()-m_p0_global.X())/m_dim))-m_global_idx;
362
440
int ymax=int(floor((vmax.Y()-m_p0_global.Y())/m_dim))-m_global_idy;
363
441
int zmax=int(floor((vmax.Z()-m_p0_global.Z())/m_dim))-m_global_idz;
364
// check for maximum outside
442
// check for maximum outside
365
443
xmax=(xmax < m_xsize) ? xmax : m_xsize-1;
366
444
ymax=(ymax < m_ysize) ? ymax : m_ysize-1;
367
445
zmax=(zmax < m_zsize) ? zmax : m_zsize-1;
368
446
369
447
return NTBlock<T>(this,xmin,xmax,ymin,ymax,zmin,zmax);
370
448
}
371
449
373
451
374
452
/*!
375
453
Return the inner region, i.e. excluding the boundary, as a NTBlock
376
*/
454
*/
377
455
template<typename T>
378
456
NTBlock<T> NeighborTable<T>::inner()
379
457
{
391
469
*/
392
470
template<typename T>
393
471
void NeighborTable<T>::addPairsToList(T_Handle<pairlist> list,int idx1,int idx2)
394
{
395
472
{
473
396
474
for(typename pointtype::iterator iter=m_array[idx1].begin();
397
475
iter!=m_array[idx1].end();
398
476
iter++){
406
484
list->push_back(make_pair(&(**iter),&(**iter2)));
407
485
} else {
408
486
list->push_back(make_pair(&(**iter2),&(**iter)));
409
}
487
}
410
488
}
411
489
}
412
490
}
435
513
list->push_back(make_pair(&(**iter),&(**iter2)));
436
514
} else {
437
515
list->push_back(make_pair(&(**iter2),&(**iter)));
438
}
516
}
439
517
}
440
518
}
441
519
}
443
521
}
444
522
445
523
/*!
446
Add pairs containing at least one flagged particle at two
447
given gridpoints to list.
524
Add pairs containing at least one flagged particle at two
525
given gridpoints to list.
448
526
449
527
\param list the pairlist
450
528
\param idx1 the index of the first gridpoint
452
530
*/
453
531
template<typename T>
454
532
void NeighborTable<T>::addPairsToListFlagged(T_Handle<pairlist> list,int idx1,int idx2)
455
{
456
533
{
534
457
535
for(typename pointtype::iterator iter=m_array[idx1].begin();
458
536
iter!=m_array[idx1].end();
459
537
iter++){
466
544
list->push_back(make_pair(&(**iter),&(**iter2)));
467
545
} else {
468
546
list->push_back(make_pair(&(**iter2),&(**iter)));
469
}
547
}
470
548
}
471
549
}
472
550
}
474
552
}
475
553
476
554
/*!
477
Add pairs containing at least one flagged particle at one
555
Add pairs containing at least one flagged particle at one
478
556
given gridpoint to list.
479
557
480
558
\param list the pairlist
496
574
list->push_back(make_pair(&(**iter),&(**iter2)));
497
575
} else {
498
576
list->push_back(make_pair(&(**iter2),&(**iter)));
499
}
577
}
500
578
}
501
579
}
502
580
}
543
621
}
544
622
545
623
/*!
546
Create a list of all pairs of neighboring particles involving a "flagged" particles
624
Create a list of all pairs of neighboring particles involving a "flagged" particles
547
625
and return handle to it.
548
626
*/
549
627
template<typename T>
550
628
T_Handle<typename NeighborTable<T>::pairlist> NeighborTable<T>::getNewList()
551
629
{
552
T_Handle<typename NeighborTable<T>::pairlist> nlist=new typename NeighborTable<T>::pairlist;
630
T_Handle<typename NeighborTable<T>::pairlist> nlist=new typename NeighborTable<T>::pairlist;
553
631
554
632
for(int ix=0;ix<m_xsize;ix++){
555
633
for(int iy=0;iy<m_ysize;iy++){
586
664
}
587
665
588
666
/*!
589
Get list of all particles along a given plane. Naive implementation, i.e. check all particles
590
for distance to plane. The plane is given by one point an the normal.
667
Get list of all particles along a given plane. Naive implementation, i.e. check all particles
668
for distance to plane. The plane is given by one point an the normal.
591
669
592
670
\param orig The origin of the plane
593
671
\param normal The normal of the plane
612
690
}
613
691
614
692
/*!
615
Get list of all particles near a given sphere.
693
Get list of all particles near a given sphere.
616
694
Naive implementation, i.e. check all particles for distance to sphere.
617
695
618
696
\param centre The centre of the sphere
647
725
T_Handle<typename NeighborTable<T>::particlelist> NeighborTable<T>::getParticlesNearTriangle(const Triangle& Tr)
648
726
{
649
727
T_Handle<typename NeighborTable<T>::particlelist> nlist=new typename NeighborTable<T>::particlelist;
650
728
651
729
// work out search area from bounding box of the triangle
652
730
Vec3 v_min=Tr.getBoundingBoxMin()-Vec3(m_dim,m_dim,m_dim);
653
731
Vec3 v_max=Tr.getBoundingBoxMax()+Vec3(m_dim,m_dim,m_dim);
654
732
// check if completely outside
655
733
if((v_min.X()<m_max_corner.X())&&(v_min.Y()<m_max_corner.Y())&&(v_min.Z()<m_max_corner.Z())&&
656
734
(v_max.X()>m_min_corner.X())&&(v_max.Y()>m_min_corner.Y())&&(v_max.Z()>m_min_corner.Z())){
657
735
658
736
NTBlock<T> TriangleBlock=block(v_min,v_max);
659
737
for(typename NTBlock<T>::iterator iter=TriangleBlock.begin();
660
738
iter!=TriangleBlock.end();
664
742
}
665
743
}
666
744
}
667
745
668
746
return nlist;
669
747
}
670
748
673
751
of the bounding box of the triangle; (could be improved - search area via Bresenham)
674
752
675
753
\param E the Edge
676
*/
754
*/
677
755
template<typename T>
678
756
T_Handle<typename NeighborTable<T>::particlelist> NeighborTable<T>::getParticlesNearEdge(const AEdge* E)
679
757
{
685
763
// check if completely outside
686
764
if((v_min.X()<m_max_corner.X())&&(v_min.Y()<m_max_corner.Y())&&(v_min.Z()<m_max_corner.Z())&&
687
765
(v_max.X()>m_min_corner.X())&&(v_max.Y()>m_min_corner.Y())&&(v_max.Z()>m_min_corner.Z())){
688
766
689
767
NTBlock<T> SearchBlock=block(v_min,v_max);
690
768
691
769
for(typename NTBlock<T>::iterator iter=SearchBlock.begin();
692
770
iter!=SearchBlock.end();
693
771
iter++){
696
774
}
697
775
}
698
776
}
699
777
700
778
return nlist;
701
779
}
702
780
703
781
/*!
704
Get list of all particles near a given point.
782
Get list of all particles near a given point.
705
783
706
784
\param p the point
707
*/
785
*/
708
786
template<typename T>
709
787
T_Handle<typename NeighborTable<T>::particlelist> NeighborTable<T>::getParticlesNearPoint(const Vec3& p)
710
788
{
716
794
// check if completely outside
717
795
if((v_min.X()<m_max_corner.X())&&(v_min.Y()<m_max_corner.Y())&&(v_min.Z()<m_max_corner.Z())&&
718
796
(v_max.X()>m_min_corner.X())&&(v_max.Y()>m_min_corner.Y())&&(v_max.Z()>m_min_corner.Z())){
719
797
720
798
NTBlock<T> SearchBlock=block(v_min,v_max);
721
799
722
800
for(typename NTBlock<T>::iterator iter=SearchBlock.begin();
723
801
iter!=SearchBlock.end();
724
802
iter++){
738
816
T_Handle<typename NeighborTable<T>::particlelist> NeighborTable<T>::getAllParticles()
739
817
{
740
818
T_Handle<typename NeighborTable<T>::particlelist> nlist=new typename NeighborTable<T>::particlelist;
741
819
742
820
for(typename list<T>::iterator iter=m_list.begin();
743
821
iter!=m_list.end();
744
822
iter++){
784
862
}
785
863
return ost;
786
864
}
865
866
867
868
/****fluid contents: begin****/
869
870
/*!
871
Adding fluid
872
*/
873
template<typename T>
874
void NeighborTable<T>::initiateFluid(
875
double Bw,
876
double Bp,
877
double Mu,
878
double alpha,
879
double flowrate,
880
double pressure,
881
Vec3 inflow,
882
Vec3 outflow,
883
int nx,
884
int ny,
885
int nz,
886
int nx_min,
887
int ny_min,
888
int nz_min
889
)
890
{
891
m_Bw=Bw;m_Bp=Bp;m_Mu=Mu;m_adjust=alpha;m_flowrate=flowrate;m_pressure=pressure;m_inflow=inflow;m_outflow=outflow;
892
m_xcell=nx;m_ycell=ny;m_zcell=nz;
893
m_global_cellidx=nx_min;m_global_cellidy=ny_min;m_global_cellidz=nz_min;
894
m_xside=m_dim*double(m_xsize-2)/double(nx-2);
895
m_yside=m_dim*double(m_ysize-2)/double(ny-2);
896
m_zside=m_dim*double(m_zsize-2)/double(nz-2);
897
/*
898
cout << "neighbor table info: "<<endl;
899
cout << " Bw=" << m_Bw << " Bp=" <<m_Bp << " Mu=" << m_Mu << " m_adjust=" << m_adjust << endl;
900
cout << " m_flowrate=" << m_flowrate << " m_pressure=" <<m_pressure << " m_inflow=" << m_inflow << " m_outflow=" << m_outflow<<endl;
901
cout << " m_xcell=" << m_xcell <<" m_ycell="<<m_ycell<<" m_zcell="<<m_zcell<< endl;
902
cout << " m_xside=" << m_xside <<" m_yside="<<m_yside<<" m_zside="<<m_zside<< endl;
903
*/
904
m_cellarray.resize(m_xcell*m_ycell*m_zcell);
905
m_cells.resize(m_xcell*m_ycell*m_zcell);
906
for(int i=0;i<m_xcell;i++){
907
for(int j=0;j<m_ycell;j++){
908
for(int k=0;k<m_zcell;k++){
909
double global_x=m_min_corner.X()+m_dim+(i-0.5)*m_xside;
910
double global_y=m_min_corner.Y()+m_dim+(j-0.5)*m_yside;
911
double global_z=m_min_corner.Z()+m_dim+(k-0.5)*m_zside;
912
Vec3 global_pos=Vec3(global_x,global_y,global_z);
913
int idx=cellindex(i,j,k);
914
m_cells[idx].setPos(global_pos);
915
m_cells[idx].setSize(Vec3(m_xside,m_yside,m_zside));
916
m_cells[idx].setIndex(Vec3(m_global_cellidx+i,m_global_cellidy+j,m_global_cellidz+k));
917
}
918
}
919
}
920
m_fluidexist=true;
921
}
922
923
924
/*!
925
Create a full list of pair of particle and cell and return handle to it
926
*/
927
template<typename T>
928
T_Handle<typename NeighborTable<T>::particlecelllist> NeighborTable<T>::getParticleCellList()
929
{
930
T_Handle<typename NeighborTable<T>::particlecelllist> nlist=new typename NeighborTable<T>::particlecelllist;
931
932
int number=0;
933
for(typename list<T>::iterator iter=m_list.begin();
934
iter!=m_list.end();
935
iter++){
936
int idx=cellindex(iter->getPos());
937
if(idx!=-1){
938
nlist->push_back(make_pair(&(*iter),&m_cells[idx]));
939
number++;
940
}
941
}
942
return nlist;
943
}
944
945
/*!
946
Create a new list of pair of particle and cell and return handle to it
947
*/
948
template<typename T>
949
T_Handle<typename NeighborTable<T>::particlecelllist> NeighborTable<T>::getNewParticleCellList()
950
{
951
T_Handle<typename NeighborTable<T>::particlecelllist> nlist=new typename NeighborTable<T>::particlecelllist;
952
for(typename list<T>::iterator iter=m_list.begin();
953
iter!=m_list.end();
954
iter++){
955
if(iter->isFlagged()){
956
int idx=cellindex(iter->getPos());
957
if(idx!=-1){
958
nlist->push_back(make_pair(&(*iter),&m_cells[idx]));
959
}
960
}
961
}
962
return nlist;
963
}
964
965
966
/*!
967
Calculate initial fluid cell porosity using measurement sphere method
968
*/
969
970
template<typename T>
971
void NeighborTable<T>::calPorositySphere0()
972
{
973
int number=0;
974
int num_cells=0;
975
double tot_d=0.0;
976
double tot_r3=0.0;
977
double r1=(m_xside+m_yside+m_zside)/6.0;
978
double cellVolume=3.1415926*4.0/3.0*r1*r1*r1;
979
for(int i=1;i<m_xcell-1;i++){
980
for(int j=1;j<m_ycell-1;j++){
981
for(int k=1;k<m_zcell-1;k++){
982
int idx=cellindex(i,j,k);
983
Vec3 pos=m_cells[idx].getPos();
984
double volume=0;
985
if(i==1||i==m_xcell-2||j==1||j==m_ycell-2||k==1||k==m_zcell-2){
986
int ix=int(floor((pos.X()-m_p0_global.X())/m_dim))-m_global_idx;
987
int iy=int(floor((pos.Y()-m_p0_global.Y())/m_dim))-m_global_idy;
988
int iz=int(floor((pos.Z()-m_p0_global.Z())/m_dim))-m_global_idz;
989
if((ix<1)||(ix>m_xsize-2)||(iy<1)||(iy>m_ysize-2)||(iz<1)||(iz>m_zsize-2)){
990
cout<<"neighbor table size does not match fluid cell size, program aborted!"<<endl;
991
abort();
992
}
993
double r=m_dim/2.0;
994
double gridVolume=3.1415926*4.0/3.0*r*r*r;
995
int ixmin=ix-1; int ixmax=ix+1;
996
int iymin=iy-1; int iymax=iy+1;
997
int izmin=iz-1; int izmax=iz+1;
998
int inumber=0;
999
double iPhi=0;
1000
for(int iix=ixmin;iix<=ixmax;iix++){
1001
for(int iiy=iymin;iiy<=iymax;iiy++){
1002
for(int iiz=izmin;iiz<=izmax;iiz++){
1003
if((iix>=1)&&(iix<=m_xsize-2)&&(iiy>=1)&&(iiy<=m_ysize-2)&&(iiz>=1)&&(iiz<=m_zsize-2)){
1004
double global_x=m_p0_global.X()+(iix+m_global_idx+0.5)*m_dim;
1005
double global_y=m_p0_global.Y()+(iiy+m_global_idy+0.5)*m_dim;
1006
double global_z=m_p0_global.Z()+(iiz+m_global_idz+0.5)*m_dim;
1007
Vec3 nt_pos=Vec3(global_x,global_y,global_z);
1008
int xmin=iix-1; int xmax=iix+1;
1009
int ymin=iiy-1; int ymax=iiy+1;
1010
int zmin=iiz-1; int zmax=iiz+1;
1011
double ivolume=0;
1012
for(int x=xmin;x<=xmax;x++){
1013
for(int y=ymin;y<=ymax;y++){
1014
for(int z=zmin;z<=zmax;z++){
1015
int nt_idx=index(x,y,z);
1016
for(typename pointtype::iterator iter=m_array[nt_idx].begin();
1017
iter!=m_array[nt_idx].end();
1018
iter++){
1019
double r2=(*iter)->getRad();
1020
double d=((*iter)->getPos()-nt_pos).norm();
1021
if(d<=r-r2){
1022
ivolume+=3.1415926*4.0/3.0*r2*r2*r2;
1023
} else if(d<r+r2){
1024
ivolume+=3.1415926/12.0/d*(r+r2-d)*(r+r2-d)*(d*d+2.0*d*(r+r2)-3.0*(r-r2)*(r-r2));
1025
}
1026
}
1027
}
1028
}
1029
}
1030
double phi=1.0-ivolume/gridVolume;
1031
if(phi<=0){phi=0.00001;};
1032
if(phi>=1){phi=0.99999;};
1033
iPhi+=phi;
1034
inumber++;
1035
}
1036
}
1037
}
1038
}
1039
1040
double newPhi=iPhi/double(inumber);
1041
if(newPhi<=0.0){newPhi=0.00001;};
1042
if(newPhi>=1.0){newPhi=0.99999;};
1043
m_cells[idx].setPhi(newPhi);
1044
m_cells[idx].setnewPhi(newPhi);
1045
m_cells[idx].seteffPhi(newPhi);
1046
} else {
1047
int xmin=i-1; int xmax=i+1;
1048
int ymin=j-1; int ymax=j+1;
1049
int zmin=k-1; int zmax=k+1;
1050
for(int x=xmin;x<=xmax;x++){
1051
for(int y=ymin;y<=ymax;y++){
1052
for(int z=zmin;z<=zmax;z++){
1053
int index=cellindex(x,y,z);
1054
for(typename pointtype::iterator iter=m_cellarray[index].begin();
1055
iter!=m_cellarray[index].end();
1056
iter++){
1057
double r2=(*iter)->getRad();
1058
double d=((*iter)->getPos()-pos).norm();
1059
if(d<=r1-r2){
1060
volume+=3.1415926*4.0/3.0*r2*r2*r2;
1061
} else if(d<r1+r2){
1062
volume+=3.1415926/12.0/d*(r1+r2-d)*(r1+r2-d)*(d*d+2.0*d*(r1+r2)-3.0*(r1-r2)*(r1-r2));
1063
}
1064
}
1065
}
1066
}
1067
}
1068
m_cells[idx].setVolume(volume);
1069
double newPhi=1.0-m_cells[idx].getVolume()/cellVolume;
1070
if(newPhi<=0.0){newPhi=0.00001;};
1071
if(newPhi>=1.0){newPhi=0.99999;};
1072
m_cells[idx].setPhi(newPhi);
1073
m_cells[idx].setnewPhi(newPhi);
1074
m_cells[idx].seteffPhi(newPhi);
1075
}
1076
int num_p=0;
1077
for(typename pointtype::iterator iter=m_cellarray[idx].begin();
1078
iter!=m_cellarray[idx].end();
1079
iter++) {
1080
double r=(*iter)->getRad();
1081
tot_d+=2.0*r;
1082
tot_r3+=pow(r,3.0);
1083
number++;
1084
num_p++;
1085
}
1086
if(num_p!=0){num_cells++;};
1087
}
1088
}
1089
}
1090
double meanK;
1091
if(number!=0){
1092
double meanD=tot_d/double(number);
1093
double tot_vp=4.0/3.0*3.1415926*tot_r3;
1094
double meanPhi=1.0-tot_vp/double(num_cells)/(m_xside*m_yside*m_zside);
1095
meanK=pow(meanD,2.0)/180.0*pow(meanPhi,3.0)/pow(1.0-meanPhi,2.0);
1096
} else{
1097
meanK=1e-7;
1098
}
1099
for(int i=1;i<m_xcell-1;i++){
1100
for(int j=1;j<m_ycell-1;j++){
1101
for(int k=1;k<m_zcell-1;k++){
1102
int idx=cellindex(i,j,k);
1103
double cell_totd=0.0;
1104
int p_in_cell=0;
1105
double newK;
1106
for(typename pointtype::iterator iter=m_cellarray[idx].begin();
1107
iter!=m_cellarray[idx].end();
1108
iter++) {
1109
double r=(*iter)->getRad();
1110
cell_totd+=2.0*r;
1111
p_in_cell++;
1112
}
1113
if(p_in_cell!=0){
1114
double D=cell_totd/double(p_in_cell);
1115
m_cells[idx].setD(D);//mean diameter of particles
1116
double Phi=m_cells[idx].getnewPhi();
1117
newK=pow(D,2.0)/180.0*pow(Phi,3.0)/pow(1.0-Phi,2.0);
1118
if(newK>meanK){newK=meanK;};
1119
} else {
1120
m_cells[idx].setD(0);
1121
newK=meanK;
1122
}
1123
m_cells[idx].setK(newK);
1124
m_cells[idx].seteffK(newK);
1125
double newBf=1.0/((1.0-m_cells[idx].getnewPhi())/m_Bp+m_cells[idx].getnewPhi()/m_Bw);
1126
m_cells[idx].setBf(newBf);
1127
m_cells[idx].seteffBf(newBf);
1128
m_cells[idx].setMu(m_Mu);
1129
}
1130
}
1131
}
1132
}
1133
1134
1135
/*!
1136
Calculate the fluid cell porosity using measurement sphere method
1137
*/
1138
1139
template<typename T>
1140
void NeighborTable<T>::calPorositySphere()
1141
{
1142
double r1=(m_xside+m_yside+m_zside)/6.0;
1143
double cellVolume=3.1415926*4.0/3.0*r1*r1*r1;
1144
1145
for(int i=1;i<m_xcell-1;i++){
1146
for(int j=1;j<m_ycell-1;j++){
1147
for(int k=1;k<m_zcell-1;k++){
1148
int idx=cellindex(i,j,k);
1149
Vec3 pos=m_cells[idx].getPos();
1150
double volume=0;
1151
1152
if(i==1||i==m_xcell-2||j==1||j==m_ycell-2||k==1||k==m_zcell-2){
1153
int ix=int(floor((pos.X()-m_p0_global.X())/m_dim))-m_global_idx;
1154
int iy=int(floor((pos.Y()-m_p0_global.Y())/m_dim))-m_global_idy;
1155
int iz=int(floor((pos.Z()-m_p0_global.Z())/m_dim))-m_global_idz;
1156
if((ix<1)||(ix>m_xsize-2)||(iy<1)||(iy>m_ysize-2)||(iz<1)||(iz>m_zsize-2)){
1157
cerr<<"neighbor table size does not match fluid cell size, program aborted!"<<endl;
1158
abort();
1159
}
1160
double r=m_dim/2.0;
1161
double gridVolume=3.1415926*4.0/3.0*r*r*r;
1162
int ixmin=ix-1; int ixmax=ix+1;
1163
int iymin=iy-1; int iymax=iy+1;
1164
int izmin=iz-1; int izmax=iz+1;
1165
int inumber=0;
1166
double iPhi=0;
1167
for(int iix=ixmin;iix<=ixmax;iix++){
1168
for(int iiy=iymin;iiy<=iymax;iiy++){
1169
for(int iiz=izmin;iiz<=izmax;iiz++){
1170
if((iix>=1)&&(iix<=m_xsize-2)&&(iiy>=1)&&(iiy<=m_ysize-2)&&(iiz>=1)&&(iiz<=m_zsize-2)){
1171
double global_x=m_p0_global.X()+(iix+m_global_idx+0.5)*m_dim;
1172
double global_y=m_p0_global.Y()+(iiy+m_global_idy+0.5)*m_dim;
1173
double global_z=m_p0_global.Z()+(iiz+m_global_idz+0.5)*m_dim;
1174
Vec3 nt_pos=Vec3(global_x,global_y,global_z);
1175
int xmin=iix-1; int xmax=iix+1;
1176
int ymin=iiy-1; int ymax=iiy+1;
1177
int zmin=iiz-1; int zmax=iiz+1;
1178
double ivolume=0;
1179
for(int x=xmin;x<=xmax;x++){
1180
for(int y=ymin;y<=ymax;y++){
1181
for(int z=zmin;z<=zmax;z++){
1182
int nt_idx=index(x,y,z);
1183
for(typename pointtype::iterator iter=m_array[nt_idx].begin();
1184
iter!=m_array[nt_idx].end();
1185
iter++){
1186
double r2=(*iter)->getRad();
1187
double d=((*iter)->getPos()-nt_pos).norm();
1188
if(d<=r-r2){
1189
ivolume+=3.1415926*4.0/3.0*r2*r2*r2;
1190
} else if(d<r+r2){
1191
ivolume+=3.1415926/12.0/d*(r+r2-d)*(r+r2-d)*(d*d+2.0*d*(r+r2)-3.0*(r-r2)*(r-r2));
1192
}
1193
}
1194
}
1195
}
1196
}
1197
double phi=1.0-ivolume/gridVolume;
1198
if(phi<=0){phi=0.00001;};
1199
if(phi>=1){phi=0.99999;};
1200
iPhi+=phi;
1201
inumber++;
1202
}
1203
}
1204
}
1205
}
1206
m_cells[idx].setPhi(m_cells[idx].getnewPhi());//porosity
1207
m_cells[idx].setnewPhi(iPhi/double(inumber));//new porosity
1208
if(m_cells[idx].getnewPhi()<=0){m_cells[idx].setnewPhi(0.00001);};
1209
if(m_cells[idx].getnewPhi()>=1){m_cells[idx].setnewPhi(0.99999);};
1210
} else {
1211
int xmin=i-1; int xmax=i+1;
1212
int ymin=j-1; int ymax=j+1;
1213
int zmin=k-1; int zmax=k+1;
1214
for(int x=xmin;x<=xmax;x++){
1215
for(int y=ymin;y<=ymax;y++){
1216
for(int z=zmin;z<=zmax;z++){
1217
int index=cellindex(x,y,z);
1218
for(typename pointtype::iterator iter=m_cellarray[index].begin();
1219
iter!=m_cellarray[index].end();
1220
iter++){
1221
double r2=(*iter)->getRad();
1222
double d=((*iter)->getPos()-pos).norm();
1223
if(d<=r1-r2){
1224
volume+=3.1415926*4.0/3.0*r2*r2*r2;
1225
} else if(d<r1+r2){
1226
volume+=3.1415926/12.0/d*(r1+r2-d)*(r1+r2-d)*(d*d+2.0*d*(r1+r2)-3.0*(r1-r2)*(r1-r2));
1227
}
1228
}
1229
}
1230
}
1231
}
1232
m_cells[idx].setVolume(volume);
1233
m_cells[idx].setPhi(m_cells[idx].getnewPhi());//porosity
1234
m_cells[idx].setnewPhi(1.0-m_cells[idx].getVolume()/cellVolume);//new porosity
1235
if(m_cells[idx].getnewPhi()<=0){m_cells[idx].setnewPhi(0.00001);};
1236
if(m_cells[idx].getnewPhi()>=1){m_cells[idx].setnewPhi(0.99999);};
1237
}
1238
1239
m_cells[idx].seteffPhi((1.0-m_adjust)*m_cells[idx].getPhi()+m_adjust*m_cells[idx].getnewPhi());//effective porosity
1240
}
1241
}
1242
}
1243
}
1244
1245
1246
/*!
1247
Update fluid cell: m_D, m_Vp, m_K, m_Bf;
1248
*/
1249
template<typename T>
1250
void NeighborTable<T>::updateFluidcells(){
1251
int number=0;
1252
int num_cells=0;
1253
double tot_d=0.0;
1254
double tot_r3=0.0;
1255
for(int i=1;i<m_xcell-1;i++){
1256
for(int j=1;j<m_ycell-1;j++){
1257
for(int k=1;k<m_zcell-1;k++){
1258
int idx=cellindex(i,j,k);
1259
int num_p=0;
1260
for(typename pointtype::iterator iter=m_cellarray[idx].begin();
1261
iter!=m_cellarray[idx].end();
1262
iter++) {
1263
double r=(*iter)->getRad();
1264
tot_d+=2.0*r;
1265
tot_r3+=pow(r,3.0);
1266
number++;
1267
num_p++;
1268
}
1269
if(num_p!=0){num_cells++;};
1270
}
1271
}
1272
}
1273
1274
double meanK;
1275
if(number!=0){
1276
double meanD=tot_d/double(number);
1277
double tot_vp=4.0/3.0*3.1415926*tot_r3;
1278
double meanPhi=1.0-tot_vp/double(num_cells)/(m_xside*m_yside*m_zside);
1279
meanK=pow(meanD,2.0)/180.0*pow(meanPhi,3.0)/pow(1.0-meanPhi,2.0);
1280
} else{
1281
meanK=1e-7;
1282
}
1283
for(int i=1;i<m_xcell-1;i++){
1284
for(int j=1;j<m_ycell-1;j++){
1285
for(int k=1;k<m_zcell-1;k++){
1286
int idx=cellindex(i,j,k);
1287
double cell_totd=0.0;
1288
int p_in_cell=0;
1289
double newK;
1290
Vec3 tot_vel;
1291
for(typename pointtype::iterator iter=m_cellarray[idx].begin();
1292
iter!=m_cellarray[idx].end();
1293
iter++) {
1294
double r=(*iter)->getRad();
1295
cell_totd+=2.0*r;
1296
p_in_cell++;
1297
}
1298
if(p_in_cell!=0){
1299
double D=cell_totd/double(p_in_cell);
1300
m_cells[idx].setD(D);//mean diameter of particles
1301
Vec3 Vp=tot_vel/double(p_in_cell);
1302
m_cells[idx].setVp(Vp);//mean particle velocity
1303
double Phi=m_cells[idx].getnewPhi();
1304
newK=pow(D,2.0)/180.0*pow(Phi,3.0)/pow(1.0-Phi,2.0);
1305
if(newK>meanK){newK=meanK;};
1306
} else {
1307
m_cells[idx].setD(0);
1308
m_cells[idx].setVp(Vec3(0.0,0.0,0.0));
1309
newK=meanK;
1310
}
1311
m_cells[idx].seteffK((1.0-m_adjust)*m_cells[idx].getK()+m_adjust*newK); //effective permeability
1312
m_cells[idx].setK(newK);//updated permeability
1313
double Phi=m_cells[idx].getnewPhi();
1314
double newBf=1.0/((1.0-Phi)/m_Bp+Phi/m_Bw);
1315
m_cells[idx].seteffBf((1.0-m_adjust)*m_cells[idx].getBf()+m_adjust*newBf);//effective bulk mudulus
1316
m_cells[idx].setBf(newBf);//updated bulk mudulus
1317
m_cells[idx].setMu(m_Mu);
1318
}
1319
}
1320
}
1321
}
1322
1323
1324
/*!
1325
Calculate coefficients for linear equations
1326
*/
1327
template<typename T>
1328
void NeighborTable<T>::calCoeffi(double timestep,int nt) {
1329
1330
double A_w, A_e, A_n, A_s, A_u, A_d, A_c, B_w, B_e, B_n, B_s, B_u, B_d, B_c;
1331
double co_w,co_e,co_n,co_s,co_u, co_d, co_c,co_b;
1332
int idx,idx_w,idx_e,idx_s,idx_n,idx_u,idx_d;
1333
double K,Kw,Ke,Ks,Kn,Ku,Kd;
1334
double P,Pw,Pe,Ps,Pn,Pu,Pd;
1335
double flowrate;
1336
double V=m_xside*m_yside*m_zside;
1337
for(int ix=1;ix<m_xcell-1;ix++){
1338
for(int iy=1;iy<m_ycell-1;iy++){
1339
for(int iz=1;iz<m_zcell-1;iz++){
1340
idx=cellindex(ix,iy,iz);
1341
idx_w=cellindex(ix-1,iy,iz); idx_e=cellindex(ix+1,iy,iz); idx_n=cellindex(ix,iy-1,iz);
1342
idx_s=cellindex(ix,iy+1,iz); idx_d=cellindex(ix,iy,iz-1); idx_u=cellindex(ix,iy,iz+1);
1343
K=m_cells[idx].geteffK();
1344
Kw=m_cells[idx_w].geteffK(); Ke=m_cells[idx_e].geteffK(); Ks=m_cells[idx_s].geteffK();
1345
Kn=m_cells[idx_n].geteffK(); Ku=m_cells[idx_u].geteffK(); Kd=m_cells[idx_d].geteffK();
1346
P=m_cells[idx].getdisP();
1347
Pw=m_cells[idx_w].getdisP(); Pe=m_cells[idx_e].getdisP(); Ps=m_cells[idx_s].getdisP();
1348
Pn=m_cells[idx_n].getdisP(); Pu=m_cells[idx_u].getdisP(); Pd=m_cells[idx_d].getdisP();
1349
double Bf=m_cells[idx].geteffBf();
1350
double Phi=m_cells[idx].geteffPhi();
1351
double dV=(m_cells[idx].getPhi()-m_cells[idx].getnewPhi())/nt; //pore volume change
1352
double dP=Bf*dV; //pressure change
1353
if(fabs(m_min_corner.X()+m_dim-m_p0_global.X())<1e-6 && ix==1){ //cell located at west boundary
1354
if(m_inflow.X()==1 || m_inflow.X()==2){ //inlet
1355
A_w=0;
1356
co_w=0;
1357
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_xside;}
1358
else{flowrate=m_flowrate;};
1359
B_w=Bf/V/Phi*flowrate*timestep*m_yside*m_zside;
1360
} else if(m_outflow.X()==-1 || m_outflow.X()==2){ //outlet
1361
A_w=K*m_xside*timestep/m_Mu;
1362
co_w=0;
1363
B_w=0;
1364
} else { //closed boundary
1365
A_w=0;co_w=0;B_w=0;
1366
}
1367
} else{ //cell located in inner area
1368
A_w=0.5*(K+Kw)*timestep*m_yside*m_zside/m_xside/m_Mu;
1369
co_w=Bf/V/Phi*A_w;
1370
B_w=(1.0-m_adjust)*co_w*Pw;
1371
}
1372
1373
if(fabs(m_max_corner.X()-m_dim-m_pmax_global.X())<1e-6 && ix==(m_xcell-2)){ //cell located at east boundary
1374
if(m_inflow.X()==-1 || m_inflow.X()==2){ //inlet
1375
A_e=0;
1376
co_e=0;
1377
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_xside;}
1378
else{flowrate=m_flowrate;};
1379
B_e=Bf/V/Phi*flowrate*timestep*m_yside*m_zside;
1380
} else if(m_outflow.X()==1 || m_outflow.X()==2){ //outlet
1381
A_e=K*m_xside*timestep/m_Mu;
1382
co_e=0;
1383
B_e=0;
1384
} else { //closed boundary
1385
A_e=0;co_e=0;B_e=0;
1386
}
1387
} else{ //cell located in inner area
1388
A_e=0.5*(K+Ke)*timestep*m_yside*m_zside/m_xside/m_Mu;
1389
co_e=Bf/V/Phi*A_e;
1390
B_e=(1.0-m_adjust)*co_e*Pe;
1391
}
1392
1393
if(fabs(m_min_corner.Y()+m_dim-m_p0_global.Y())<1e-6 && iy==1){ //cell located at north boundary
1394
if(m_inflow.Y()==1 || m_inflow.Y()==2){ //inlet
1395
A_n=0;
1396
co_n=0;
1397
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_yside;}
1398
else{flowrate=m_flowrate;};
1399
B_n=Bf/V/Phi*flowrate*timestep*m_xside*m_zside;
1400
} else if(m_outflow.Y()==-1 || m_outflow.Y()==2){ //outlet
1401
A_n=K*m_yside*timestep/m_Mu;
1402
co_n=0;
1403
B_n=0;
1404
} else { //closed boundary
1405
A_n=0;co_n=0;B_n=0;
1406
}
1407
} else{ //cell located in inner area
1408
A_n=0.5*(K+Kn)*timestep*m_xside*m_zside/m_yside/m_Mu;
1409
co_n=Bf/V/Phi*A_n;
1410
B_n=(1.0-m_adjust)*co_n*Pn;
1411
}
1412
1413
if(fabs(m_max_corner.Y()-m_dim-m_pmax_global.Y())<1e-6 && iy==(m_ycell-2)){ //cell located at south boundary
1414
if(m_inflow.Y()==-1 || m_inflow.Y()==2){ //inlet
1415
A_s=0;
1416
co_s=0;
1417
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_yside;}
1418
else{flowrate=m_flowrate;};
1419
B_s=Bf/V/Phi*flowrate*timestep*m_xside*m_zside;
1420
} else if(m_outflow.Y()==1 || m_outflow.Y()==2){ //outlet
1421
A_s=K*m_yside*timestep/m_Mu;
1422
co_s=0;
1423
B_s=0;
1424
} else { //closed boundary
1425
A_s=0;co_s=0;B_s=0;
1426
}
1427
} else{ //cell located in inner area
1428
A_s=0.5*(K+Ks)*timestep*m_xside*m_zside/m_yside/m_Mu;
1429
co_s=Bf/V/Phi*A_s;
1430
B_s=(1.0-m_adjust)*co_s*Ps;
1431
}
1432
1433
if(fabs(m_min_corner.Z()+m_dim-m_p0_global.Z())<1e-6 && iz==1){ //cell located at down/lower boundary
1434
if(m_inflow.Z()==1 || m_inflow.Z()==2){ //inlet
1435
A_d=0;
1436
co_d=0;
1437
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_zside;}
1438
else{flowrate=m_flowrate;};
1439
B_d=Bf/V/Phi*flowrate*timestep*m_xside*m_yside;
1440
} else if(m_outflow.Z()==-1 || m_outflow.Z()==2){ //outlet
1441
A_d=K*m_zside*timestep/m_Mu;
1442
co_d=0;
1443
B_d=0;
1444
} else { //closed boundary
1445
A_d=0;co_d=0;B_d=0;
1446
}
1447
} else{ //cell located in inner area
1448
A_d=0.5*(K+Kd)*timestep*m_xside*m_yside/m_zside/m_Mu;
1449
co_d=Bf/V/Phi*A_d;
1450
B_d=(1.0-m_adjust)*co_d*Pd;
1451
}
1452
1453
if(fabs(m_max_corner.Z()-m_dim-m_pmax_global.Z())<1e-6 && iz==(m_zcell-2)){ //cell located at upper boundary
1454
if(m_inflow.Z()==-1 || m_inflow.Z()==2){ //inlet
1455
A_u=0;
1456
co_u=0;
1457
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_zside;}
1458
else{flowrate=m_flowrate;};
1459
B_u=Bf/V/Phi*flowrate*timestep*m_xside*m_yside;
1460
} else if(m_outflow.Z()==1 || m_outflow.Z()==2){ //outlet
1461
A_u=K*m_zside*timestep/m_Mu;
1462
co_u=0;
1463
B_u=0;
1464
} else { //closed boundary
1465
A_u=0;co_u=0;B_u=0;
1466
}
1467
} else{ //cell located in inner area
1468
A_u=0.5*(K+Ku)*timestep*m_xside*m_yside/m_zside/m_Mu;
1469
co_u=Bf/V/Phi*A_u;
1470
B_u=(1.0-m_adjust)*co_u*Pu;
1471
}
1472
A_c=-(A_w+A_e+A_n+A_s+A_d+A_u);
1473
co_c=m_adjust*Bf/V/Phi*A_c-1.0;
1474
B_c=((1.0-m_adjust)*A_c*Bf/V/Phi+1.0)*P; //coefficent for central cell
1475
co_b=-(B_w+B_e+B_n+B_s+B_d+B_u+B_c+dP); //Constant
1476
m_cells[idx].setc_W(m_adjust*co_w);
1477
m_cells[idx].setc_E(m_adjust*co_e);
1478
m_cells[idx].setc_N(m_adjust*co_n);
1479
m_cells[idx].setc_S(m_adjust*co_s);
1480
m_cells[idx].setc_D(m_adjust*co_d);
1481
m_cells[idx].setc_U(m_adjust*co_u);
1482
m_cells[idx].setc_C(co_c);
1483
m_cells[idx].setc_B(co_b);
1484
}
1485
}
1486
}
1487
}
1488
1489
1490
/*!
1491
setting distributed pore pressures;
1492
*/
1493
template<typename T>
1494
void NeighborTable<T>::setPressure(vector<pair<Vec3,double> > pressure)
1495
{
1496
for(int ix=1;ix<m_xcell-1;ix++){
1497
for(int iy=1;iy<m_ycell-1;iy++){
1498
for(int iz=1;iz<m_zcell-1;iz++){
1499
int idx=cellindex(ix,iy,iz);
1500
Vec3 global_index=m_cells[idx].getIndex();
1501
for(vector<pair<Vec3,double> >::iterator iter=pressure.begin();
1502
iter!=pressure.end();
1503
iter++){
1504
if(global_index==iter->first){
1505
1506
if(iter->second<0){
1507
m_cells[idx].setdisP(0);
1508
} else {
1509
m_cells[idx].setdisP(iter->second);
1510
};
1511
break;
1512
};
1513
}
1514
}
1515
}
1516
}
1517
calVelocity();
1518
}
1519
1520
1521
/*!
1522
calculate fluid velocity and pressure gradient
1523
*/
1524
template<typename T>
1525
void NeighborTable<T>::calVelocity()
1526
{
1527
int idx,idx_w,idx_e,idx_s,idx_n,idx_u,idx_d;
1528
double K,Kw,Ke,Ks,Kn,Ku,Kd;
1529
double P,Pw,Pe,Ps,Pn,Pu,Pd;
1530
double Vw,Ve,Vn,Vs,Vd,Vu;
1531
double flowrate;
1532
1533
for(int ix=1;ix<m_xcell-1;ix++){
1534
for(int iy=1;iy<m_ycell-1;iy++){
1535
for(int iz=1;iz<m_zcell-1;iz++){
1536
idx=cellindex(ix,iy,iz);
1537
idx_w=cellindex(ix-1,iy,iz); idx_e=cellindex(ix+1,iy,iz); idx_n=cellindex(ix,iy-1,iz);
1538
idx_s=cellindex(ix,iy+1,iz); idx_d=cellindex(ix,iy,iz-1); idx_u=cellindex(ix,iy,iz+1);
1539
K=m_cells[idx].getK();
1540
Kw=m_cells[idx_w].getK(); Ke=m_cells[idx_e].getK(); Ks=m_cells[idx_s].getK();
1541
Kn=m_cells[idx_n].getK(); Ku=m_cells[idx_u].getK(); Kd=m_cells[idx_d].getK();
1542
P=m_cells[idx].getdisP();
1543
Pw=m_cells[idx_w].getdisP(); Pe=m_cells[idx_e].getdisP(); Ps=m_cells[idx_s].getdisP();
1544
Pn=m_cells[idx_n].getdisP(); Pu=m_cells[idx_u].getdisP(); Pd=m_cells[idx_d].getdisP();
1545
double Phi=m_cells[idx].getnewPhi();
1546
//calculate fluid velocity:
1547
if(fabs(m_min_corner.X()+m_dim-m_p0_global.X())<1e-6 && ix==1){ //cell located at west boundary
1548
if(m_inflow.X()==1 || m_inflow.X()==2){
1549
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_xside;}
1550
else{flowrate=m_flowrate;};
1551
Vw=flowrate/Phi;Pw=flowrate*m_Mu*m_xside/K+P;
1552
} //inlet
1553
else if(m_outflow.X()==-1 || m_outflow.X()==2){
1554
Vw=K*(0-P)/m_xside/m_Mu/Phi;Pw=0;
1555
} //outlet
1556
else {Vw=999999;Pw=999999;} //closed boundary
1557
} else {
1558
Vw=0.5*(K+Kw)*(P-Pw)/m_xside/m_Mu/Phi;
1559
}
1560
1561
if(fabs(m_max_corner.X()-m_dim-m_pmax_global.X())<1e-6 && ix==(m_xcell-2)){ //cell located at east boundary
1562
if(m_inflow.X()==-1 || m_inflow.X()==2){
1563
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_xside;}
1564
else{flowrate=m_flowrate;};
1565
Ve=flowrate/Phi;Pe=flowrate*m_Mu*m_xside/K+P;
1566
} //inlet
1567
else if(m_outflow.X()==1 || m_outflow.X()==2){
1568
Ve=K*(0-P)/m_xside/m_Mu/Phi;Pe=0;
1569
} //outlet
1570
else {Ve=999999;Pe=999999;} //closed boundary
1571
} else {
1572
Ve=0.5*(K+Ke)*(P-Pe)/m_xside/m_Mu/Phi;
1573
}
1574
1575
if(fabs(m_min_corner.Y()+m_dim-m_p0_global.Y())<1e-6 && iy==1){ //cell located at north boundary
1576
if(m_inflow.Y()==1 || m_inflow.Y()==2){
1577
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_yside;}
1578
else{flowrate=m_flowrate;};
1579
Vn=flowrate/Phi;Pn=flowrate*m_Mu*m_yside/K+P;
1580
} //inlet
1581
else if(m_outflow.Y()==-1 || m_outflow.Y()==2){
1582
Vn=K*(0-P)/m_yside/m_Mu/Phi;Pn=0;
1583
} //outlet
1584
else {Vn=999999;Pn=999999;} //closed boundary
1585
} else {
1586
Vn=0.5*(K+Kn)*(P-Pn)/m_yside/m_Mu/Phi;
1587
}
1588
1589
if(fabs(m_max_corner.Y()-m_dim-m_pmax_global.Y())<1e-6 && iy==(m_ycell-2)){ //cell located at south boundary
1590
if(m_inflow.Y()==-1 || m_inflow.Y()==2){
1591
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_yside;}
1592
else{flowrate=m_flowrate;};
1593
Vs=flowrate/Phi;Ps=flowrate*m_Mu*m_yside/K+P;
1594
} //inlet
1595
else if(m_outflow.Y()==1 || m_outflow.Y()==2){
1596
Vs=K*(0-P)/m_yside/m_Mu/Phi;Ps=0;
1597
} //outlet
1598
else {Vs=999999;Ps=999999;} //closed boundary
1599
} else {
1600
Vs=0.5*(K+Ks)*(P-Ps)/m_yside/m_Mu/Phi;
1601
}
1602
1603
if(fabs(m_min_corner.Z()+m_dim-m_p0_global.Z())<1e-6 && iz==1){ //cell located at lower boundary
1604
if(m_inflow.Z()==1 || m_inflow.Z()==2){
1605
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_zside;}
1606
else{flowrate=m_flowrate;};
1607
Vd=flowrate/Phi;Pd=flowrate*m_Mu*m_zside/K+P;
1608
} //inlet
1609
else if(m_outflow.Z()==-1 || m_outflow.Z()==2){
1610
Vd=K*(0-P)/m_zside/m_Mu/Phi;Pd=0;
1611
} //outlet
1612
else {Vd=999999;Pd=999999;} //closed boundary
1613
} else {
1614
Vd=0.5*(K+Kd)*(P-Pd)/m_zside/m_Mu/Phi;
1615
}
1616
1617
if(fabs(m_max_corner.Z()-m_dim-m_pmax_global.Z())<1e-6 && iz==(m_zcell-2)){ //cell located at upper boundary
1618
if(m_inflow.Z()==-1 || m_inflow.Z()==2){
1619
if(m_flowrate==0){flowrate=K*m_pressure/m_Mu/m_zside;}
1620
else{flowrate=m_flowrate;};
1621
Vu=flowrate/Phi;Pu=flowrate*m_Mu*m_zside/K+P;
1622
} //inlet
1623
else if(m_outflow.Z()==1 || m_outflow.Z()==2){
1624
Vu=K*(0-P)/m_zside/m_Mu/Phi;Pu=0;
1625
} //outlet
1626
else {Vu=999999;Pu=999999;} //closed boundary
1627
} else {
1628
Vu=0.5*(K+Ku)*(P-Pu)/m_zside/m_Mu/Phi;
1629
}
1630
1631
double Vx,Vy,Vz,Px,Py,Pz;
1632
if(Vw==999999 && Ve==999999){
1633
Vx=0;Px=0;
1634
} else if(Vw==999999){
1635
Vx=Ve;Px=-Pe/m_xside;
1636
} else if(Ve==999999){
1637
Vx=-Vw;Px=Pw/m_xside;
1638
} else {
1639
Vx=(-Vw+Ve)/2.0;Px=(Pw-Pe)/2.0/m_xside;
1640
}
1641
1642
if(Vn==999999 && Vs==999999){
1643
Vy=0;Py=0;
1644
} else if(Vn==999999){
1645
Vy=Vs;Py=-Ps/m_yside;
1646
} else if(Vs==999999){
1647
Vy=-Vn;Py=Pn/m_yside;
1648
} else {
1649
Vy=(-Vn+Vs)/2.0;Py=(Pn-Ps)/2.0/m_yside;
1650
}
1651
1652
if(Vd==999999 && Vu==999999){
1653
Vz=0;Pz=0;
1654
} else if(Vd==999999){
1655
Vz=Vu;Pz=-Pu/m_zside;
1656
} else if(Vu==999999){
1657
Vz=-Vd;Pz=Pd/m_zside;
1658
} else {
1659
Vz=(-Vd+Vu)/2.0;Pz=(Pd-Pu)/2.0/m_zside;
1660
}
1661
m_cells[idx].setVf(Vec3(Vx,Vy,Vz));//fluid velocity
1662
m_cells[idx].setPg(Vec3(Px,Py,Pz));//pressure gradient
1663
}
1664
}
1665
}
1666
}
1667
1668
1669
template<typename T>
1670
vector<CFluidCell> NeighborTable<T>::get_yz_cellslab(int x)
1671
{
1672
vector<CFluidCell> res;
1673
for(int iy=1;iy<m_ycell-1;iy++){
1674
for(int iz=1;iz<m_zcell-1;iz++){
1675
int idx=cellindex(x,iy,iz);
1676
res.push_back(m_cells[idx]);
1677
}
1678
}
1679
return res;
1680
}
1681
1682
1683
template<typename T>
1684
vector<CFluidCell> NeighborTable<T>::get_xz_cellslab(int y)
1685
{
1686
vector<CFluidCell> res;;
1687
for(int ix=1;ix<m_xcell-1;ix++){
1688
for(int iz=1;iz<m_zcell-1;iz++){
1689
int idx=cellindex(ix,y,iz);
1690
res.push_back(m_cells[idx]);
1691
}
1692
}
1693
return res;
1694
}
1695
1696
1697
template<typename T>
1698
vector<CFluidCell> NeighborTable<T>::get_xy_cellslab(int z)
1699
{
1700
vector<CFluidCell> res;
1701
for(int ix=1;ix<m_xcell-1;ix++){
1702
for(int iy=1;iy<m_ycell-1;iy++){
1703
int idx=cellindex(ix,iy,z);
1704
res.push_back(m_cells[idx]);
1705
}
1706
}
1707
return res;
1708
}
1709
1710
1711
1712
template<typename T>
1713
void NeighborTable<T>::set_yz_cellslab(int x, vector<CFluidCell> recv_slab)
1714
{
1715
if(recv_slab.size()==(m_ycell-2)*(m_zcell-2)){
1716
int number=0;
1717
for(int iy=1;iy<m_ycell-1;iy++){
1718
for(int iz=1;iz<m_zcell-1;iz++){
1719
int idx=cellindex(x,iy,iz);
1720
m_cells[idx]=recv_slab[number];
1721
number++;
1722
}
1723
}
1724
} else{
1725
//cerr << "FluidcellTable<T>::set_yz_slab ERROR: size mismatch in received data, recv_slab.size()!=m_ycell*m_zcell - (" << recv_slab.size() << "!=" << m_ycell*m_zcell << ")"<< std::endl;
1726
}
1727
}
1728
1729
template<typename T>
1730
void NeighborTable<T>::set_xz_cellslab(int y, vector<CFluidCell> recv_slab)
1731
{
1732
if(recv_slab.size()==(m_xcell-2)*(m_zcell-2)){
1733
int number=0;
1734
for(int ix=1;ix<m_ycell-1;ix++){
1735
for(int iz=1;iz<m_zcell-1;iz++){
1736
int idx=cellindex(ix,y,iz);
1737
m_cells[idx]=recv_slab[number];
1738
number++;
1739
}
1740
}
1741
} else{
1742
//cerr << "FluidcellTable<T>::set_xz_slab ERROR: size mismatch in received data, recv_slab.size()!=m_xcell*m_zcell - (" << recv_slab.size() << "!=" << m_xcell*m_zcell << ")"<< std::endl;
1743
}
1744
}
1745
1746
template<typename T>
1747
void NeighborTable<T>::set_xy_cellslab(int z, vector<CFluidCell> recv_slab)
1748
{
1749
if(recv_slab.size()==(m_xcell-2)*(m_ycell-2)){
1750
int number=0;
1751
for(int ix=1;ix<m_xcell-1;ix++){
1752
for(int iy=1;iy<m_ycell-1;iy++){
1753
int idx=cellindex(ix,iy,z);
1754
m_cells[idx]=recv_slab[number];
1755
number++;
1756
}
1757
}
1758
} else{
1759
//cerr << "FluidcellTable<T>::set_xy_slab ERROR: size mismatch in received data, recv_slab.size()!=m_xcell*m_ycell - (" << recv_slab.size() << "!=" << m_xcell*m_ycell << ")"<< std::endl;
1760
}
1761
}
1762
1763
1764
/*!
1765
Get a value for each fluid cell using a fluid cell member function and return a vector
1766
of values, with position of cell.
1767
1768
\param rdf the fluid cell member function
1769
*/
1770
template<typename T>
1771
template<typename P>
1772
vector<pair<Vec3,P> > NeighborTable<T>::forAllInnerCellsGet(P (CFluidCell::*rdf)() const)
1773
{
1774
vector<pair<Vec3,P> > res;
1775
for(int ix=1;ix<m_xcell-1;ix++){
1776
for(int iy=1;iy<m_ycell-1;iy++){
1777
for(int iz=1;iz<m_zcell-1;iz++){
1778
int idx=cellindex(ix,iy,iz);
1779
Vec3 global_pos=m_cells[idx].getPos();
1780
res.push_back(make_pair(global_pos,(m_cells[idx].*rdf)()));
1781
}
1782
}
1783
}
1784
return res;
1785
}
1786
1787
1788
/*!
1789
Get a value for each fluid cell using a fluid cell member function and return a vector
1790
of values,with global index of cell
1791
1792
\param rdf the fluid cell member function
1793
*/
1794
template<typename T>
1795
template<typename P>
1796
vector<pair<Vec3,P> > NeighborTable<T>::forAllInnerCellsGetIndexed(P (CFluidCell::*rdf)() const)
1797
{
1798
vector<pair<Vec3,P> > res;
1799
for(int ix=1;ix<m_xcell-1;ix++){
1800
for(int iy=1;iy<m_ycell-1;iy++){
1801
for(int iz=1;iz<m_zcell-1;iz++){
1802
int idx=cellindex(ix,iy,iz);
1803
Vec3 global_index=m_cells[idx].getIndex();
1804
res.push_back(make_pair(global_index,(m_cells[idx].*rdf)()));
1805
}
1806
}
1807
}
1808
return res;
1809
}
1810
1811
1812
template<typename T>
1813
template<typename P>
1814
vector<P> NeighborTable<T>::forAllInnerCellsGetSum(P (CFluidCell::*rdf)() const)
1815
{
1816
vector<P> res;
1817
for(int ix=1;ix<m_xcell-1;ix++){
1818
for(int iy=1;iy<m_ycell-1;iy++){
1819
for(int iz=1;iz<m_zcell-1;iz++){
1820
int idx=cellindex(ix,iy,iz);
1821
res.push_back((m_cells[idx].*rdf)());
1822
}
1823
}
1824
}
1825
return res;
1826
}
1827
1828
/****fluid contents: end****/
1829
Loggerhead is a web-based interface for Breezy
Version: 2.0.1