-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathin.flow.couette.3d.u2_bon
340 lines (224 loc) · 11.4 KB
/
in.flow.couette.3d.u2_bon
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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
################################################################################
# 3-d LJ flow simulation
# Author: D. Chibouti
################################################################################
#-------------------------------------------------------------------------------
# settings variables
#-------------------------------------------------------------------------------
variable x equal 10
variable y equal 10
variable z equal 1
variable dx equal floor(${x}/9)
variable dy equal floor(${y}/9)
variable dz equal floor(${z}/9)
variable rho equal 0.6 # fluid density
variable RO equal 0.8 # wall density
variable T equal 1.5 # temperature
variable Ti equal 8.0 # wall temperature !!!
variable V equal vol # volume
variable v0 equal 2.0 # velocity
variable rc equal 2.5 # cutoff radius for shifted LJ-potential
variable dt equal 0.006 # timestep delta t (pas de temps)
variable p equal 100 # correlation length (step)
variable s equal 5 # sample interval (echantillon)
variable d equal $p*$s # dump interval (timesteps interval)
variable tprod equal 5000 # simulation time for the production run
variable Nprod equal floor(${tprod}/${dt}) # total number of timesteps
#-------------------------------------------------------------------------------
# settings
#-------------------------------------------------------------------------------
# options used for fix ave/time; sample the quantities every 10 steps
variable Niter0 equal 3000
variable Niter1 equal 5000
variable Niter equal 10000
variable Nsamp equal 10
variable Nrepeat equal floor(${Nprod}/${Nsamp})
variable Nevery equal ${Nsamp}*${Nrepeat}
#-------------------------------------------------------------------------------
# convert from LAMMPS real units to SI
#-------------------------------------------------------------------------------
variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable atm2Pa equal 101325.0
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal ${atm2Pa}*${atm2Pa}*${fs2s}*${A2m}*${A2m}*${A2m}
#-------------------------------------------------------------------------------
# setup problem
#-------------------------------------------------------------------------------
units lj # define units
dimension 3 # dimension : 2d, 3d
boundary p p p # specify periodic boundary conditions
atom_style atomic
neighbor 0.3 bin # specify parameters for neighbor list
# rnbr = rcut + 0.3
neigh_modify delay 5
#-------------------------------------------------------------------------------
#
# create geometry : regions, groups ...
#
#-------------------------------------------------------------------------------
lattice fcc ${RO} # fcc = face centered cubic (atomic structure)
# density : {RO} for solid
#*******************************************************************************
# define regions
#*******************************************************************************
region simbox block 0.0 15.0 0.0 15.0 0.0 30.0
region upper_wall block 0.0 14.5 0.0 14.5 20.0 25.0
region lower_wall block 0.0 14.5 0.0 14.5 5.0 10.0
region upper_space block 0.0 14.5 0.0 14.5 25.0 30.0
region lower_space block 0.0 14.5 0.0 14.5 0.0 5.0
region bulk block 1.0 13.0 1.0 13.0 12.0 18.0
create_box 2 simbox
# create atoms for regions "upper_wall" & "lower_wall"
create_atoms 1 region upper_wall
create_atoms 1 region lower_wall
mass 1 1.0
#-------------------------------------------------------------------------------
# define groups
#-------------------------------------------------------------------------------
# create the group "wall"
group upper region upper_wall
group lower region lower_wall
group wall union upper lower
set group wall type 1
# create the group "wall_space"
group upper_space region upper_space
group lower_space region lower_space
group wall_space union upper_space lower_space
set group wall_space type 1
lattice fcc ${rho} # fcc = face centered cubic (atomic structure)
# density : {rho} for fluid
create_atoms 2 region bulk
mass 2 1.0
# create the group "fluid"
group fluid subtract all wall wall_space
set group fluid type 2
# create the group "l_wall_no_space" and "u_wall_no_space"
group l_wall_no_space subtract all wall_space upper fluid
group u_wall_no_space subtract all wall_space lower fluid
#-------------------------------------------------------------------------------
#
# Interaction potential
#
#-------------------------------------------------------------------------------
pair_style lj/cut ${rc}
pair_coeff 1 1 50.0 1.0 ${rc} # ... between atoms of type 1 and 1
pair_coeff 2 2 1.0 1.0 ${rc} # ... between atoms of type 2 and 2
pair_coeff 1 2 1.0 1.0 ${rc} # ... between atoms of type 1 and 2
pair_modify shift yes
#-------------------------------------------------------------------------------
#
# SIMULATION
#
#-------------------------------------------------------------------------------
fix force wall_space setforce 0.0 0.0 0.0
#-------------------------------------------------------------------------------
# initial velocities
#-------------------------------------------------------------------------------
compute mobile fluid temp/profile 1 0 0 z 20
velocity fluid create $T 482748 temp mobile
fix 1 all nve
fix 2 fluid temp/rescale 200 $T $T 0.02 1.0
fix_modify 2 temp mobile
#-------------------------------------------------------------------------------
#
# SIMULATION - 1st run computation
#
#-------------------------------------------------------------------------------
# langevin thermostat :
fix lang_upper u_wall_no_space langevin ${T} ${T} 0.5 6586
fix lang_lower l_wall_no_space langevin $T $T 0.5 3215
# compute temp : # compute mobile fluid temp
compute temp_fluid fluid temp
compute temp_lower lower temp
compute temp_upper upper temp
# set variables
variable T_fluid equal c_temp_fluid
variable T_lower equal c_temp_lower
variable T_upper equal c_temp_upper
# setting run parameters :
timestep ${dt}
thermo_style custom step temp c_temp_fluid c_temp_lower c_temp_upper pe ke etotal
thermo ${d}
run ${Niter0}
#*******************************************************************************
# save process !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#*******************************************************************************
write_restart restart.save
#reset_timestep 0
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# outputs
#-------------------------------------------------------------------------------
# compute per-atom potencial and kinetic energy
compute ke all ke/atom
compute pe all pe/atom
variable Temp1 atom c_ke/1.5 # temperature !!!
variable NiterA equal ${Niter0}+${Niter1}
print "***************************************************************"
print " run after equilibration "
print "***************************************************************"
# run simulation
run ${Niter1}
#-------------------------------------------------------------------------------
# couette flow inialisation
#-------------------------------------------------------------------------------
# initial velocities : # velocity fluid create $T 482748 temp mobile
velocity fluid create $T 482748 temp mobile
fix 1 all nve
fix 2 fluid temp/rescale 200 $T $T 0.02 1.0
fix_modify 2 temp mobile
#velocity fluid create $T 12345
fix lang_upper u_wall_no_space langevin ${T} ${T} 0.5 6586
fix lang_lower l_wall_no_space langevin $T $T 0.5 3215
#-------------------------------------------------------------------------------
#
# Couette flow
#
#-------------------------------------------------------------------------------
# velocities
# condition 1 : magnifique temperature profile {but not velocity :( }
#velocity l_wall_no_space set 0.0 0.0 0.0 units box # stationary wall
#velocity u_wall_no_space set ${v0} 0.0 0.0 units box # mobile wall
# condition 2 : magnifique velocity profile {but not temperature :( }
velocity lower set 0.0 0.0 0.0 units box
velocity upper set ${v0} 0.0 0.0 units box
# condition 3 : !! can't fix it !!!!!!!
variable x0 equal vdisplace(0.0,${v0})
#fix d1 upper move variable v_x0 NULL NULL v_v0 NULL NULL # mobile wall
# rescaling temp !!!
#fix 20 upper temp/rescale 200 $T $T 0.02 1.0
#*******************************************************************************
# Maxwellian reflection ??????????
# upper wall velocity = 2.0, lower wall fixed,
# accommodation coefficient = 0.2 for both walls
#fix zwalls all wall/reflect/stochastic maxwell 29839 &
# zlo EDGE $T 0 0 0 0.2 zhi EDGE $T 2.0 0 0 0.2
# Deffusive reflection ??????????
# upper wall velocity = 2.0, lower wall fixed
#fix zwalls all wall/reflect/stochastic diffusive 23424 &
zlo EDGE 1.5 0 0 0 zhi EDGE 1.5 2.0 0 0
#*******************************************************************************
#fix 3 wall setforce 0.0 0.0 0.0
#fix 4 all addforce 0.037 0.0 0.0 # for poiseuille flow !!!
#------------------------------------------------------------------------------
# save to file
#------------------------------------------------------------------------------
variable NiterB equal ${Niter0}+${Niter1}+${Niter}
compute layer2 all chunk/atom bin/1d z lower 0.0125 units reduced # all
compute layer3 fluid chunk/atom bin/1d z lower 0.0125 units reduced # fluid
fix fsave2 all ave/chunk 1 200 ${NiterB} layer2 vx vy vz v_Temp1 density/mass &
norm sample file profile_vel_Temp.all
fix fsave22 all ave/chunk 1 200 ${NiterB} layer2 vx vy vz temp density/mass &
norm sample file profile_vel_Temp.all_temp
fix fsave3 fluid ave/chunk 1 200 ${NiterB} layer3 vx vy vz v_Temp1 density/mass &
norm sample file profile_vel_Temp.fluid
fix fsave33 fluid ave/chunk 1 200 ${NiterB} layer3 vx vy vz temp density/mass &
norm sample file profile_vel_Temp.fluid_temp
#-------------------------------------------------------------------------------
print "***************************************************************"
print " final run "
print "***************************************************************"
# run simulation
run ${Niter}
#-------------------------------------------------------------------------------