← Back to branch summary

~fluidity-core/fluidity/sea-ice-coupling

« back to all changes in this revision

Viewing changes to tests/mphase_wen_yu_drag_correlation/mphase_wen_yu_drag_correlation.flml

  • Committer: Simon Mouradian
  • Date: 2012-10-09 20:00:15 UTC
  • Revision ID: mouradian@gmail.com-20121009200015-7hztzrng2bkqgz1t
revert last trunk merge, something broke

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/mphase_wen_yu_drag_correlation/mphase_wen_yu_drag_correlation.flml
1
 
<?xml version='1.0' encoding='utf-8'?>
2
 
<fluidity_options>
3
 
  <simulation_name>
4
 
    <string_value lines="1">mphase_wen_yu_drag_correlation</string_value>
5
 
  </simulation_name>
6
 
  <problem_type>
7
 
    <string_value lines="1">multiphase</string_value>
8
 
  </problem_type>
9
 
  <geometry>
10
 
    <dimension>
11
 
      <integer_value rank="0">2</integer_value>
12
 
    </dimension>
13
 
    <mesh name="CoordinateMesh">
14
 
      <from_file file_name="mphase_wen_yu_drag_correlation">
15
 
        <format name="triangle"/>
16
 
        <stat>
17
 
          <include_in_stat/>
18
 
        </stat>
19
 
      </from_file>
20
 
    </mesh>
21
 
    <mesh name="VelocityMesh">
22
 
      <from_mesh>
23
 
        <mesh name="CoordinateMesh"/>
24
 
        <stat>
25
 
          <exclude_from_stat/>
26
 
        </stat>
27
 
      </from_mesh>
28
 
    </mesh>
29
 
    <mesh name="PressureMesh">
30
 
      <from_mesh>
31
 
        <mesh name="CoordinateMesh"/>
32
 
        <stat>
33
 
          <exclude_from_stat/>
34
 
        </stat>
35
 
      </from_mesh>
36
 
    </mesh>
37
 
    <quadrature>
38
 
      <degree>
39
 
        <integer_value rank="0">3</integer_value>
40
 
      </degree>
41
 
    </quadrature>
42
 
  </geometry>
43
 
  <io>
44
 
    <dump_format>
45
 
      <string_value>vtk</string_value>
46
 
    </dump_format>
47
 
    <dump_period>
48
 
      <constant>
49
 
        <real_value rank="0">0.0</real_value>
50
 
      </constant>
51
 
    </dump_period>
52
 
    <output_mesh name="VelocityMesh"/>
53
 
    <stat/>
54
 
  </io>
55
 
  <timestepping>
56
 
    <current_time>
57
 
      <real_value rank="0">0</real_value>
58
 
    </current_time>
59
 
    <timestep>
60
 
      <real_value rank="0">0.001</real_value>
61
 
    </timestep>
62
 
    <finish_time>
63
 
      <real_value rank="0">5.0</real_value>
64
 
    </finish_time>
65
 
    <steady_state>
66
 
      <tolerance>
67
 
        <real_value rank="0">1.0e-7</real_value>
68
 
        <infinity_norm/>
69
 
      </tolerance>
70
 
      <steady_state_file>
71
 
        <binary_output/>
72
 
      </steady_state_file>
73
 
    </steady_state>
74
 
  </timestepping>
75
 
  <physical_parameters>
76
 
    <gravity>
77
 
      <magnitude>
78
 
        <real_value rank="0">9.8</real_value>
79
 
      </magnitude>
80
 
      <vector_field name="GravityDirection" rank="1">
81
 
        <prescribed>
82
 
          <mesh name="CoordinateMesh"/>
83
 
          <value name="WholeMesh">
84
 
            <constant>
85
 
              <real_value shape="2" dim1="dim" rank="1">0.0 -1.0</real_value>
86
 
            </constant>
87
 
          </value>
88
 
          <output/>
89
 
          <stat>
90
 
            <include_in_stat/>
91
 
          </stat>
92
 
          <detectors>
93
 
            <exclude_from_detectors/>
94
 
          </detectors>
95
 
        </prescribed>
96
 
      </vector_field>
97
 
    </gravity>
98
 
  </physical_parameters>
99
 
  <material_phase name="Water">
100
 
    <equation_of_state>
101
 
      <fluids>
102
 
        <linear>
103
 
          <reference_density>
104
 
            <real_value rank="0">1000.0</real_value>
105
 
          </reference_density>
106
 
        </linear>
107
 
      </fluids>
108
 
    </equation_of_state>
109
 
    <scalar_field name="Pressure" rank="0">
110
 
      <aliased material_phase_name="Tephra" field_name="Pressure"/>
111
 
    </scalar_field>
112
 
    <scalar_field name="Density" rank="0">
113
 
      <diagnostic>
114
 
        <algorithm name="Internal" material_phase_support="multiple"/>
115
 
        <mesh name="VelocityMesh"/>
116
 
        <output/>
117
 
        <stat/>
118
 
        <convergence>
119
 
          <include_in_convergence/>
120
 
        </convergence>
121
 
        <detectors>
122
 
          <include_in_detectors/>
123
 
        </detectors>
124
 
        <steady_state>
125
 
          <include_in_steady_state/>
126
 
        </steady_state>
127
 
      </diagnostic>
128
 
    </scalar_field>
129
 
    <vector_field name="Velocity" rank="1">
130
 
      <prognostic>
131
 
        <mesh name="VelocityMesh"/>
132
 
        <equation name="LinearMomentum"/>
133
 
        <spatial_discretisation>
134
 
          <continuous_galerkin>
135
 
            <stabilisation>
136
 
              <streamline_upwind>
137
 
                <nu_bar_optimal/>
138
 
                <nu_scale name="0.5">
139
 
                  <real_value shape="1" rank="0">0.5</real_value>
140
 
                </nu_scale>
141
 
              </streamline_upwind>
142
 
            </stabilisation>
143
 
            <mass_terms>
144
 
              <lump_mass_matrix/>
145
 
            </mass_terms>
146
 
            <advection_terms/>
147
 
            <stress_terms>
148
 
              <tensor_form/>
149
 
            </stress_terms>
150
 
          </continuous_galerkin>
151
 
          <conservative_advection>
152
 
            <real_value rank="0">0.0</real_value>
153
 
          </conservative_advection>
154
 
        </spatial_discretisation>
155
 
        <temporal_discretisation>
156
 
          <theta>
157
 
            <real_value rank="0">1.0</real_value>
158
 
          </theta>
159
 
          <relaxation>
160
 
            <real_value rank="0">0.5</real_value>
161
 
          </relaxation>
162
 
        </temporal_discretisation>
163
 
        <solver>
164
 
          <iterative_method name="gmres">
165
 
            <restart>
166
 
              <integer_value rank="0">30</integer_value>
167
 
            </restart>
168
 
          </iterative_method>
169
 
          <preconditioner name="sor"/>
170
 
          <relative_error>
171
 
            <real_value rank="0">1.0e-6</real_value>
172
 
          </relative_error>
173
 
          <max_iterations>
174
 
            <integer_value rank="0">1000</integer_value>
175
 
          </max_iterations>
176
 
          <never_ignore_solver_failures/>
177
 
          <diagnostics>
178
 
            <monitors/>
179
 
          </diagnostics>
180
 
        </solver>
181
 
        <initial_condition name="WholeMesh">
182
 
          <constant>
183
 
            <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value>
184
 
          </constant>
185
 
        </initial_condition>
186
 
        <boundary_conditions name="Sides">
187
 
          <surface_ids>
188
 
            <integer_value shape="1" rank="1">666</integer_value>
189
 
          </surface_ids>
190
 
          <type name="no_normal_flow"/>
191
 
        </boundary_conditions>
192
 
        <boundary_conditions name="Bottom">
193
 
          <surface_ids>
194
 
            <integer_value shape="1" rank="1">444</integer_value>
195
 
          </surface_ids>
196
 
          <type name="no_normal_flow"/>
197
 
        </boundary_conditions>
198
 
        <boundary_conditions name="Top">
199
 
          <surface_ids>
200
 
            <integer_value shape="1" rank="1">333</integer_value>
201
 
          </surface_ids>
202
 
          <type name="no_normal_flow"/>
203
 
        </boundary_conditions>
204
 
        <tensor_field name="Viscosity" rank="2">
205
 
          <prescribed>
206
 
            <value name="WholeMesh">
207
 
              <isotropic>
208
 
                <constant>
209
 
                  <real_value rank="0">0.001</real_value>
210
 
                </constant>
211
 
              </isotropic>
212
 
            </value>
213
 
            <output/>
214
 
          </prescribed>
215
 
        </tensor_field>
216
 
        <output/>
217
 
        <stat>
218
 
          <include_in_stat/>
219
 
          <previous_time_step>
220
 
            <exclude_from_stat/>
221
 
          </previous_time_step>
222
 
          <nonlinear_field>
223
 
            <exclude_from_stat/>
224
 
          </nonlinear_field>
225
 
        </stat>
226
 
        <convergence>
227
 
          <include_in_convergence/>
228
 
        </convergence>
229
 
        <detectors>
230
 
          <include_in_detectors/>
231
 
        </detectors>
232
 
        <steady_state>
233
 
          <include_in_steady_state/>
234
 
        </steady_state>
235
 
        <consistent_interpolation/>
236
 
      </prognostic>
237
 
    </vector_field>
238
 
    <scalar_field name="PhaseVolumeFraction" rank="0">
239
 
      <diagnostic>
240
 
        <mesh name="VelocityMesh"/>
241
 
        <algorithm name="Internal" material_phase_support="multiple"/>
242
 
        <output/>
243
 
        <stat/>
244
 
        <detectors>
245
 
          <exclude_from_detectors/>
246
 
        </detectors>
247
 
        <cap_values>
248
 
          <upper_cap>
249
 
            <real_value rank="0">1.0</real_value>
250
 
          </upper_cap>
251
 
          <lower_cap>
252
 
            <real_value rank="0">0.0</real_value>
253
 
          </lower_cap>
254
 
        </cap_values>
255
 
      </diagnostic>
256
 
    </scalar_field>
257
 
  </material_phase>
258
 
  <material_phase name="Tephra">
259
 
    <equation_of_state>
260
 
      <fluids>
261
 
        <linear>
262
 
          <reference_density>
263
 
            <real_value rank="0">2340.0</real_value>
264
 
          </reference_density>
265
 
        </linear>
266
 
      </fluids>
267
 
    </equation_of_state>
268
 
    <scalar_field name="Pressure" rank="0">
269
 
      <prognostic>
270
 
        <mesh name="PressureMesh"/>
271
 
        <spatial_discretisation>
272
 
          <continuous_galerkin>
273
 
            <remove_stabilisation_term/>
274
 
            <integrate_continuity_by_parts/>
275
 
          </continuous_galerkin>
276
 
        </spatial_discretisation>
277
 
        <reference_node>
278
 
          <integer_value rank="0">1</integer_value>
279
 
        </reference_node>
280
 
        <scheme>
281
 
          <poisson_pressure_solution>
282
 
            <string_value lines="1">only first timestep</string_value>
283
 
          </poisson_pressure_solution>
284
 
          <use_projection_method/>
285
 
        </scheme>
286
 
        <solver>
287
 
          <iterative_method name="cg"/>
288
 
          <preconditioner name="mg"/>
289
 
          <relative_error>
290
 
            <real_value rank="0">1.0e-6</real_value>
291
 
          </relative_error>
292
 
          <max_iterations>
293
 
            <integer_value rank="0">1000</integer_value>
294
 
          </max_iterations>
295
 
          <never_ignore_solver_failures/>
296
 
          <diagnostics>
297
 
            <monitors/>
298
 
          </diagnostics>
299
 
        </solver>
300
 
        <output/>
301
 
        <stat/>
302
 
        <convergence>
303
 
          <include_in_convergence/>
304
 
        </convergence>
305
 
        <detectors>
306
 
          <exclude_from_detectors/>
307
 
        </detectors>
308
 
        <steady_state>
309
 
          <exclude_from_steady_state/>
310
 
        </steady_state>
311
 
        <no_interpolation/>
312
 
      </prognostic>
313
 
    </scalar_field>
314
 
    <scalar_field name="Density" rank="0">
315
 
      <diagnostic>
316
 
        <algorithm name="Internal" material_phase_support="multiple"/>
317
 
        <mesh name="VelocityMesh"/>
318
 
        <output/>
319
 
        <stat/>
320
 
        <convergence>
321
 
          <include_in_convergence/>
322
 
        </convergence>
323
 
        <detectors>
324
 
          <include_in_detectors/>
325
 
        </detectors>
326
 
        <steady_state>
327
 
          <include_in_steady_state/>
328
 
        </steady_state>
329
 
      </diagnostic>
330
 
    </scalar_field>
331
 
    <vector_field name="Velocity" rank="1">
332
 
      <prognostic>
333
 
        <mesh name="VelocityMesh"/>
334
 
        <equation name="LinearMomentum"/>
335
 
        <spatial_discretisation>
336
 
          <continuous_galerkin>
337
 
            <stabilisation>
338
 
              <streamline_upwind>
339
 
                <nu_bar_optimal/>
340
 
                <nu_scale name="0.5">
341
 
                  <real_value shape="1" rank="0">0.5</real_value>
342
 
                </nu_scale>
343
 
              </streamline_upwind>
344
 
            </stabilisation>
345
 
            <mass_terms>
346
 
              <lump_mass_matrix/>
347
 
            </mass_terms>
348
 
            <advection_terms/>
349
 
            <stress_terms>
350
 
              <tensor_form/>
351
 
            </stress_terms>
352
 
          </continuous_galerkin>
353
 
          <conservative_advection>
354
 
            <real_value rank="0">0.0</real_value>
355
 
          </conservative_advection>
356
 
        </spatial_discretisation>
357
 
        <temporal_discretisation>
358
 
          <theta>
359
 
            <real_value rank="0">1.0</real_value>
360
 
          </theta>
361
 
          <relaxation>
362
 
            <real_value rank="0">0.5</real_value>
363
 
          </relaxation>
364
 
        </temporal_discretisation>
365
 
        <solver>
366
 
          <iterative_method name="gmres">
367
 
            <restart>
368
 
              <integer_value rank="0">30</integer_value>
369
 
            </restart>
370
 
          </iterative_method>
371
 
          <preconditioner name="sor"/>
372
 
          <relative_error>
373
 
            <real_value rank="0">1.0e-6</real_value>
374
 
          </relative_error>
375
 
          <max_iterations>
376
 
            <integer_value rank="0">1000</integer_value>
377
 
          </max_iterations>
378
 
          <never_ignore_solver_failures/>
379
 
          <diagnostics>
380
 
            <monitors/>
381
 
          </diagnostics>
382
 
        </solver>
383
 
        <initial_condition name="WholeMesh">
384
 
          <constant>
385
 
            <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value>
386
 
          </constant>
387
 
        </initial_condition>
388
 
        <boundary_conditions name="Sides">
389
 
          <surface_ids>
390
 
            <integer_value shape="1" rank="1">666</integer_value>
391
 
          </surface_ids>
392
 
          <type name="no_normal_flow"/>
393
 
        </boundary_conditions>
394
 
        <boundary_conditions name="Bottom">
395
 
          <surface_ids>
396
 
            <integer_value shape="1" rank="1">444</integer_value>
397
 
          </surface_ids>
398
 
          <type name="no_normal_flow"/>
399
 
        </boundary_conditions>
400
 
        <boundary_conditions name="Top">
401
 
          <surface_ids>
402
 
            <integer_value shape="1" rank="1">333</integer_value>
403
 
          </surface_ids>
404
 
          <type name="no_normal_flow"/>
405
 
        </boundary_conditions>
406
 
        <tensor_field name="Viscosity" rank="2">
407
 
          <prescribed>
408
 
            <value name="WholeMesh">
409
 
              <isotropic>
410
 
                <constant>
411
 
                  <real_value rank="0">0.001</real_value>
412
 
                </constant>
413
 
              </isotropic>
414
 
            </value>
415
 
            <output/>
416
 
          </prescribed>
417
 
        </tensor_field>
418
 
        <output/>
419
 
        <stat>
420
 
          <include_in_stat/>
421
 
          <previous_time_step>
422
 
            <exclude_from_stat/>
423
 
          </previous_time_step>
424
 
          <nonlinear_field>
425
 
            <exclude_from_stat/>
426
 
          </nonlinear_field>
427
 
        </stat>
428
 
        <convergence>
429
 
          <include_in_convergence/>
430
 
        </convergence>
431
 
        <detectors>
432
 
          <include_in_detectors/>
433
 
        </detectors>
434
 
        <steady_state>
435
 
          <include_in_steady_state/>
436
 
        </steady_state>
437
 
        <consistent_interpolation/>
438
 
      </prognostic>
439
 
    </vector_field>
440
 
    <scalar_field name="PhaseVolumeFraction" rank="0">
441
 
      <prognostic>
442
 
        <mesh name="VelocityMesh"/>
443
 
        <equation name="AdvectionDiffusion"/>
444
 
        <spatial_discretisation>
445
 
          <control_volumes>
446
 
            <face_value name="FirstOrderUpwind"/>
447
 
            <diffusion_scheme name="ElementGradient"/>
448
 
          </control_volumes>
449
 
          <conservative_advection>
450
 
            <real_value rank="0">1.0</real_value>
451
 
          </conservative_advection>
452
 
        </spatial_discretisation>
453
 
        <temporal_discretisation>
454
 
          <theta>
455
 
            <real_value rank="0">1.0</real_value>
456
 
          </theta>
457
 
        </temporal_discretisation>
458
 
        <solver>
459
 
          <iterative_method name="gmres">
460
 
            <restart>
461
 
              <integer_value rank="0">30</integer_value>
462
 
            </restart>
463
 
          </iterative_method>
464
 
          <preconditioner name="sor"/>
465
 
          <relative_error>
466
 
            <real_value rank="0">1.0e-7</real_value>
467
 
          </relative_error>
468
 
          <max_iterations>
469
 
            <integer_value rank="0">1000</integer_value>
470
 
          </max_iterations>
471
 
          <never_ignore_solver_failures/>
472
 
          <diagnostics>
473
 
            <monitors/>
474
 
          </diagnostics>
475
 
        </solver>
476
 
        <initial_condition name="WholeMesh">
477
 
          <constant>
478
 
            <real_value rank="0">0.05</real_value>
479
 
          </constant>
480
 
        </initial_condition>
481
 
        <output/>
482
 
        <stat/>
483
 
        <convergence>
484
 
          <include_in_convergence/>
485
 
        </convergence>
486
 
        <detectors>
487
 
          <include_in_detectors/>
488
 
        </detectors>
489
 
        <steady_state>
490
 
          <include_in_steady_state/>
491
 
        </steady_state>
492
 
        <consistent_interpolation/>
493
 
      </prognostic>
494
 
    </scalar_field>
495
 
    <scalar_field name="ParticleReynoldsNumber" rank="0">
496
 
      <diagnostic>
497
 
        <algorithm name="particle_reynolds_number" material_phase_support="multiple">
498
 
          <depends>
499
 
            <string_value lines="1">Water::Velocity,Tephra::Velocity,Water::Density</string_value>
500
 
          </depends>
501
 
          <particle_diameter>
502
 
            <real_value rank="0">1e-3</real_value>
503
 
          </particle_diameter>
504
 
          <continuous_phase_name>Water</continuous_phase_name>
505
 
        </algorithm>
506
 
        <mesh name="VelocityMesh"/>
507
 
        <output/>
508
 
        <stat/>
509
 
        <convergence>
510
 
          <include_in_convergence/>
511
 
        </convergence>
512
 
        <detectors>
513
 
          <include_in_detectors/>
514
 
        </detectors>
515
 
        <steady_state>
516
 
          <include_in_steady_state/>
517
 
        </steady_state>
518
 
      </diagnostic>
519
 
    </scalar_field>
520
 
    <scalar_field name="DragCoefficient" rank="0">
521
 
      <diagnostic>
522
 
        <algorithm name="scalar_python_diagnostic" material_phase_support="multiple">
523
 
          <string_value type="code" lines="20" language="python">from math import sqrt, pi
524
 
vfrac_tephra = states["Tephra"].scalar_fields["PhaseVolumeFraction"]
525
 
vfrac_water = states["Water"].scalar_fields["PhaseVolumeFraction"]
526
 
 
527
 
velocity_tephra = states["Tephra"].vector_fields["Velocity"]
528
 
velocity_water = states["Water"].vector_fields["Velocity"]
529
 
 
530
 
density_tephra = states["Tephra"].scalar_fields["Density"]
531
 
density_water = states["Water"].scalar_fields["Density"]
532
 
 
533
 
viscosity_water = states["Water"].tensor_fields["Viscosity"]
534
 
 
535
 
# Particle diameter
536
 
d = 1e-3
537
 
 
538
 
for i in range(field.node_count):
539
 
  # Reynolds number
540
 
  Re = (vfrac_water.node_val(i)*density_water.node_val(i)*sqrt((velocity_water.node_val(i)[0] - velocity_tephra.node_val(i)[0])**2 + (velocity_water.node_val(i)[1] - velocity_tephra.node_val(i)[1])**2)*d)/viscosity_water.node_val(0)[0][0]
541
 
  
542
 
  # Drag coefficient
543
 
  if(Re &lt; 1000):
544
 
    C = (24.0/Re)*(1.0 + 0.15*Re**0.687)
545
 
  else:
546
 
    C = 0.44
547
 
    
548
 
  # Don't let C be NaN
549
 
  if(Re == 0):
550
 
    C = 1e16
551
 
 
552
 
  field.set(i, C)</string_value>
553
 
        </algorithm>
554
 
        <mesh name="VelocityMesh"/>
555
 
        <output/>
556
 
        <stat/>
557
 
        <convergence>
558
 
          <include_in_convergence/>
559
 
        </convergence>
560
 
        <detectors>
561
 
          <include_in_detectors/>
562
 
        </detectors>
563
 
        <steady_state>
564
 
          <include_in_steady_state/>
565
 
        </steady_state>
566
 
      </diagnostic>
567
 
    </scalar_field>
568
 
    <vector_field name="DragForce" rank="1">
569
 
      <diagnostic>
570
 
        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
571
 
          <string_value lines="20" type="code" language="python">## Calculates the fluid-particle drag force interaction term
572
 
from math import sqrt, pi
573
 
 
574
 
vfrac_tephra = states["Tephra"].scalar_fields["PhaseVolumeFraction"]
575
 
vfrac_water = states["Water"].scalar_fields["PhaseVolumeFraction"]
576
 
 
577
 
velocity_tephra = states["Tephra"].vector_fields["Velocity"]
578
 
velocity_water = states["Water"].vector_fields["Velocity"]
579
 
 
580
 
density_tephra = states["Tephra"].scalar_fields["Density"]
581
 
density_water = states["Water"].scalar_fields["Density"]
582
 
 
583
 
viscosity_water = states["Water"].tensor_fields["Viscosity"]
584
 
 
585
 
# Particle diameter
586
 
d = 1e-3
587
 
 
588
 
for i in range(field.node_count):
589
 
 
590
 
  # Reynolds number
591
 
  Re = (vfrac_water.node_val(i)*density_water.node_val(i)*sqrt((velocity_water.node_val(i)[0] - velocity_tephra.node_val(i)[0])**2 + (velocity_water.node_val(i)[1] - velocity_tephra.node_val(i)[1])**2)*d)/viscosity_water.node_val(0)[0][0]
592
 
  
593
 
  # Drag coefficient
594
 
  if(Re &lt; 1000):
595
 
    C = (24.0/Re)*(1.0 + 0.15*Re**0.687)
596
 
  else:
597
 
    C = 0.44
598
 
    
599
 
  # Don't let C be NaN
600
 
  if(Re == 0):
601
 
    C = 1e16
602
 
 
603
 
  magnitude = sqrt((velocity_water.node_val(i)[0] - velocity_tephra.node_val(i)[0])**2 + (velocity_water.node_val(i)[1] - velocity_tephra.node_val(i)[1])**2)
604
 
  
605
 
  K = (3.0/4.0)*C*(vfrac_tephra.node_val(i)*vfrac_water.node_val(i)*density_water.node_val(i)*magnitude)/(d*vfrac_water.node_val(i)**2.7)
606
 
 
607
 
  Force = K*(velocity_tephra.node_val(i) - velocity_water.node_val(i))
608
 
 
609
 
  field.set(i, -Force)</string_value>
610
 
          <depends>
611
 
            <string_value lines="1">Tephra::PhaseVolumeFraction,Tephra::Density,Tephra::Velocity,Water::PhaseVolumeFraction,Water::Density,Water::Velocity</string_value>
612
 
          </depends>
613
 
        </algorithm>
614
 
        <mesh name="VelocityMesh"/>
615
 
        <output/>
616
 
        <stat>
617
 
          <include_in_stat/>
618
 
        </stat>
619
 
        <convergence>
620
 
          <include_in_convergence/>
621
 
        </convergence>
622
 
        <detectors>
623
 
          <include_in_detectors/>
624
 
        </detectors>
625
 
        <steady_state>
626
 
          <include_in_steady_state/>
627
 
        </steady_state>
628
 
      </diagnostic>
629
 
    </vector_field>
630
 
    <vector_field name="BuoyancyForce" rank="1">
631
 
      <diagnostic>
632
 
        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
633
 
          <string_value lines="20" type="code" language="python">## Calculates the fluid-particle drag force interaction term
634
 
from math import sqrt, pi
635
 
 
636
 
vfrac_tephra = states["Tephra"].scalar_fields["PhaseVolumeFraction"]
637
 
vfrac_water = states["Water"].scalar_fields["PhaseVolumeFraction"]
638
 
 
639
 
velocity_tephra = states["Tephra"].vector_fields["Velocity"]
640
 
velocity_water = states["Water"].vector_fields["Velocity"]
641
 
 
642
 
density_tephra = states["Tephra"].scalar_fields["Density"]
643
 
density_water = states["Water"].scalar_fields["Density"]
644
 
 
645
 
viscosity_water = states["Water"].tensor_fields["Viscosity"]
646
 
 
647
 
# Particle diameter
648
 
d = 1e-3
649
 
 
650
 
for i in range(field.node_count):
651
 
  
652
 
  Force = vfrac_water.node_val(i)*9.8*vfrac_tephra.node_val(i)*(density_tephra.node_val(i) - density_water.node_val(i))
653
 
 
654
 
  field.set(i, [0.0, -Force])</string_value>
655
 
          <depends>
656
 
            <string_value lines="1">Water::Velocity,Tephra::Velocity,Water::PhaseVolumeFraction,Tephra::PhaseVolumeFraction,Water::Density,Tephra::Density</string_value>
657
 
          </depends>
658
 
        </algorithm>
659
 
        <mesh name="VelocityMesh"/>
660
 
        <output/>
661
 
        <stat>
662
 
          <include_in_stat/>
663
 
        </stat>
664
 
        <convergence>
665
 
          <include_in_convergence/>
666
 
        </convergence>
667
 
        <detectors>
668
 
          <include_in_detectors/>
669
 
        </detectors>
670
 
        <steady_state>
671
 
          <include_in_steady_state/>
672
 
        </steady_state>
673
 
      </diagnostic>
674
 
    </vector_field>
675
 
    <vector_field name="TerminalVelocityWenYu" rank="1">
676
 
      <diagnostic>
677
 
        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
678
 
          <string_value lines="20" type="code" language="python">## Calculates the fluid-particle drag force interaction term
679
 
from math import sqrt, pi
680
 
 
681
 
vfrac_tephra = states["Tephra"].scalar_fields["PhaseVolumeFraction"]
682
 
vfrac_water = states["Water"].scalar_fields["PhaseVolumeFraction"]
683
 
 
684
 
velocity_tephra = states["Tephra"].vector_fields["Velocity"]
685
 
velocity_water = states["Water"].vector_fields["Velocity"]
686
 
 
687
 
density_tephra = states["Tephra"].scalar_fields["Density"]
688
 
density_water = states["Water"].scalar_fields["Density"]
689
 
 
690
 
viscosity_water = states["Water"].tensor_fields["Viscosity"]
691
 
 
692
 
# Particle diameter
693
 
d = 1e-3
694
 
 
695
 
for i in range(field.node_count):
696
 
 
697
 
  # Reynolds number
698
 
  Re = (vfrac_water.node_val(i)*density_water.node_val(i)*sqrt((velocity_water.node_val(i)[0] - velocity_tephra.node_val(i)[0])**2 + (velocity_water.node_val(i)[1] - velocity_tephra.node_val(i)[1])**2)*d)/viscosity_water.node_val(0)[0][0]
699
 
  
700
 
  # Drag coefficient
701
 
  if(Re &lt; 1000):
702
 
    C = (24.0/Re)*(1.0 + 0.15*Re**0.687)
703
 
  else:
704
 
    C = 0.44
705
 
    
706
 
  # Don't let C be NaN
707
 
  if(Re == 0):
708
 
    C = 1e16
709
 
 
710
 
  magnitude = sqrt((velocity_water.node_val(i)[0] - velocity_tephra.node_val(i)[0])**2 + (velocity_water.node_val(i)[1] - velocity_tephra.node_val(i)[1])**2)
711
 
  
712
 
  K = (3.0/4.0)*C*(vfrac_tephra.node_val(i)*vfrac_water.node_val(i)*density_water.node_val(i)*magnitude)/(d*vfrac_water.node_val(i)**2.7)
713
 
    
714
 
  Drag = K
715
 
  BuoyancyForce = vfrac_water.node_val(i)*9.8*vfrac_tephra.node_val(i)*(density_tephra.node_val(i) - density_water.node_val(i))
716
 
  
717
 
  field.set(i, [0.0, -BuoyancyForce/Drag])</string_value>
718
 
          <depends>
719
 
            <string_value lines="1">Water::Velocity,Tephra::Velocity,Water::PhaseVolumeFraction,Tephra::PhaseVolumeFraction,Water::Density,Tephra::Density</string_value>
720
 
          </depends>
721
 
        </algorithm>
722
 
        <mesh name="VelocityMesh"/>
723
 
        <output/>
724
 
        <stat>
725
 
          <include_in_stat/>
726
 
        </stat>
727
 
        <convergence>
728
 
          <include_in_convergence/>
729
 
        </convergence>
730
 
        <detectors>
731
 
          <include_in_detectors/>
732
 
        </detectors>
733
 
        <steady_state>
734
 
          <include_in_steady_state/>
735
 
        </steady_state>
736
 
      </diagnostic>
737
 
    </vector_field>
738
 
    <vector_field name="TerminalVelocityStokes" rank="1">
739
 
      <diagnostic>
740
 
        <algorithm name="vector_python_diagnostic" material_phase_support="multiple">
741
 
          <string_value lines="20" type="code" language="python">## Calculates the fluid-particle drag force interaction term
742
 
from math import sqrt, pi
743
 
 
744
 
vfrac_tephra = states["Tephra"].scalar_fields["PhaseVolumeFraction"]
745
 
vfrac_water = states["Water"].scalar_fields["PhaseVolumeFraction"]
746
 
 
747
 
velocity_tephra = states["Tephra"].vector_fields["Velocity"]
748
 
velocity_water = states["Water"].vector_fields["Velocity"]
749
 
 
750
 
density_tephra = states["Tephra"].scalar_fields["Density"]
751
 
density_water = states["Water"].scalar_fields["Density"]
752
 
 
753
 
viscosity_water = states["Water"].tensor_fields["Viscosity"]
754
 
 
755
 
# Particle diameter
756
 
d = 1e-3
757
 
 
758
 
for i in range(field.node_count):
759
 
  relative_terminal_velocity = vfrac_water.node_val(i)*9.8*(density_tephra.node_val(i) - density_water.node_val(i))*d*d/(18.0*viscosity_water.node_val(0)[0][0])
760
 
  
761
 
  field.set(i, [0.0, -relative_terminal_velocity])</string_value>
762
 
          <depends>
763
 
            <string_value lines="1">Water::Velocity,Tephra::Velocity,Water::PhaseVolumeFraction,Tephra::PhaseVolumeFraction,Water::Density,Tephra::Density</string_value>
764
 
          </depends>
765
 
        </algorithm>
766
 
        <mesh name="VelocityMesh"/>
767
 
        <output/>
768
 
        <stat>
769
 
          <include_in_stat/>
770
 
        </stat>
771
 
        <convergence>
772
 
          <include_in_convergence/>
773
 
        </convergence>
774
 
        <detectors>
775
 
          <include_in_detectors/>
776
 
        </detectors>
777
 
        <steady_state>
778
 
          <include_in_steady_state/>
779
 
        </steady_state>
780
 
      </diagnostic>
781
 
    </vector_field>
782
 
    <multiphase_properties>
783
 
      <particle_diameter>
784
 
        <real_value rank="0">1e-3</real_value>
785
 
      </particle_diameter>
786
 
    </multiphase_properties>
787
 
  </material_phase>
788
 
  <multiphase_interaction>
789
 
    <fluid_particle_drag>
790
 
      <drag_correlation name="wen_yu"/>
791
 
    </fluid_particle_drag>
792
 
  </multiphase_interaction>
793
 
</fluidity_options>