← Back to branch summary

~fluidity-core/fluidity/excise-fldecomp

« back to all changes in this revision

Viewing changes to assemble/Multiphase.F90

  • Committer: Mark Filipiak
  • Date: 2012-08-13 11:42:30 UTC
  • mfrom: (4003.1.23 dev-trunk)
  • Revision ID: mjf@staffmail.ed.ac.uk-20120813114230-wzoyf2gi4p4oxeh4
Merge in of the latest trunk.  To try to cure non-flredecomp tests that are passing at EPCC but failing in buildbot.

Show diffs side-by-side

added added

removed removed

Lines of Context:
assemble/Multiphase.F90
46
46
      private
47
47
      public :: get_phase_submaterials, get_nonlinear_volume_fraction, &
48
48
                calculate_diagnostic_phase_volume_fraction, &
49
 
                add_fluid_particle_drag
 
49
                add_fluid_particle_drag, add_heat_transfer
50
50
 
51
51
   contains
52
52
 
268
268
         integer :: ele
269
269
         type(element_type) :: test_function
270
270
         type(element_type), pointer :: u_shape
271
 
         integer, dimension(:), pointer :: u_ele
 
271
         integer, dimension(:), pointer :: u_nodes
272
272
         logical :: dg
273
273
             
274
274
         type(vector_field), pointer :: velocity_fluid
281
281
         integer :: i, dim
282
282
         logical :: not_found ! Error flag. Have we found the fluid phase?
283
283
         integer :: istate_fluid, istate_particle
284
 
         
 
284
   
 
285
         ! Types of drag correlation
 
286
         integer, parameter :: DRAG_CORRELATION_TYPE_STOKES = 1, DRAG_CORRELATION_TYPE_WEN_YU = 2, DRAG_CORRELATION_TYPE_ERGUN = 3
285
287
         
286
288
         ewrite(1, *) "Entering add_fluid_particle_drag"
 
289
         
 
290
         ! Let's check whether we actually have at least one particle phase.
 
291
         if(option_count("/material_phase/multiphase_properties/particle_diameter") == 0) then
 
292
            FLExit("Fluid-particle drag enabled but no particle_diameter has been specified for the particle phase(s).")
 
293
         end if
287
294
               
288
295
         ! Get the timestepping options
289
296
         call get_option("/timestepping/timestep", dt)
290
 
         call get_option(trim(u%option_path)//"/prognostic/temporal_discretisation/theta", &
291
 
                        theta)
 
297
         call get_option(trim(u%option_path)//"/prognostic/temporal_discretisation/theta", theta)
292
298
                                          
293
299
         ! For the big_m matrix. Controls whether the off diagonal entries are used    
294
300
         block_mask = .false.
300
306
         dg = continuity(u) < 0
301
307
 
302
308
         ! Is this phase a particle phase?
303
 
         is_particle_phase = have_option("/material_phase["//int2str(istate-1)//&
304
 
                          &"]/multiphase_properties/particle_diameter")
 
309
         is_particle_phase = have_option(trim(state(istate)%option_path)//"/multiphase_properties/particle_diameter")
305
310
              
306
311
         ! Retrieve the index of the fluid phase in the state array.
307
312
         not_found = .true.
308
313
         if(is_particle_phase) then    
309
314
            do i = 1, size(state)
310
 
               if(.not.have_option("/material_phase["//int2str(i-1)//&
311
 
                        &"]/multiphase_properties/particle_diameter")) then
 
315
               if(.not.have_option(trim(state(i)%option_path)//"/multiphase_properties/particle_diameter")) then
312
316
 
313
317
                  velocity_fluid => extract_vector_field(state(i), "Velocity")
314
318
                  ! Aliased material_phases will also not have a particle_diameter,
359
363
               type(tensor_field), pointer :: viscosity_fluid
360
364
               type(scalar_field) :: nvfrac_fluid, nvfrac_particle
361
365
               real :: d ! Particle diameter
 
366
               character(len=OPTION_PATH_LEN) :: drag_correlation_name
 
367
               integer :: drag_correlation
362
368
 
363
369
               ! Get the necessary fields to calculate the drag force
364
370
               velocity_fluid => extract_vector_field(state(istate_fluid), "Velocity")
371
377
                  density_particle => extract_scalar_field(state(istate_particle), "Density")
372
378
                  viscosity_fluid => extract_tensor_field(state(istate_fluid), "Viscosity")
373
379
         
374
 
                  call get_option("/material_phase["//int2str(istate_particle-1)//&
375
 
                           &"]/multiphase_properties/particle_diameter", d)
 
380
                  call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter", d)
376
381
 
377
382
                  ! Calculate the non-linear approximation to the PhaseVolumeFractions
378
383
                  call allocate(nvfrac_fluid, vfrac_fluid%mesh, "NonlinearPhaseVolumeFraction")
387
392
                  nu_particle => extract_vector_field(state(istate_particle), "NonlinearVelocity")
388
393
                  oldu_fluid => extract_vector_field(state(istate_fluid), "OldVelocity")
389
394
                  oldu_particle => extract_vector_field(state(istate_particle), "OldVelocity")
390
 
      
 
395
                  
 
396
                  call get_option("/multiphase_interaction/fluid_particle_drag/drag_correlation/name", drag_correlation_name)
 
397
                  select case(trim(drag_correlation_name))
 
398
                     case("stokes")
 
399
                        drag_correlation = DRAG_CORRELATION_TYPE_STOKES
 
400
                     case("wen_yu")
 
401
                        drag_correlation = DRAG_CORRELATION_TYPE_WEN_YU
 
402
                     case("ergun")
 
403
                        drag_correlation = DRAG_CORRELATION_TYPE_ERGUN
 
404
                     case("default")
 
405
                        FLAbort("Unknown correlation for fluid-particle drag")
 
406
                  end select
 
407
                  
391
408
                  ! ----- Volume integrals over elements -------------           
392
409
                  call profiler_tic(u, "element_loop")
393
410
                  element_loop: do ele = 1, element_count(u)
394
411
 
395
412
                     if(.not.dg .or. (dg .and. element_owned(u,ele))) then
396
 
                        u_ele=>ele_nodes(u, ele)
 
413
                        u_nodes => ele_nodes(u, ele)
397
414
                        u_shape => ele_shape(u, ele)
398
415
                        test_function = u_shape                         
399
416
 
403
420
                                                            density_fluid, density_particle, &
404
421
                                                            nu_fluid, nu_particle, &
405
422
                                                            oldu_fluid, oldu_particle, &
406
 
                                                            viscosity_fluid, d)
 
423
                                                            viscosity_fluid, d, drag_correlation)
407
424
                     end if
408
425
 
409
426
                  end do element_loop
421
438
                                                      density_fluid, density_particle, &
422
439
                                                      nu_fluid, nu_particle, &
423
440
                                                      oldu_fluid, oldu_particle, &
424
 
                                                      viscosity_fluid, d)
 
441
                                                      viscosity_fluid, d, drag_correlation)
425
442
                                                         
426
443
               integer, intent(in) :: ele
427
444
               type(element_type), intent(in) :: test_function
436
453
               type(vector_field), intent(in) :: oldu_fluid, oldu_particle
437
454
               type(tensor_field), intent(in) :: viscosity_fluid    
438
455
               real, intent(in) :: d ! Particle diameter 
 
456
               integer, intent(in) :: drag_correlation
439
457
               
440
458
               ! Local variables
441
459
               real, dimension(ele_ngi(u,ele)) :: vfrac_fluid_gi, vfrac_particle_gi
447
465
               
448
466
               real, dimension(ele_loc(u, ele), ele_ngi(u, ele), x%dim) :: du_t
449
467
               real, dimension(ele_ngi(u, ele)) :: detwei
450
 
               real, dimension(u%dim, ele_loc(u,ele)) :: interaction_rhs_mat
 
468
               real, dimension(u%dim, ele_loc(u,ele)) :: interaction_rhs
451
469
               real, dimension(ele_loc(u, ele), ele_loc(u, ele)) :: interaction_big_m_mat
452
470
               real, dimension(u%dim, ele_loc(u,ele)) :: rhs_addto
453
471
               real, dimension(u%dim, u%dim, ele_loc(u,ele), ele_loc(u,ele)) :: big_m_tensor_addto
454
472
               
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|
 
473
               real, dimension(ele_ngi(u,ele)) :: particle_re_gi ! Particle Reynolds number
 
474
               real, dimension(ele_ngi(u,ele)) :: drag_coefficient_gi
 
475
               real, dimension(ele_ngi(u,ele)) :: magnitude_gi ! |v_f - v_p|
458
476
               
459
477
               real, dimension(ele_ngi(u,ele)) :: K
460
478
               real, dimension(ele_ngi(u,ele)) :: drag_force_big_m
476
494
         
477
495
               ! Compute the magnitude of the relative velocity
478
496
               do gi = 1, ele_ngi(u,ele)
479
 
                  magnitude(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi))
 
497
                  magnitude_gi(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi))
480
498
               end do
481
499
 
482
500
               ! Compute the particle Reynolds number
483
501
               ! (Assumes isotropic viscosity for now)
484
 
               particle_re = (vfrac_fluid_gi*density_fluid_gi*magnitude*d) / viscosity_fluid_gi(1,1,:)
 
502
               particle_re_gi = (vfrac_fluid_gi*density_fluid_gi*magnitude_gi*d) / viscosity_fluid_gi(1,1,:)
485
503
           
486
504
               ! 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
 
505
               select case(drag_correlation)
 
506
                  case(DRAG_CORRELATION_TYPE_STOKES)
 
507
                     ! Stokes drag correlation
 
508
                     do gi = 1, ele_ngi(u,ele)
 
509
                        if(particle_re_gi(gi) < 1000) then
 
510
                           drag_coefficient_gi(gi) = (24.0/particle_re_gi(gi))
 
511
                        else
 
512
                           drag_coefficient_gi(gi) = 0.44
 
513
                        end if
 
514
                     end do
 
515
                     
 
516
                  case(DRAG_CORRELATION_TYPE_WEN_YU)
 
517
                     ! Wen & Yu (1966) drag correlation
 
518
                     do gi = 1, ele_ngi(u,ele)
 
519
                        if(particle_re_gi(gi) < 1000) then
 
520
                           drag_coefficient_gi(gi) = (24.0/particle_re_gi(gi))*(1.0+0.15*particle_re_gi(gi)**0.687)
 
521
                        else
 
522
                           drag_coefficient_gi(gi) = 0.44
 
523
                        end if
 
524
                     end do
 
525
                     
 
526
                  case(DRAG_CORRELATION_TYPE_ERGUN)
 
527
                     ! No drag coefficient is needed here.                  
 
528
               end select
496
529
                      
497
 
               ! Don't let the drag_coefficient be NaN
 
530
               ! Don't let the drag_coefficient_gi be NaN
498
531
               do gi = 1, ele_ngi(u,ele)
499
 
                  if(particle_re(gi) == 0) then
500
 
                     drag_coefficient(gi) = 1e16
 
532
                  if(particle_re_gi(gi) == 0) then
 
533
                     drag_coefficient_gi(gi) = 1e16
501
534
                  end if
502
535
               end do
503
536
           
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)
 
537
               select case(drag_correlation)
 
538
                  case(DRAG_CORRELATION_TYPE_STOKES)
 
539
                     K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(vfrac_fluid_gi*density_fluid_gi*magnitude_gi)/(d)
 
540
                  case(DRAG_CORRELATION_TYPE_WEN_YU)
 
541
                     ! Wen & Yu (1966) drag term
 
542
                     K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(vfrac_fluid_gi*density_fluid_gi*magnitude_gi)/(d*vfrac_fluid_gi**2.7)
 
543
                  case(DRAG_CORRELATION_TYPE_ERGUN)
 
544
                     K = 150.0*((vfrac_particle_gi**2)*viscosity_fluid_gi(1,1,:))/(vfrac_fluid_gi*(d**2)) + 1.75*(vfrac_particle_gi*density_fluid_gi*magnitude_gi/d)
 
545
               end select               
507
546
               
508
547
               if(is_particle_phase) then
509
548
                  drag_force_big_m = -K
519
558
               
520
559
               ! Form the element interaction/drag matrix
521
560
               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)
 
561
               interaction_rhs = shape_vector_rhs(test_function, drag_force_rhs, detwei)
523
562
              
524
563
               ! Add contribution  
525
564
               big_m_tensor_addto = 0.0            
531
570
               end if
532
571
               do dim = 1, u%dim
533
572
                  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,:)
 
573
                  rhs_addto(dim,:) = rhs_addto(dim,:) + matmul(interaction_big_m_mat, oldu_val(dim,:)) + interaction_rhs(dim,:)
535
574
               end do
536
575
               
537
576
               ! Add the contribution to mom_rhs
538
 
               call addto(mom_rhs, u_ele, rhs_addto) 
 
577
               call addto(mom_rhs, u_nodes, rhs_addto) 
539
578
               ! Add to the big_m matrix
540
 
               call addto(big_m, u_ele, u_ele, big_m_tensor_addto, block_mask=block_mask)
 
579
               call addto(big_m, u_nodes, u_nodes, big_m_tensor_addto, block_mask=block_mask)
541
580
 
542
581
            end subroutine add_fluid_particle_drag_element
543
582
 
544
583
      end subroutine add_fluid_particle_drag
545
584
 
 
585
      
 
586
      !! Multiphase energy interaction term Q_i
 
587
      !! to be added to the RHS of the internal energy equation.
 
588
      subroutine add_heat_transfer(state, istate, internal_energy, matrix, rhs)
 
589
         !!< This computes the inter-phase heat transfer term.
 
590
         !!< Only between fluid and particle phase pairs.
 
591
         !!< Uses the empirical correlation by Gunn (1978).
 
592
         
 
593
         type(state_type), dimension(:), intent(inout) :: state     
 
594
         integer, intent(in) :: istate
 
595
         type(scalar_field), intent(in) :: internal_energy         
 
596
         type(csr_matrix), intent(inout) :: matrix
 
597
         type(scalar_field), intent(inout) :: rhs
 
598
         
 
599
         ! Local variables              
 
600
         integer :: ele
 
601
         type(element_type) :: test_function
 
602
         type(element_type), pointer :: internal_energy_shape
 
603
         integer, dimension(:), pointer :: internal_energy_nodes
 
604
         logical :: dg
 
605
             
 
606
         type(vector_field), pointer :: x, velocity_fluid
 
607
         
 
608
         real :: dt, theta
 
609
         
 
610
         logical :: is_particle_phase
 
611
         logical :: not_found ! Error flag. Have we found the fluid phase?
 
612
         integer :: i, istate_fluid, istate_particle
 
613
         
 
614
         
 
615
         ewrite(1, *) "Entering add_heat_transfer"
 
616
               
 
617
         ! Get the timestepping options
 
618
         call get_option("/timestepping/timestep", dt)
 
619
         call get_option(trim(internal_energy%option_path)//"/prognostic/temporal_discretisation/theta", &
 
620
                        theta)
 
621
         
 
622
         ! Get the coordinate field from state(istate)
 
623
         x => extract_vector_field(state(istate), "Coordinate")
 
624
         ewrite_minmax(x)
 
625
         assert(x%dim == mesh_dim(internal_energy))
 
626
         assert(ele_count(x) == ele_count(internal_energy))
 
627
 
 
628
         ! Are we using a discontinuous Galerkin discretisation?
 
629
         dg = continuity(internal_energy) < 0
 
630
 
 
631
         ! Is this phase a particle phase?
 
632
         is_particle_phase = have_option(trim(state(istate)%option_path)//"/multiphase_properties/particle_diameter")
 
633
              
 
634
         ! Retrieve the index of the fluid phase in the state array.
 
635
         not_found = .true.
 
636
         if(is_particle_phase) then    
 
637
            do i = 1, size(state)
 
638
               if(.not.have_option(trim(state(i)%option_path)//"/multiphase_properties/particle_diameter")) then
 
639
 
 
640
                  velocity_fluid => extract_vector_field(state(i), "Velocity")
 
641
                  ! Aliased material_phases will also not have a particle_diameter,
 
642
                  ! so here we make sure that we don't count these as the fluid phase
 
643
                  if(.not.aliased(velocity_fluid)) then
 
644
                     istate_fluid = i
 
645
                     if(.not.not_found) then
 
646
                        FLExit("Heat transfer term does not currently support more than one fluid phase.")
 
647
                     end if
 
648
                     not_found = .false.
 
649
                  end if
 
650
 
 
651
               end if
 
652
            end do
 
653
         else
 
654
            istate_fluid = istate
 
655
            not_found = .false.
 
656
         end if
 
657
 
 
658
         if(not_found) then
 
659
            FLExit("No fluid phase found for the heat transfer term.")
 
660
         end if
 
661
 
 
662
         ! If we have a fluid-particle pair, then assemble the heat transfer term
 
663
         if(is_particle_phase) then
 
664
            call assemble_heat_transfer(istate_fluid, istate)
 
665
         else
 
666
            state_loop: do i = 1, size(state)
 
667
               if(i /= istate_fluid) then
 
668
                  call assemble_heat_transfer(istate_fluid, i)
 
669
               end if
 
670
            end do state_loop
 
671
         end if
 
672
         
 
673
         ewrite(1, *) "Exiting add_heat_transfer"
 
674
         
 
675
         contains
 
676
 
 
677
            subroutine assemble_heat_transfer(istate_fluid, istate_particle)
 
678
 
 
679
               integer, intent(in) :: istate_fluid, istate_particle
 
680
 
 
681
               type(scalar_field), pointer :: vfrac_fluid, vfrac_particle
 
682
               type(scalar_field), pointer :: density_fluid, density_particle
 
683
               type(vector_field), pointer :: velocity_fluid, velocity_particle
 
684
               type(scalar_field), pointer :: internal_energy_fluid, internal_energy_particle
 
685
               type(scalar_field), pointer :: old_internal_energy_fluid, old_internal_energy_particle
 
686
               type(vector_field), pointer :: nu_fluid, nu_particle ! Non-linear approximation to the Velocities
 
687
               type(tensor_field), pointer :: viscosity_fluid
 
688
               type(scalar_field) :: nvfrac_fluid, nvfrac_particle
 
689
               real :: d ! Particle diameter
 
690
               real :: k ! Effective gas conductivity
 
691
               real :: C_fluid, C_particle ! Specific heat of the fluid and particle phases at constant volume
 
692
               real :: gamma ! Ratio of specific heats for compressible phase
 
693
               integer :: kstat, cstat_fluid, cstat_particle, gstat
 
694
 
 
695
               ! Get the necessary fields to calculate the heat transfer term
 
696
               velocity_fluid => extract_vector_field(state(istate_fluid), "Velocity")
 
697
               velocity_particle => extract_vector_field(state(istate_particle), "Velocity")
 
698
               if(.not.aliased(velocity_particle)) then ! Don't count the aliased material_phases
 
699
                  
 
700
                  vfrac_fluid => extract_scalar_field(state(istate_fluid), "PhaseVolumeFraction")
 
701
                  vfrac_particle => extract_scalar_field(state(istate_particle), "PhaseVolumeFraction")
 
702
                  density_fluid => extract_scalar_field(state(istate_fluid), "Density")
 
703
                  density_particle => extract_scalar_field(state(istate_particle), "Density")
 
704
                  viscosity_fluid => extract_tensor_field(state(istate_fluid), "Viscosity")
 
705
         
 
706
                  call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter", d)
 
707
                           
 
708
                  call get_option(trim(state(istate_fluid)%option_path)//"/multiphase_properties/effective_conductivity", k, kstat)
 
709
                           
 
710
                  call get_option(trim(state(istate_fluid)%option_path)//"/multiphase_properties/specific_heat", C_fluid, cstat_fluid)
 
711
                  
 
712
                  call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/specific_heat", C_particle, cstat_particle)
 
713
                  
 
714
                  call get_option(trim(state(istate_fluid)%option_path)//"/equation_of_state/compressible/stiffened_gas/ratio_specific_heats", gamma, gstat)
 
715
                           
 
716
                  if(kstat /= 0) then
 
717
                     FLExit("For inter-phase heat transfer, an effective_conductivity needs to be specified for the fluid phase.")
 
718
                  end if                  
 
719
                  if(cstat_fluid /= 0 .or. cstat_particle /= 0) then
 
720
                     FLExit("For inter-phase heat transfer, a specific_heat needs to be specified for each phase.")
 
721
                  end if
 
722
                  if(gstat /= 0) then
 
723
                     FLExit("For inter-phase heat transfer, ratio_specific_heats needs to be specified for the compressible phase.")
 
724
                  end if
 
725
 
 
726
                  ! Calculate the non-linear approximation to the PhaseVolumeFractions
 
727
                  call allocate(nvfrac_fluid, vfrac_fluid%mesh, "NonlinearPhaseVolumeFraction")
 
728
                  call allocate(nvfrac_particle, vfrac_particle%mesh, "NonlinearPhaseVolumeFraction")
 
729
                  call zero(nvfrac_fluid)
 
730
                  call zero(nvfrac_particle)
 
731
                  call get_nonlinear_volume_fraction(state(istate_fluid), nvfrac_fluid)
 
732
                  call get_nonlinear_volume_fraction(state(istate_particle), nvfrac_particle)
 
733
                  
 
734
                  ! Get the non-linear approximation to the Velocities
 
735
                  nu_fluid => extract_vector_field(state(istate_fluid), "NonlinearVelocity")
 
736
                  nu_particle => extract_vector_field(state(istate_particle), "NonlinearVelocity")
 
737
                  
 
738
                  ! Get the current and old internal energy fields
 
739
                  internal_energy_fluid => extract_scalar_field(state(istate_fluid), "InternalEnergy")
 
740
                  internal_energy_particle => extract_scalar_field(state(istate_particle), "InternalEnergy")
 
741
                  old_internal_energy_fluid => extract_scalar_field(state(istate_fluid), "OldInternalEnergy")
 
742
                  old_internal_energy_particle => extract_scalar_field(state(istate_particle), "OldInternalEnergy")
 
743
      
 
744
                  ! ----- Volume integrals over elements -------------           
 
745
                  call profiler_tic(internal_energy, "element_loop")
 
746
                  element_loop: do ele = 1, element_count(internal_energy)
 
747
 
 
748
                     if(.not.dg .or. (dg .and. element_owned(internal_energy,ele))) then
 
749
                        internal_energy_nodes => ele_nodes(internal_energy, ele)
 
750
                        internal_energy_shape => ele_shape(internal_energy, ele)
 
751
                        test_function = internal_energy_shape                         
 
752
 
 
753
                        call add_heat_transfer_element(ele, test_function, internal_energy_shape, &
 
754
                                                      x, internal_energy, matrix, rhs, &
 
755
                                                      nvfrac_fluid, nvfrac_particle, &
 
756
                                                      density_fluid, density_particle, &
 
757
                                                      nu_fluid, nu_particle, &
 
758
                                                      internal_energy_fluid, &
 
759
                                                      internal_energy_particle, &
 
760
                                                      old_internal_energy_fluid, &
 
761
                                                      old_internal_energy_particle, &
 
762
                                                      viscosity_fluid, d, k, C_fluid, &
 
763
                                                      C_particle, gamma)
 
764
                     end if
 
765
 
 
766
                  end do element_loop
 
767
                  call profiler_toc(internal_energy, "element_loop")
 
768
 
 
769
                  call deallocate(nvfrac_fluid)
 
770
                  call deallocate(nvfrac_particle)
 
771
               end if
 
772
 
 
773
            end subroutine assemble_heat_transfer
 
774
         
 
775
            subroutine add_heat_transfer_element(ele, test_function, internal_energy_shape, &
 
776
                                                x, internal_energy, matrix, rhs, &
 
777
                                                vfrac_fluid, vfrac_particle, &
 
778
                                                density_fluid, density_particle, &
 
779
                                                nu_fluid, nu_particle, &
 
780
                                                internal_energy_fluid, &
 
781
                                                internal_energy_particle, &
 
782
                                                old_internal_energy_fluid, &
 
783
                                                old_internal_energy_particle, &
 
784
                                                viscosity_fluid, d, k, C_fluid, &
 
785
                                                C_particle, gamma)
 
786
                                                         
 
787
               integer, intent(in) :: ele
 
788
               type(element_type), intent(in) :: test_function
 
789
               type(element_type), intent(in) :: internal_energy_shape
 
790
               type(vector_field), intent(in) :: x
 
791
               type(scalar_field), intent(in) :: internal_energy
 
792
               type(csr_matrix), intent(inout) :: matrix
 
793
               type(scalar_field), intent(inout) :: rhs
 
794
                    
 
795
               type(scalar_field), intent(in) :: vfrac_fluid, vfrac_particle
 
796
               type(scalar_field), intent(in) :: density_fluid, density_particle
 
797
               type(vector_field), intent(in) :: nu_fluid, nu_particle
 
798
               type(scalar_field), intent(in) :: internal_energy_fluid, internal_energy_particle
 
799
               type(scalar_field), intent(in) :: old_internal_energy_fluid, old_internal_energy_particle
 
800
               type(tensor_field), intent(in) :: viscosity_fluid    
 
801
               real, intent(in) :: d, k, C_fluid, C_particle, gamma
 
802
               
 
803
               ! Local variables
 
804
               real, dimension(ele_ngi(x,ele)) :: internal_energy_fluid_gi, internal_energy_particle_gi
 
805
               real, dimension(ele_ngi(x,ele)) :: vfrac_fluid_gi, vfrac_particle_gi
 
806
               real, dimension(ele_ngi(x,ele)) :: density_fluid_gi, density_particle_gi
 
807
               real, dimension(x%dim, ele_ngi(x,ele)) :: nu_fluid_gi, nu_particle_gi
 
808
               real, dimension(x%dim, x%dim, ele_ngi(x,ele)) :: viscosity_fluid_gi
 
809
               
 
810
               real, dimension(ele_loc(internal_energy,ele)) :: old_internal_energy_val
 
811
               
 
812
               real, dimension(ele_ngi(x,ele)) :: detwei
 
813
               real, dimension(ele_ngi(x,ele)) :: coefficient_for_matrix, coefficient_for_rhs
 
814
               real, dimension(ele_loc(internal_energy,ele)) :: rhs_addto
 
815
               real, dimension(ele_loc(internal_energy,ele), ele_loc(internal_energy,ele)) :: matrix_addto
 
816
               
 
817
               real, dimension(ele_ngi(x,ele)) :: particle_Re ! Particle Reynolds number
 
818
               real, dimension(ele_ngi(x,ele)) :: Pr ! Prandtl number
 
819
               real, dimension(ele_ngi(x,ele)) :: particle_Nu ! Particle Nusselt number
 
820
               real, dimension(ele_ngi(x,ele)) :: velocity_magnitude ! |v_f - v_p|
 
821
               
 
822
               real, dimension(ele_ngi(x,ele)) :: Q ! heat transfer term = Q*(T_p - T_f)
 
823
               real, dimension(ele_loc(internal_energy,ele), ele_loc(internal_energy,ele)) :: heat_transfer_matrix
 
824
               real, dimension(ele_loc(internal_energy,ele)) :: heat_transfer_rhs
 
825
               
 
826
               integer ::  gi
 
827
               
 
828
               ! Compute detwei
 
829
               call transform_to_physical(x, ele, detwei)
 
830
               
 
831
               ! Get the values of the necessary fields at the Gauss points
 
832
               vfrac_fluid_gi = ele_val_at_quad(vfrac_fluid, ele)
 
833
               vfrac_particle_gi = ele_val_at_quad(vfrac_particle, ele)
 
834
               density_fluid_gi = ele_val_at_quad(density_fluid, ele)
 
835
               density_particle_gi = ele_val_at_quad(density_particle, ele)
 
836
               nu_fluid_gi = ele_val_at_quad(nu_fluid, ele)
 
837
               nu_particle_gi = ele_val_at_quad(nu_particle, ele)
 
838
               viscosity_fluid_gi = ele_val_at_quad(viscosity_fluid, ele)  
 
839
               internal_energy_fluid_gi = ele_val_at_quad(internal_energy_fluid, ele)
 
840
               internal_energy_particle_gi = ele_val_at_quad(internal_energy_particle, ele)
 
841
         
 
842
               ! Compute the magnitude of the relative velocity
 
843
               do gi = 1, ele_ngi(x,ele)
 
844
                  velocity_magnitude(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi))
 
845
               end do
 
846
 
 
847
               ! Compute the particle Reynolds number
 
848
               ! (Assumes isotropic viscosity for now)
 
849
               particle_Re = (density_fluid_gi*velocity_magnitude*d) / viscosity_fluid_gi(1,1,:)
 
850
           
 
851
               ! Compute the Prandtl number
 
852
               ! (Assumes isotropic viscosity for now)
 
853
               ! Note: C_fluid (at constant volume) multiplied by gamma = C_fluid at constant pressure
 
854
               Pr = C_fluid*gamma*viscosity_fluid_gi(1,1,:)/k
 
855
               
 
856
               particle_Nu = (7.0 - 10.0*vfrac_fluid_gi + 5.0*vfrac_fluid_gi**2)*(1.0 + 0.7*(particle_Re**0.2)*(Pr**(1.0/3.0))) + &
 
857
                              (1.33 - 2.4*vfrac_fluid_gi + 1.2*vfrac_fluid_gi**2)*(particle_Re**0.7)*(Pr**(1.0/3.0))
 
858
 
 
859
               Q = (6.0*k*vfrac_particle_gi*particle_Nu)/(d**2)
 
860
               
 
861
               ! Note that the transfer term is defined in terms of temperatures (T_fluid and T_particle)
 
862
               ! Let's convert the temperatures to internal energy (E) using E_i = C_i*T_i,
 
863
               ! where C is the specific heat of phase i at constant volume.
 
864
               if(is_particle_phase) then
 
865
                  coefficient_for_matrix = -Q/C_particle
 
866
                  coefficient_for_rhs = -Q*(-internal_energy_fluid_gi/C_fluid)
 
867
               else
 
868
                  coefficient_for_matrix = -Q/C_fluid
 
869
                  coefficient_for_rhs = Q*(internal_energy_particle_gi/C_particle)
 
870
               end if
 
871
               
 
872
               ! Form the element heat transfer matrix and RHS
 
873
               heat_transfer_matrix = shape_shape(test_function, internal_energy_shape, detwei*coefficient_for_matrix)
 
874
               heat_transfer_rhs = shape_rhs(test_function, coefficient_for_rhs*detwei)
 
875
              
 
876
               ! Add contribution  
 
877
               matrix_addto = 0.0            
 
878
               rhs_addto = 0.0
 
879
               if(is_particle_phase) then
 
880
                  old_internal_energy_val = ele_val(old_internal_energy_particle, ele)
 
881
               else
 
882
                  old_internal_energy_val = ele_val(old_internal_energy_fluid, ele)
 
883
               end if
 
884
 
 
885
               matrix_addto = matrix_addto - dt*theta*heat_transfer_matrix
 
886
               rhs_addto = rhs_addto + matmul(heat_transfer_matrix, old_internal_energy_val) + heat_transfer_rhs
 
887
               
 
888
               ! Add to the internal energy equation's RHS
 
889
               call addto(rhs, internal_energy_nodes, rhs_addto) 
 
890
               ! Add to the matrix
 
891
               call addto(matrix, internal_energy_nodes, internal_energy_nodes, matrix_addto)
 
892
 
 
893
            end subroutine add_heat_transfer_element
 
894
 
 
895
      end subroutine add_heat_transfer
 
896
      
546
897
   end module multiphase_module