modification to user_mark & raw data saving mehtod

ProgramDefinedGeometry.f90

@@ -1,44 +1,10 @@
-!-------------- cemfDEM code: a dem simulation code ----------------------------
-! D C enter of
-! D M E ngineering and
-! D M M ultiscale modeling of
-! D M F luid flow 
-! EEEEEEEEM .ir
-!------------------------------------------------------------------------------
-! Copyright (C): cemf
-! website: www.cemf.ir
-!------------------------------------------------------------------------------ 
-! This file is part of cemfDEM code. It is a free software for simulating 
-! granular flow. You can redistribute it and/or modify it under the terms of
-! GNU General Public License version 3 or any other later versions. 
-! 
-! cemfDEM code is distributed to help others in their research in the field 
-! of granular flow, but WITHOUT ANY WARRANTY; without even the implied 
-! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-!------------------------------------------------------------------------------
 
 subroutine ProgramDefinedGeometry( Geom )
  use g_Geometry
  use g_Line
  implicit none
  class(Geometry),intent(in):: Geom
 
- !// locals 
- type(real3) p1, p2, p3, p4
- type(PlaneWall) plane
- type(CylinderWall) cyl,cyl1,cyl2
- logical res
-
- ! cylinder shell, the radius is 0.1 m 
- res = cyl%CreateCylinder( 0.1_RK, 0.1_RK, p_line( real3(0.0, 0.0, 0.0), real3(0.0,0.0, 0.03) ),24, 1, 1 )
- call Geom%add_Cylinder( cyl ) 
-
- ! rear wall 
- res = cyl1%CreateCylinder( 0.001_RK, 0.1_RK, p_line( real3(0.0, 0.0, -0.00001), real3(0.0,0.0, 0.0) ),24, 1, 1 )
- call Geom%add_Cylinder( cyl1 )
-
- ! front wall 
- res = cyl2%CreateCylinder( 0.1_RK, 0.001_RK, p_line( real3(0.0, 0.0, 0.03), real3(0.0,0.0, 0.03001) ),24, 1, 1 )
- call Geom%add_Cylinder( cyl2 )
+
 
 end subroutine 

User_Mark.f90

@@ -1,55 +1,17 @@
-!-------------- cemfDEM code: a dem simulation code ----------------------------
-! D C enter of
-! D M E ngineering and
-! D M M ultiscale modeling of
-! D M F luid flow 
-! EEEEEEEEM .ir
-!------------------------------------------------------------------------------
-! Copyright (C): cemf
-! website: www.cemf.ir
-!------------------------------------------------------------------------------ 
-! This file is part of cemfDEM code. It is a free software for simulating 
-! granular flow. You can redistribute it and/or modify it under the terms of
-! GNU General Public License version 3 or any other later versions. 
-! 
-! cemfDEM code is distributed to help others in their research in the field 
-! of granular flow, but WITHOUT ANY WARRANTY; without even the implied 
-! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-!------------------------------------------------------------------------------
-
-function User_Mark( id, flag, ptype, dpos ) result (mark)
+function User_Mark( id, flag, ptype, dpos, oldMark ) result (mark)
  use g_TypeDef
+ use g_Prtcl_DefaultValues
  implicit none 
 
- integer(IK),intent(in):: id, flag, ptype
+ integer(IK),intent(in):: id, flag, ptype, oldMark
  type(real4),intent(in):: dpos 
  integer(IK) mark
-
+
+
  !// locals
- real minx, maxx, dx
- integer(IK) numx
- real(RK) dist
- real(RK) d1, d2, d3
- type(real3) pos, point1, point2, point3
 
- ! position and diameter of regions which are marked as tracer
- pos = dpos
- d1 = 0.025
- d2 = 0.025
- d3 = 0.025
- point1 = real3(0.0, -0.018, 0.0)
- point2 = real3(0.02, -0.042, 0.0)
- point3 = real3(0.043, -0.066, 0.0)
+ mark = oldMark
+
+ return 
 
- ! finds particles in each region
- if( sqrt((pos%x- point1%x)**2 + (pos%y- point1%y)**2 ) .le. d1/2) then
- mark = 1
- elseif( sqrt((pos%x- point2%x)**2 + (pos%y- point2%y)**2 ) .le. d2/2) then
- mark = 2
- elseif( sqrt((pos%x- point3%x)**2 + (pos%y- point3%y)**2 ) .le. d3/2) then
- mark = 3
- else
- mark = 0
- end if
-
 end function

main.f90

@@ -1,24 +1,3 @@
-!-------------- cemfDEM code: a dem simulation code ----------------------------
-! D C enter of
-! D M E ngineering and
-! D M M ultiscale modeling of
-! D M F luid flow 
-! EEEEEEEEM .ir
-!------------------------------------------------------------------------------
-! Copyright (C): cemf
-! website: www.cemf.ir
-!------------------------------------------------------------------------------ 
-! This file is part of cemfDEM code. It is a free software for simulating 
-! granular flow. You can redistribute it and/or modify it under the terms of
-! GNU General Public License version 3 or any other later versions. 
-! 
-! cemfDEM code is distributed to help others in their research in the field 
-! of granular flow, but WITHOUT ANY WARRANTY; without even the implied 
-! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-!------------------------------------------------------------------------------
-
-program main
-
  use g_DEMSystem
  use g_stl_reader
  use g_Prtcl_PureProperty
@@ -42,14 +21,14 @@ program main
 !//// Initial settings 
 
  ! options and settings for DEM simulation
- DEM_opt%gravity = real3( 0.0, -9.8, 0.0 )
+ DEM_opt%gravity = real3( 0.0, 0.0, -9.8 )
  DEM_opt%SaveFreq = 1000
- DEM_opt%RunName = "mixingRotatingDrumRolling"
+ DEM_opt%RunName = "Continuous_Mixer"
  DEM_opt%Res_Dir = "./Results/"
- DEM_opt%OutputFileType = OP_Type_VTK
- DEM_opt%CT_model = CTM_ConstantTorque
- DEM_opt%CF_Type = CFT_nLin_nl
- DEM_opt%CS_Method = CSM_NBS_Munjiza !_Hrchl !
+ DEM_opt%OutputFileType = OP_Type_VTK ! vtk output
+ DEM_opt%CT_model = CTM_ConstantTorque ! constant torque model
+ DEM_opt%CF_Type = CFT_nLin_nl ! non-linear model
+ DEM_opt%CS_Method = CSM_NBS_Munjiza_Hrchl !for polydisperse system
  DEM_opt%PI_Method = PIM_AB4AM5
  DEM_opt%PRI_Method= PIM_AB4AM5
 
@@ -61,58 +40,69 @@ program main
  !//// particles
  allocate(PSD)
 
- ! 19000 particles with diameter of 3 mm
- PSD = PS_Distribution( 19000 , PSD_Uniform, 1 , 0.003_RK, 0.003_RK)
+ ! tri-modal distribution with size range between 3 to 6 mm 
+ PSD = PS_Distribution( 600000 , PSD_Uniform, 1 , 0.006_RK, 0.006_RK)
 
  ! associates particle size distribution with property
  allocate(PSDP)
  PSDP = PSD_Property( PSD , prop )
  deallocate(PSD)
 
- ! positions particles inside the drum
- allocate(Pos)
- Pos = PSDP_Position( PSDP )
- call Pos%PositionOrdered( real3(-0.071 , -0.071, 0.0), real3(0.071, 0.071, 0.03) )
-
+
 !// main components of DEM system 
  !//// geometry
  allocate(geom)
- call ProgramDefinedGeometry(geom)
+ ! reads geometry from stl files 
+ call geom%add_stl_file("./body.stl", 1, 1, .true. ) ! shell, user_id = 1
+ call geom%add_stl_file("./impeller.stl", 2, 1, .true. ) ! screw, user_id = 2
+ call geom%add_stl_file("./cylender.stl", 3, 1, .true. ) ! cylender, user_id = 3
+
+ ! separated output for geometry
+ ! user_id = 1: separate files with a name containing "body"
+ ! user_id = 2: separate files with a name containing "impeller"
+ call geom%setWallOutputName( (/1,2,3/) , (/"body","impe","cyle"/) )
 
- !//// Property
+ !//// Property 
+ !!!DEM_Options ,dynamic friction, rolling friction, coefficient of normal restitution, coefficient of tangential restitution 
  allocate( Property )
  call Property%ParticleProperty( PSDP ) ! particles 
  call Property%WallProperty(1 , (/prop/) ) ! walls 
- call Property%PP_BinaryProp(DEM_opt, 0.5_RK, 0.1_RK, 0.8_RK, 0.8_RK) ! binary pp
- call Property%PW_BinaryProp(DEM_opt, 0.7_RK, 0.1_RK, 0.8_RK, 0.8_RK) ! binary pw
+ call Property%PP_BinaryProp(DEM_opt, 0.3_RK, 0.1_RK, 0.7_RK, 0.7_RK) ! binary pp !!! inaro chi bezaram??
+ call Property%PW_BinaryProp(DEM_opt, 0.4_RK, 0.1_RK, 0.7_RK, 0.7_RK) ! binary pw
 
  !//// DEM system with particle insertion
  ! simulation domain
- minDomain = real3(-0.11, -0.11, -0.01)
- maxDomain = real3( 0.11, 0.11 , 0.04)
-
-
- ! initializes DEM system, dt = 0.00002 sec. 
- call DEM%Initialize( 0.00002_RK, Pos ,geom, Property, minDomain, maxDomain, DEM_opt ) 
-
-
-!//// iteration loop for 0.2 seconds. This lets particles settle under the gravity 
- call DEM%iterate(10000)
-
-
- !// sets the rotating velocity of drum, it is 10 RPM
- call geom%setWallVelocity(1, real3(0,0,0), &
- RPMtoRAD_S(10.0_RK), &
- p_line( real3(0,0,0) , real3(0,0,0.2) ) )
-
- !//// iteration for 0.8 more seconds to reach steady condition
- call DEM%iterate(39999)
- !// marks three groups of particles as tracers 
- call DEM%User_prtclMark()
+ minDomain = real3(-0.11 , -0.01, -0.13)
+ maxDomain = real3( 0.11, 0.65 , 0.12)
+
+ ! insertion plane near the inlet gate of conveyor
+ p1 = real3( -0.050 , 0.065, 0.101)
+ p2 = real3( 0.050, 0.065,0.101)
+ p3 = real3( 0.050 , 0.010,0.101)
+ p4 = real3( -0.050 , 0.010, 0.101)
+ res = insPlane%CreateWall_nv(p1,p2,p3,p4)
+ !!!----------------------------------------------------------------------------------------------------------------------------
+ ! initializes DEM system, dt = 0.00001 s, particles should be inserted during 2 !!!in chie?
+ ! million iterations (equivalent to 20 seconds)
+ call DEM%Initialize( 0.00001_RK, PSDP, InsPlane, 6000000, 0.3_RK ,geom, Property, minDomain, maxDomain, DEM_opt )
+ !initial velocity m/s
+
+!//// iteration loop for 0.1 seconds 
  call DEM%iterate(1)
-
-
- !/// iteration for 20 more seconds 
- call DEM%iterate(1000000)
+
+ ! sets the velocity of the impeller (user_id = 2)
+ ! rotation speed is 60 RPM and rotation axis is the central shaft 
+ call geom%setWallVelocity( 2, real3(0,0,0), &
+ RPMtoRAD_S(-120.0_RK), &
+ p_line( real3(0,0,0) , real3(0,0.75,0.0) ) )
+
+!/// iteration for 20 seconds
+ call DEM%User_prtclMark(); 
+
+ do i = 1,20000000
+ call DEM%User_prtclMark(); 
+ call DEM%iterate(10)
+ end do
+
 
 end program 

src/g_DEMSystem.f90

@@ -58,9 +58,9 @@ subroutine ProgramDefinedGeometry( Geom )
  end interface
 
  interface 
- function User_Mark( id, flag, ptype, dpos ) result (mark)
+ function User_Mark( id, flag, ptype, dpos, oldMark ) result (mark)
  import IK, real4
- integer(IK),intent(in):: id, flag, ptype
+ integer(IK),intent(in):: id, flag, ptype, oldMark
  type(real4),intent(in):: dpos 
  integer(IK) mark
  end function
@@ -163,6 +163,8 @@ function User_Mark( id, flag, ptype, dpos ) result (mark)
  procedure:: write_prtcl => DEMS_write_prtcl
  procedure:: write_geom => DEMS_write_geom
 
+ procedure:: writeRawData => DEMS_WriteRawData
+
  ! VTK file output for particles
  procedure:: PrtclToVTK => DEMS_PrtclToVTK
 
@@ -690,6 +692,49 @@ subroutine DEMS_write_geom( this, num )
 end subroutine 
 
 
+!**********************************************************************
+! writing rawData bed condition of particles in an external file
+!**********************************************************************
+subroutine DEMS_WriteRawData(this, num)
+ implicit none
+ class(DEMSystem) this
+ integer(IK),intent(in):: num
+
+ !// locals 
+ integer(IK):: nUnit = d_wrt_prtc_unit
+ integer(IK) i, num_InDomain
+ character(256) chFile, ch
+
+ !// body
+
+ ! openning file
+ write(ch,*) this%m_DEM_Options%base_Number + num
+ chFile = trim(this%m_DEM_Options%Res_Dir)//'RawData'//trim(adjustl(ch))//"for"//trim(this%m_DEM_Options%RunName)//".dat"
+ open( file = chFile, unit = nUnit )
+
+ num_InDomain = this%numInDomain()
+
+ ! first writing number of particles
+ write( nUnit , * ) "numberOfParticles" , num_inDomain 
+ write( nUnit , * ) "id type posision velocity mark"
+ ! writing particle data 
+ do i=1,this%numPrtcl
+
+ if( this%prtcl_flag(i) >= Pflg_inDomain ) then
+
+ write( nUnit , * ) this%prtcl_ids(i) , &
+ this%prtcl_type(i), &
+ this%prtcl_dpos(i), &
+ this%prtcl_vel(i) , &
+ this%prtcl_usr_mark(i)
+ end if
+
+ end do
+
+ close( nUnit ) 
+
+end subroutine 
+
 !**********************************************************************
 ! writing particle data to a VTK file
 !**********************************************************************
@@ -817,6 +862,10 @@ subroutine DEMS_outputFile(this )
  !call this%write_prtcl(this%iterNumber)
  print*, "*********** Binary file is not specified *****************"
  end if
+
+ if( iAND(this%m_DEM_options%OutputFileType,OP_Type_Bin) == OP_Type_Raw ) then
+ call this%writeRawData(this%iterNumber)
+ end if
 
  endif
 
@@ -1172,7 +1221,12 @@ subroutine DEMS_User_prtclMark( this )
 
  do i=1, this%max_nPrtcl
 
- this%prtcl_usr_mark(i) = User_Mark( this%prtcl_ids(i), this%prtcl_flag(i), this%prtcl_type(i), this%prtcl_dpos(i))
+ this%prtcl_usr_mark(i) = User_Mark &
+ ( this%prtcl_ids(i), &
+ this%prtcl_flag(i), &
+ this%prtcl_type(i), &
+ this%prtcl_dpos(i), &
+ this%prtcl_usr_mark(i))
 
  end do
 
@@ -1378,6 +1432,8 @@ subroutine DEMS_report_contacts( this )
 end subroutine
 
 
+
+
 end module
 
 

src/g_Prtcl_DefaultValues.f90

@@ -105,6 +105,7 @@ module g_Prtcl_DefaultValues
  integer(IK),parameter:: OP_Type_VTK = 1_IK
  integer(IK),parameter:: OP_Type_Tec = 2_IK
  integer(IK),parameter:: OP_Type_Bin = 4_IK
+ integer(IK),parameter:: OP_Type_Raw = 8_IK
 
 
  ! default values