← Back to branch summary

~fluidity-core/fluidity/ocean_forcing

« back to all changes in this revision

Viewing changes to assemble/Multiphase.F90

  • Committer: Christian Jacobs
  • Date: 2011-08-31 20:38:52 UTC
  • mfrom: (3408.1.124 ctjacobs-multiphase)
  • Revision ID: c.jacobs10@imperial.ac.uk-20110831203852-4cburyv2422y6ms1
Merging ctjacobs-multiphase branch revisions (from r3452 to r3532 inclusive) into trunk. Branch queue passes all unit, short and medium tests on buildbot.

Summary of changes:
===================

- Implemented the fluid-particle drag term in Multiphase.F90. The mphase_stokes_law and mphase_tephra_settling test cases have been updated to use this new implementation, rather than a Python diagnostic Source field. An MMS test which includes this fluid-particle drag term is available in longtests/mphase_mms_p1dgp2_fpdrag, and converges at the expected order.

- Implemented a new boundary condition type, called 'flux', for control volume discretisations. Updated the mphase_tephra_settling test case to use this, which allows to flux in through the boundary at the correct rate. Also added documentation regarding this boundary condition type to the manual.

- When checking the tfield_bc_type in Field_Equations_CV.F90, we now compare against the more meaningful integer parameters BC_TYPE_WEAKDIRICHLET (=1), BC_TYPE_NEUMANN (=2), etc.

- Added the 'tephra_settling' example to the examples directory, along with some documentation (background, simulation setup, results, pretty pictures) in manual/examples.tex.

- In manual/model_equations.tex and manual/configuring_fluidity.tex: Added documentation on the multi-phase flow model equations, setting up multi-phase flow simulations, available inter-phase interaction terms, and current limitations of the model.

- Corrected typos in various places in the manual.

- Made the mphase_tephra_settling test more realistic by adding in random perturbations to the PhaseVolumeFraction influx. Also made a few minor tweaks (e.g. finish_time, timestep), and added pass_tests.

- Added a three-phase test called mphase_three_layers. This checks that our model can simulate more than one dispersed phase, allowing us to use several different particle sizes in our tephra settling simulation later on.

- Moved the mphase_tephra_settling_3d test case to the longtests directory.

Show diffs side-by-side

added added

removed removed

Lines of Context:
assemble/Multiphase.F90
24
24
!    License along with this library; if not, write to the Free Software
25
25
!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
26
26
!    USA
 
27
 
27
28
#include "fdebug.h"
28
29
 
29
30
   module multiphase_module
36
37
      use global_parameters, only: OPTION_PATH_LEN
37
38
      use field_priority_lists
38
39
      use field_options
 
40
      use fetools
 
41
      use sparse_tools_petsc
 
42
      use profiler
 
43
 
39
44
      implicit none
40
45
 
41
46
      private
42
47
      public :: get_phase_submaterials, get_nonlinear_volume_fraction, &
43
 
                calculate_diagnostic_phase_volume_fraction
 
48
                calculate_diagnostic_phase_volume_fraction, &
 
49
                add_fluid_particle_drag
44
50
 
45
51
   contains
46
52
 
245
251
         ewrite(1,*) 'Exiting calculate_diagnostic_phase_volume_fraction'
246
252
 
247
253
      end subroutine calculate_diagnostic_phase_volume_fraction
 
254
      
 
255
 
 
256
      !! Multiphase interaction terms, F_i
 
257
      subroutine add_fluid_particle_drag(state, istate, u, x, big_m, mom_rhs)
 
258
         !!< This computes the fluid-particle drag force term.
 
259
         !!< Note that this assumes only one fluid phase, and one or more particle phases.
 
260
         
 
261
         type(state_type), dimension(:), intent(inout) :: state     
 
262
         integer, intent(in) :: istate
 
263
         type(vector_field), intent(in) :: u, x
 
264
         type(petsc_csr_matrix), intent(inout) :: big_m
 
265
         type(vector_field), intent(inout) :: mom_rhs
 
266
         
 
267
         ! Local variables              
 
268
         integer :: ele
 
269
         type(element_type) :: test_function
 
270
         type(element_type), pointer :: u_shape
 
271
         integer, dimension(:), pointer :: u_ele
 
272
         logical :: dg
 
273
             
 
274
         type(vector_field), pointer :: velocity_fluid
 
275
                 
 
276
         logical :: is_particle_phase
 
277
         
 
278
         real :: dt, theta
 
279
         logical, dimension(u%dim, u%dim) :: block_mask ! Control whether the off diagonal entries are used
 
280
                 
 
281
         integer :: i, dim
 
282
         logical :: not_found ! Error flag. Have we found the fluid phase?
 
283
         integer :: istate_fluid, istate_particle
 
284
         
 
285
         
 
286
         ewrite(1, *) "Entering add_fluid_particle_drag"
 
287
               
 
288
         ! Get the timestepping options
 
289
         call get_option("/timestepping/timestep", dt)
 
290
         call get_option(trim(u%option_path)//"/prognostic/temporal_discretisation/theta", &
 
291
                        theta)
 
292
                                          
 
293
         ! For the big_m matrix. Controls whether the off diagonal entries are used    
 
294
         block_mask = .false.
 
295
         do dim = 1, u%dim
 
296
            block_mask(dim, dim) = .true.
 
297
         end do
 
298
 
 
299
         ! Are we using a discontinuous Galerkin discretisation?
 
300
         dg = continuity(u) < 0
 
301
 
 
302
         ! Is this phase a particle phase?
 
303
         is_particle_phase = have_option("/material_phase["//int2str(istate-1)//&
 
304
                          &"]/multiphase_properties/particle_diameter")
 
305
              
 
306
         ! Retrieve the index of the fluid phase in the state array.
 
307
         not_found = .true.
 
308
         if(is_particle_phase) then    
 
309
            do i = 1, size(state)
 
310
               if(.not.have_option("/material_phase["//int2str(i-1)//&
 
311
                        &"]/multiphase_properties/particle_diameter")) then
 
312
 
 
313
                  velocity_fluid => extract_vector_field(state(i), "Velocity")
 
314
                  ! Aliased material_phases will also not have a particle_diameter,
 
315
                  ! so here we make sure that we don't count these as the fluid phase
 
316
                  if(.not.aliased(velocity_fluid)) then
 
317
                     istate_fluid = i
 
318
                     if(.not.not_found) then
 
319
                        FLExit("Fluid-particle drag does not currently support more than one fluid phase.")
 
320
                     end if
 
321
                     not_found = .false.
 
322
                  end if
 
323
 
 
324
               end if
 
325
            end do
 
326
         else
 
327
            istate_fluid = istate
 
328
            not_found = .false.
 
329
         end if
 
330
 
 
331
         if(not_found) then
 
332
            FLAbort("No fluid phase found for the fluid-particle drag.")
 
333
         end if
 
334
 
 
335
         ! If we have a fluid-particle pair, then assemble the drag term
 
336
         if(is_particle_phase) then
 
337
            call assemble_fluid_particle_drag(istate_fluid, istate)
 
338
         else
 
339
            state_loop: do i = 1, size(state)
 
340
               if(i /= istate_fluid) then
 
341
                  call assemble_fluid_particle_drag(istate_fluid, i)
 
342
               end if
 
343
            end do state_loop
 
344
         end if
 
345
         
 
346
         ewrite(1, *) "Exiting add_fluid_particle_drag"
 
347
         
 
348
         contains
 
349
 
 
350
            subroutine assemble_fluid_particle_drag(istate_fluid, istate_particle)
 
351
 
 
352
               integer, intent(in) :: istate_fluid, istate_particle
 
353
 
 
354
               type(scalar_field), pointer :: vfrac_fluid, vfrac_particle
 
355
               type(scalar_field), pointer :: density_fluid, density_particle
 
356
               type(vector_field), pointer :: velocity_fluid, velocity_particle
 
357
               type(vector_field), pointer :: oldu_fluid, oldu_particle
 
358
               type(vector_field), pointer :: nu_fluid, nu_particle ! Non-linear approximation to the Velocities
 
359
               type(tensor_field), pointer :: viscosity_fluid
 
360
               type(scalar_field) :: nvfrac_fluid, nvfrac_particle
 
361
               real :: d ! Particle diameter
 
362
 
 
363
               ! Get the necessary fields to calculate the drag force
 
364
               velocity_fluid => extract_vector_field(state(istate_fluid), "Velocity")
 
365
               velocity_particle => extract_vector_field(state(istate_particle), "Velocity")
 
366
               if(.not.aliased(velocity_particle)) then ! Don't count the aliased material_phases
 
367
                  
 
368
                  vfrac_fluid => extract_scalar_field(state(istate_fluid), "PhaseVolumeFraction")
 
369
                  vfrac_particle => extract_scalar_field(state(istate_particle), "PhaseVolumeFraction")
 
370
                  density_fluid => extract_scalar_field(state(istate_fluid), "Density")
 
371
                  density_particle => extract_scalar_field(state(istate_particle), "Density")
 
372
                  viscosity_fluid => extract_tensor_field(state(istate_fluid), "Viscosity")
 
373
         
 
374
                  call get_option("/material_phase["//int2str(istate_particle-1)//&
 
375
                           &"]/multiphase_properties/particle_diameter", d)
 
376
 
 
377
                  ! Calculate the non-linear approximation to the PhaseVolumeFractions
 
378
                  call allocate(nvfrac_fluid, vfrac_fluid%mesh, "NonlinearPhaseVolumeFraction")
 
379
                  call allocate(nvfrac_particle, vfrac_particle%mesh, "NonlinearPhaseVolumeFraction")
 
380
                  call zero(nvfrac_fluid)
 
381
                  call zero(nvfrac_particle)
 
382
                  call get_nonlinear_volume_fraction(state(istate_fluid), nvfrac_fluid)
 
383
                  call get_nonlinear_volume_fraction(state(istate_particle), nvfrac_particle)
 
384
                  
 
385
                  ! Get the non-linear approximation to the Velocities
 
386
                  nu_fluid => extract_vector_field(state(istate_fluid), "NonlinearVelocity")
 
387
                  nu_particle => extract_vector_field(state(istate_particle), "NonlinearVelocity")
 
388
                  oldu_fluid => extract_vector_field(state(istate_fluid), "OldVelocity")
 
389
                  oldu_particle => extract_vector_field(state(istate_particle), "OldVelocity")
 
390
      
 
391
                  ! ----- Volume integrals over elements -------------           
 
392
                  call profiler_tic(u, "element_loop")
 
393
                  element_loop: do ele = 1, element_count(u)
 
394
 
 
395
                     if(.not.dg .or. (dg .and. element_owned(u,ele))) then
 
396
                        u_ele=>ele_nodes(u, ele)
 
397
                        u_shape => ele_shape(u, ele)
 
398
                        test_function = u_shape                         
 
399
 
 
400
                        call add_fluid_particle_drag_element(ele, test_function, u_shape, &
 
401
                                                            x, u, big_m, mom_rhs, &
 
402
                                                            nvfrac_fluid, nvfrac_particle, &
 
403
                                                            density_fluid, density_particle, &
 
404
                                                            nu_fluid, nu_particle, &
 
405
                                                            oldu_fluid, oldu_particle, &
 
406
                                                            viscosity_fluid, d)
 
407
                     end if
 
408
 
 
409
                  end do element_loop
 
410
                  call profiler_toc(u, "element_loop")
 
411
 
 
412
                  call deallocate(nvfrac_fluid)
 
413
                  call deallocate(nvfrac_particle)
 
414
               end if
 
415
 
 
416
            end subroutine assemble_fluid_particle_drag
 
417
         
 
418
            subroutine add_fluid_particle_drag_element(ele, test_function, u_shape, &
 
419
                                                      x, u, big_m, mom_rhs, &
 
420
                                                      vfrac_fluid, vfrac_particle, &
 
421
                                                      density_fluid, density_particle, &
 
422
                                                      nu_fluid, nu_particle, &
 
423
                                                      oldu_fluid, oldu_particle, &
 
424
                                                      viscosity_fluid, d)
 
425
                                                         
 
426
               integer, intent(in) :: ele
 
427
               type(element_type), intent(in) :: test_function
 
428
               type(element_type), intent(in) :: u_shape
 
429
               type(vector_field), intent(in) :: u, x
 
430
               type(petsc_csr_matrix), intent(inout) :: big_m
 
431
               type(vector_field), intent(inout) :: mom_rhs
 
432
                    
 
433
               type(scalar_field), intent(in) :: vfrac_fluid, vfrac_particle
 
434
               type(scalar_field), intent(in) :: density_fluid, density_particle
 
435
               type(vector_field), intent(in) :: nu_fluid, nu_particle
 
436
               type(vector_field), intent(in) :: oldu_fluid, oldu_particle
 
437
               type(tensor_field), intent(in) :: viscosity_fluid    
 
438
               real, intent(in) :: d ! Particle diameter 
 
439
               
 
440
               ! Local variables
 
441
               real, dimension(ele_ngi(u,ele)) :: vfrac_fluid_gi, vfrac_particle_gi
 
442
               real, dimension(ele_ngi(u,ele)) :: density_fluid_gi, density_particle_gi
 
443
               real, dimension(u%dim, ele_ngi(u,ele)) :: nu_fluid_gi, nu_particle_gi
 
444
               real, dimension(u%dim, u%dim, ele_ngi(u,ele)) :: viscosity_fluid_gi
 
445
               
 
446
               real, dimension(u%dim, ele_loc(u, ele)) :: oldu_val
 
447
               
 
448
               real, dimension(ele_loc(u, ele), ele_ngi(u, ele), x%dim) :: du_t
 
449
               real, dimension(ele_ngi(u, ele)) :: detwei
 
450
               real, dimension(u%dim, ele_loc(u,ele)) :: interaction_rhs_mat
 
451
               real, dimension(ele_loc(u, ele), ele_loc(u, ele)) :: interaction_big_m_mat
 
452
               real, dimension(u%dim, ele_loc(u,ele)) :: rhs_addto
 
453
               real, dimension(u%dim, u%dim, ele_loc(u,ele), ele_loc(u,ele)) :: big_m_tensor_addto
 
454
               
 
455
               real, dimension(ele_ngi(u,ele)) :: particle_re ! Particle Reynolds number
 
456
               real, dimension(ele_ngi(u,ele)) :: drag_coefficient
 
457
               real, dimension(ele_ngi(u,ele)) :: magnitude ! |v_f - v_p|
 
458
               
 
459
               real, dimension(ele_ngi(u,ele)) :: K
 
460
               real, dimension(ele_ngi(u,ele)) :: drag_force_big_m
 
461
               real, dimension(u%dim, ele_ngi(u,ele)) :: drag_force_rhs ! drag_force = K*(v_f - v_p)
 
462
                             
 
463
               integer :: dim, gi
 
464
               
 
465
               ! Compute detwei
 
466
               call transform_to_physical(x, ele, u_shape, dshape=du_t, detwei=detwei)
 
467
               
 
468
               ! Get the values of the necessary fields at the Gauss points
 
469
               vfrac_fluid_gi = ele_val_at_quad(vfrac_fluid, ele)
 
470
               vfrac_particle_gi = ele_val_at_quad(vfrac_particle, ele)
 
471
               density_fluid_gi = ele_val_at_quad(density_fluid, ele)
 
472
               density_particle_gi = ele_val_at_quad(density_particle, ele)
 
473
               nu_fluid_gi = ele_val_at_quad(nu_fluid, ele)
 
474
               nu_particle_gi = ele_val_at_quad(nu_particle, ele)
 
475
               viscosity_fluid_gi = ele_val_at_quad(viscosity_fluid, ele)               
 
476
         
 
477
               ! Compute the magnitude of the relative velocity
 
478
               do gi = 1, ele_ngi(u,ele)
 
479
                  magnitude(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi))
 
480
               end do
 
481
 
 
482
               ! Compute the particle Reynolds number
 
483
               ! (Assumes isotropic viscosity for now)
 
484
               particle_re = (vfrac_fluid_gi*density_fluid_gi*magnitude*d) / viscosity_fluid_gi(1,1,:)
 
485
           
 
486
               ! Compute the drag coefficient
 
487
               do gi = 1, ele_ngi(u,ele)
 
488
                  if(particle_re(gi) < 1000) then
 
489
                     drag_coefficient(gi) = (24.0/particle_re(gi))
 
490
                     ! For the Wen & Yu (1966) drag term, this should be:
 
491
                     ! drag_coefficient(gi) = (24.0/particle_re(gi))*(1.0+0.15*particle_re(gi)**0.687)
 
492
                  else
 
493
                     drag_coefficient(gi) = 0.44
 
494
                  end if
 
495
               end do
 
496
                      
 
497
               ! Don't let the drag_coefficient be NaN
 
498
               do gi = 1, ele_ngi(u,ele)
 
499
                  if(particle_re(gi) == 0) then
 
500
                     drag_coefficient(gi) = 1e16
 
501
                  end if
 
502
               end do
 
503
           
 
504
               ! For the Wen & Yu (1966) drag term, this should be:
 
505
               ! K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient*(vfrac_fluid_gi*density_fluid_gi*magnitude)/(d*vfrac_fluid_gi**2.7)
 
506
               K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient*(vfrac_fluid_gi*density_fluid_gi*magnitude)/(d)
 
507
               
 
508
               if(is_particle_phase) then
 
509
                  drag_force_big_m = -K
 
510
                  do dim = 1, u%dim
 
511
                     drag_force_rhs(dim,:) = -K*(-nu_fluid_gi(dim,:))
 
512
                  end do
 
513
               else
 
514
                  drag_force_big_m = -K
 
515
                  do dim = 1, u%dim
 
516
                     drag_force_rhs(dim,:) = K*(nu_particle_gi(dim,:))
 
517
                  end do
 
518
               end if
 
519
               
 
520
               ! Form the element interaction/drag matrix
 
521
               interaction_big_m_mat = shape_shape(test_function, u_shape, detwei*drag_force_big_m)
 
522
               interaction_rhs_mat = shape_vector_rhs(test_function, drag_force_rhs, detwei)
 
523
              
 
524
               ! Add contribution  
 
525
               big_m_tensor_addto = 0.0            
 
526
               rhs_addto = 0.0
 
527
               if(is_particle_phase) then
 
528
                  oldu_val = ele_val(oldu_particle, ele)
 
529
               else
 
530
                  oldu_val = ele_val(oldu_fluid, ele)
 
531
               end if
 
532
               do dim = 1, u%dim
 
533
                  big_m_tensor_addto(dim, dim, :, :) = big_m_tensor_addto(dim, dim, :, :) - dt*theta*interaction_big_m_mat
 
534
                  rhs_addto(dim,:) = rhs_addto(dim,:) + matmul(interaction_big_m_mat, oldu_val(dim,:)) + interaction_rhs_mat(dim,:)
 
535
               end do
 
536
               
 
537
               ! Add the contribution to mom_rhs
 
538
               call addto(mom_rhs, u_ele, rhs_addto) 
 
539
               ! Add to the big_m matrix
 
540
               call addto(big_m, u_ele, u_ele, big_m_tensor_addto, block_mask=block_mask)
 
541
 
 
542
            end subroutine add_fluid_particle_drag_element
 
543
 
 
544
      end subroutine add_fluid_particle_drag
248
545
 
249
546
   end module multiphase_module