File: N:\mfix\model\des\pic\apply_wall_bc_pic.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: APPLY_WALL_BC_PIC                                       !
4     !  Author: R. Garg                                                     !
5     !                                                                      !
6     !  Purpose: Detect collisions between PIC particles and STLs.          !
7     !                                                                      !
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9           SUBROUTINE APPLY_WALL_BC_PIC
10     
11     ! Global Variables:
12     !---------------------------------------------------------------------//
13     ! Particle postions: Current/Previous
14           use discretelement, only: DES_POS_NEW
15     ! Paricle velocities
16           use discretelement, only: DES_VEL_NEW
17     ! Particle radius
18           use discretelement, only: DES_RADIUS
19     ! Max number of particles on this process
20           use discretelement, only: MAX_PIP
21     ! Map from particle to DES grid cell.
22           use discretelement, only: DG_PIJK, DTSOLID
23     ! Flag indicating that the index cell contains no STLs
24           use stl, only: FACETS_AT_DG
25     ! STL Vertices
26           use stl, only: VERTEX
27     ! STL Facet normals
28           use stl, only: NORM_FACE
29     ! Minimum velocity to offset gravitational forces
30           use pic_bc, only: minVEL, minVEL_MAG, OoMinVEL_MAG
31     ! Solids time step size
32           use discretelement, only: DTSOLID
33     ! Gravitational force vector
34           use discretelement, only: GRAV
35     
36     ! Global Parameters:
37     !---------------------------------------------------------------------//
38           use param1, only: ZERO, SMALL_NUMBER, ONE, UNDEFINED
39           use functions, only: IS_NONEXISTENT
40     
41     ! Module Procedures:
42     !---------------------------------------------------------------------//
43           use stl_functions_des
44     
45           use error_manager
46     
47           IMPLICIT NONE
48     
49     ! Local Variables:
50     !---------------------------------------------------------------------//
51     ! Loop counters.
52           INTEGER NF, NP, LC1
53     ! STL normal vector
54           DOUBLE PRECISION :: NORM_PLANE(3)
55     ! line is parameterized as p = p_ref + t * dir_line, t is line_param
56           DOUBLE PRECISION :: LINE_T
57     
58     ! distance parcel travelled out of domain along the facet normal
59           DOUBLE PRECISION :: DIST(3), tPOS(3)
60     
61           DOUBLE PRECISION :: DT_RBND, RADSQ
62     !......................................................................!
63     
64           double precision :: vv(3)
65     
66     ! Minimum velocity needed to offset gravity.
67           minVEL = -DTSOLID*GRAV(:)
68           minVEL_MAG = dot_product(minVEL,minVEL)
69           OoMinVEL_MAG = 1.0d0/minVEL_MAG
70     
71           DT_RBND = 1.01d0*DTSOLID
72     
73           DO NP = 1, MAX_PIP
74     
75     ! Skip non-existent particles
76              IF(IS_NONEXISTENT(NP)) CYCLE
77     
78     ! If no neighboring facet in the surrounding 27 cells, then exit
79              IF(FACETS_AT_DG(DG_PIJK(NP))%COUNT < 1) CYCLE
80     
81              RADSQ = DES_RADIUS(NP)*DES_RADIUS(NP)
82     
83     ! Loop through the STLs associated with the DES grid cell.
84              LC1_LP: DO LC1=1, FACETS_AT_DG(DG_PIJK(NP))%COUNT
85     
86     ! Get the index of the neighbor facet
87                 NF = FACETS_AT_DG(DG_PIJK(NP))%ID(LC1)
88     
89                 NORM_PLANE = NORM_FACE(:,NF)
90     
91     !            IF(dot_product(NORM_PLANE, DES_VEL_NEW(NP,:)) > 0.d0) &
92     !               CYCLE LC1_LP
93     
94                 vv = VERTEX(1,:,NF)
95                 CALL INTERSECTLNPLANE(DES_POS_NEW(NP,:), DES_VEL_NEW(NP,:),&
96                     vv, NORM_FACE(:,NF), LINE_T)
97     
98                 IF(abs(LINE_T) <= DT_RBND) THEN
99     
100     ! Project the parcel onto the facet.
101                    DIST = LINE_T*DES_VEL_NEW(NP,:)
102                    tPOS = DES_POS_NEW(NP,:) + DIST
103     
104     ! Avoid collisions with STLs that are not next to the parcel
105                    IF(HIT_FACET(VERTEX(:,:,NF), tPOS)) THEN
106     
107     ! Correct the position of particles found too close to the STL.
108                       IF(dot_product(DIST,DIST) <= RADSQ .OR. &
109                          LINE_T <= ZERO) THEN
110                          DES_POS_NEW(NP,:) = DES_POS_NEW(NP,:) + DIST(:) + &
111                             DES_RADIUS(NP)*NORM_PLANE
112                       ENDIF
113     
114     ! Reflect the parcel.
115                       CALL PIC_REFLECT_PART(NP, NORM_PLANE(:))
116                    ENDIF
117                 ENDIF
118     
119              ENDDO LC1_LP
120           END DO
121     
122           RETURN
123     
124           contains
125     
126     
127     !......................................................................!
128           LOGICAL FUNCTION HIT_FACET(VERTS, POINT)
129     
130           DOUBLE PRECISION, INTENT(IN) :: VERTS(3,3)
131           DOUBLE PRECISION, INTENT(IN) :: POINT(3)
132     
133           DOUBLE PRECISION :: V0(3), V1(3), V2(3)
134           DOUBLE PRECISION :: d00, d01, d02, d11, d12
135     
136           DOUBLE PRECISION :: OoDEN, u, v
137     
138           V0 = VERTS(3,:) - VERTS(1,:)
139           V1 = VERTS(2,:) - VERTS(1,:)
140           V2 = POINT - VERTS(1,:)
141     
142           d00 = dot_product(V0,V0)
143           d01 = dot_product(V0,V1)
144           d02 = dot_product(V0,V2)
145     
146           d11 = dot_product(V1,V1)
147           d12 = dot_product(V1,V2)
148     
149           OoDEN = 1.0d0/(d00*d11 - d01*d01)
150           u = (d11*d02 - d01*d12)*OoDEN
151           v = (d00*d12 - d01*d02)*OoDEN
152     
153           HIT_FACET = ((u>=0.0d0) .AND. (v>=0.0d0) .AND. (u+v < 1.0d0))
154     
155           RETURN
156           END FUNCTION HIT_FACET
157     
158           END SUBROUTINE APPLY_WALL_BC_PIC
159     
160     
161     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
162     !                                                                      C
163     !  Subroutine: PIC_REFLECT_PARTICLE                                    C
164     !  Purpose:                                                            C
165     !                                                                      C
166     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
167           SUBROUTINE PIC_REFLECT_PART(LL, WALL_NORM)
168     
169           USE discretelement, only : dimn, DES_VEL_NEW
170           USE mfix_pic, only : MPPIC_COEFF_EN_WALL, MPPIC_COEFF_ET_WALL
171     
172     ! Minimum velocity to offset gravitational forces
173           use pic_bc, only: minVEL, minVEL_MAG, OoMinVEL_MAG
174     
175     
176           use param1, only: SMALL_NUMBER
177     
178           IMPLICIT NONE
179     !-----------------------------------------------
180     ! Dummy arguments
181     !-----------------------------------------------
182           INTEGER, INTENT(IN) :: LL
183           DOUBLE PRECISION, INTENT(IN) ::  WALL_NORM(DIMN)
184     !-----------------------------------------------
185     ! Local variables
186     !-----------------------------------------------
187           !magnitude of pre-collisional normal and tangential velocity components
188           DOUBLE PRECISION :: VEL_NORMMAG_APP
189     
190           !pre collisional normal and tangential velocity components in vector format
191           !APP ==> approach
192           DOUBLE PRECISION :: VEL_NORM_APP(DIMN), VEL_TANG_APP(DIMN)
193     
194     
195           !post collisional normal and tangential velocity components in vector format
196           !SEP ==> separation
197           DOUBLE PRECISION :: VEL_NORM_SEP(DIMN), VEL_TANG_SEP(DIMN)
198     
199           DOUBLE PRECISION :: COEFF_REST_EN, COEFF_REST_ET
200     
201     ! Minimum particle velocity with respect to gravity. This ensures that
202     ! walls "push back" to offset gravitational forces.
203           DOUBLE PRECISION :: projGRAV(3)
204     
205           INTEGER :: LC
206     
207     !-----------------------------------------------
208     
209     
210     
211     
212     
213           VEL_NORMMAG_APP = DOT_PRODUCT(WALL_NORM(:), DES_VEL_NEW(LL,:))
214     
215     ! Currently assuming that wall is at rest. Needs improvement for moving wall
216     
217           VEL_NORM_APP(:) = VEL_NORMMAG_APP*WALL_NORM(:)
218           VEL_TANG_APP(:) = DES_VEL_NEW(LL,:) - VEL_NORM_APP(:)
219     
220     
221     !post collisional velocities
222     
223           COEFF_REST_EN = MPPIC_COEFF_EN_WALL
224           COEFF_REST_ET = MPPIC_COEFF_ET_WALL
225     
226           VEL_NORM_SEP(:) = -COEFF_REST_EN*VEL_NORM_APP(:)
227           VEL_TANG_SEP(:) =  COEFF_REST_ET*VEL_TANG_APP(:)
228     
229           DES_VEL_NEW(LL,:) = VEL_NORM_SEP(:) + VEL_TANG_SEP(:)
230     
231           IF(dot_product(WALL_NORM, minVEL) > 0.0d0)THEN
232     
233     ! Projection of rebound velocity onto mininum velocity.
234              projGRAV = dot_product(minVEL, DES_VEL_NEW(LL,:))*OoMinVEL_MAG
235              projGRAV = minVEL*projGRAV
236     
237              IF(dot_product(projGRAV, projGRAV) < minVEL_MAG) THEN
238                 DO LC=1,3
239                    IF(minVEL(LC) > SMALL_NUMBER) THEN
240                       DES_VEL_NEW(LL,LC)=max(DES_VEL_NEW(LL,LC),minVEL(LC))
241                    ELSEIF(minVEL(LC) < -SMALL_NUMBER) THEN
242                       DES_VEL_NEW(LL,LC)=min(DES_VEL_NEW(LL,LC),minVEL(LC))
243                    ENDIF
244                 ENDDO
245     
246              ENDIF
247           ENDIF
248     
249           RETURN
250           END SUBROUTINE PIC_REFLECT_PART
251     
252     
253     !
254           SUBROUTINE OUT_THIS_STL(LC)
255     
256           Use usr
257     
258     ! STL Vertices
259           use stl, only: VERTEX
260     ! STL Facet normals
261           use stl, only: NORM_FACE
262     
263     
264           IMPLICIT NONE
265     !-----------------------------------------------
266           integer, INTENT(IN) :: lc
267     
268           logical :: lExists
269           character(len=8) :: IDX
270     
271     
272           write(idx,"(I8.8)") LC
273           inquire(file='geo_'//idx//'.stl', EXIST=lExists)
274     
275           IF(lExists) RETURN
276     
277           open(unit=555, file='geo_'//idx//'.stl', status='UNKNOWN')
278           write(555,*) 'solid vcg'
279           write(555,*) '   facet normal ', NORM_FACE(:,LC)
280           write(555,*) '      outer loop'
281           write(555,*) '         vertex ', VERTEX(1,1:3,LC)
282           write(555,*) '         vertex ', VERTEX(2,1:3,LC)
283           write(555,*) '         vertex ', VERTEX(3,1:3,LC)
284           write(555,*) '      endloop'
285           write(555,*) '   endfacet'
286           close(555)
287     
288           RETURN
289           END SUBROUTINE OUT_THIS_STL
290     
291     
292     
293     
294     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
295     !                                                                      C
296     !  SUBROUTINE: CHECK_IF_PARCEL_OVERLAPS_STL                            C
297     !  Authors: Rahul Garg                               Date: 21-Mar-2014 C
298     !                                                                      C
299     !  Purpose: This subroutine is special written to check if a particle  C
300     !          overlaps any of the STL faces. The routine exits on         C
301     !          detecting an overlap. It is called after initial            C
302     !          generation of lattice configuration to remove out of domain C
303     !          particles                                                   C
304     !                                                                      C
305     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
306           SUBROUTINE CHECK_IF_PARCEL_OVERLAPS_STL(POS, OVERLAP_EXISTS)
307     
308           USE discretelement, only: MAX_RADIUS
309           USE geometry, only: do_k
310     
311           USE stl_functions_des
312           USE desgrid
313           USE stl
314     
315           Implicit none
316     
317           DOUBLE PRECISION, INTENT(IN) :: POS(3)
318           LOGICAL, INTENT(OUT) :: OVERLAP_EXISTS
319     
320           INTEGER I, J, K, IJK, NF, LC
321     
322     ! line is parameterized as p = p_ref + t * dir_line, t is line_param
323           DOUBLE PRECISION :: LINE_T
324           DOUBLE PRECISION :: DIST(3), RADSQ
325           double precision :: vv(3)
326     
327           K = 1
328           IF(DO_K) K = min(DG_KEND2,max(DG_KSTART2,KOFPOS(POS(3))))
329           J = min(DG_JEND2,max(DG_JSTART2,JOFPOS(POS(2))))
330           I = min(DG_IEND2,max(DG_ISTART2,IOFPOS(POS(1))))
331     
332           IJK = DG_FUNIJK(I,J,K)
333           IF(FACETS_AT_DG(IJK)%COUNT < 1) RETURN
334     
335           OVERLAP_EXISTS = .TRUE.
336     
337     ! Pad the max radius to keep parcels sufficiently away from walls.
338           RADSQ = (1.1d0*MAX_RADIUS)**2
339     
340     ! Parametrize a line as p = p_0 + t normal and intersect with the
341     ! triangular plane.
342     
343     ! If t>0, then point is on the non-fluid side of the plane, if the
344     ! plane normal is assumed to point toward the fluid side.
345           DO LC = 1, FACETS_AT_DG(IJK)%COUNT
346              NF = FACETS_AT_DG(IJK)%ID(LC)
347     
348              vv = VERTEX(1,:,NF)
349              CALL INTERSECTLNPLANE(POS, NORM_FACE(:,NF), &
350                 vv, NORM_FACE(:,NF), LINE_T)
351     
352     ! Orthogonal projection puts the point outside of the domain or less than
353     ! one particle radius to the facet.
354              DIST = LINE_T*NORM_FACE(:,NF)
355              IF(LINE_T > ZERO .OR. dot_product(DIST,DIST) <= RADSQ) RETURN
356     
357           ENDDO
358     
359           OVERLAP_EXISTS = .FALSE.
360     
361           RETURN
362           END SUBROUTINE CHECK_IF_PARCEL_OVERLAPS_STL
363     
364     
365     !----------------------------------------------------------------------!
366     !                                                                      !
367     !                                                                      !
368     !                                                                      !
369     !                                                                      !
370     !----------------------------------------------------------------------!
371           SUBROUTINE write_this_facet_and_parcel(FID, position, velocity)
372           USE run
373           USE param1
374           USE discretelement
375           USE geometry
376           USE compar
377           USE constant
378           USE cutcell
379           USE funits
380           USE indices
381           USE physprop
382           USE parallel
383           USE stl
384           USE stl_functions_des
385           Implicit none
386           !facet id and particle id
387           double precision, intent(in), dimension(dimn) :: position, velocity
388           Integer, intent(in) :: fid
389           Integer :: stl_unit, vtp_unit , k
390           CHARACTER(LEN=100) :: stl_fname, vtp_fname
391           real :: temp_array(3)
392     
393           stl_unit = 1001
394           vtp_unit = 1002
395     
396           WRITE(vtp_fname,'(A,"_OFFENDING_PARTICLE",".vtp")') TRIM(RUN_NAME)
397           WRITE(stl_fname,'(A,"_STL_FACE",".stl")') TRIM(RUN_NAME)
398     
399           open(vtp_unit, file = vtp_fname, form='formatted',convert='big_endian')
400           open(stl_unit, file = stl_fname, form='formatted',convert='big_endian')
401     
402           write(vtp_unit,"(a)") '<?xml version="1.0"?>'
403           write(vtp_unit,"(a,es24.16,a)") '<!-- time =',s_time,'s -->'
404           write(vtp_unit,"(a,a)") '<VTKFile type="PolyData"',&
405                ' version="0.1" byte_order="LittleEndian">'
406           write(vtp_unit,"(3x,a)") '<PolyData>'
407           write(vtp_unit,"(6x,a,i10.10,a,a)")&
408                '<Piece NumberOfPoints="',1,'" NumberOfVerts="0" ',&
409                'NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0">'
410           write(vtp_unit,"(9x,a)")&
411                '<PointData Scalars="Diameter" Vectors="Velocity">'
412           write(vtp_unit,"(12x,a)")&
413                '<DataArray type="Float32" Name="Diameter" format="ascii">'
414           write (vtp_unit,"(15x,es13.6)") (1.d0)
415           write(vtp_unit,"(12x,a)") '</DataArray>'
416     
417           temp_array = zero
418           temp_array(1:DIMN) = velocity(1:dimn)
419           write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
420                'Name="Velocity" NumberOfComponents="3" format="ascii">'
421           write (vtp_unit,"(15x,3(es13.6,3x))")&
422                ((temp_array(k)),k=1,3)
423           write(vtp_unit,"(12x,a,/9x,a)") '</DataArray>','</PointData>'
424           ! skip cell data
425           write(vtp_unit,"(9x,a)") '<CellData></CellData>'
426     
427           temp_array = zero
428           temp_array(1:dimn) = position(1:dimn)
429           write(vtp_unit,"(9x,a)") '<Points>'
430           write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
431                'Name="Position" NumberOfComponents="3" format="ascii">'
432           write (vtp_unit,"(15x,3(es13.6,3x))")&
433                ((temp_array(k)),k=1,3)
434           write(vtp_unit,"(12x,a,/9x,a)")'</DataArray>','</Points>'
435           ! Write tags for data not included (vtp format style)
436           write(vtp_unit,"(9x,a,/9x,a,/9x,a,/9x,a)")'<Verts></Verts>',&
437                '<Lines></Lines>','<Strips></Strips>','<Polys></Polys>'
438           write(vtp_unit,"(6x,a,/3x,a,/a)")&
439                '</Piece>','</PolyData>','</VTKFile>'
440     
441           !Now write the facet info
442     
443           write(stl_unit,*)'solid vcg'
444     
445           write(stl_unit,*) '   facet normal ', NORM_FACE(1:3,FID)
446           write(stl_unit,*) '      outer loop'
447           write(stl_unit,*) '         vertex ', VERTEX(1,1:3,FID)
448           write(stl_unit,*) '         vertex ', VERTEX(2,1:3,FID)
449           write(stl_unit,*) '         vertex ', VERTEX(3,1:3,FID)
450           write(stl_unit,*) '      endloop'
451           write(stl_unit,*) '   endfacet'
452     
453           write(stl_unit,*)'endsolid vcg'
454     
455           close(vtp_unit, status = 'keep')
456           close(stl_unit, status = 'keep')
457           write(*,*) 'wrote a facet and a parcel. now waiting'
458           read(*,*)
459           end SUBROUTINE write_this_facet_and_parcel
460     
461