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,:)
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)
542
581
end subroutine add_fluid_particle_drag_element
544
583
end subroutine add_fluid_particle_drag
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).
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
601
type(element_type) :: test_function
602
type(element_type), pointer :: internal_energy_shape
603
integer, dimension(:), pointer :: internal_energy_nodes
606
type(vector_field), pointer :: x, velocity_fluid
610
logical :: is_particle_phase
611
logical :: not_found ! Error flag. Have we found the fluid phase?
612
integer :: i, istate_fluid, istate_particle
615
ewrite(1, *) "Entering add_heat_transfer"
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", &
622
! Get the coordinate field from state(istate)
623
x => extract_vector_field(state(istate), "Coordinate")
625
assert(x%dim == mesh_dim(internal_energy))
626
assert(ele_count(x) == ele_count(internal_energy))
628
! Are we using a discontinuous Galerkin discretisation?
629
dg = continuity(internal_energy) < 0
631
! Is this phase a particle phase?
632
is_particle_phase = have_option(trim(state(istate)%option_path)//"/multiphase_properties/particle_diameter")
634
! Retrieve the index of the fluid phase in the state array.
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
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
645
if(.not.not_found) then
646
FLExit("Heat transfer term does not currently support more than one fluid phase.")
654
istate_fluid = istate
659
FLExit("No fluid phase found for the heat transfer term.")
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)
666
state_loop: do i = 1, size(state)
667
if(i /= istate_fluid) then
668
call assemble_heat_transfer(istate_fluid, i)
673
ewrite(1, *) "Exiting add_heat_transfer"
677
subroutine assemble_heat_transfer(istate_fluid, istate_particle)
679
integer, intent(in) :: istate_fluid, istate_particle
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
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
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")
706
call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter", d)
708
call get_option(trim(state(istate_fluid)%option_path)//"/multiphase_properties/effective_conductivity", k, kstat)
710
call get_option(trim(state(istate_fluid)%option_path)//"/multiphase_properties/specific_heat", C_fluid, cstat_fluid)
712
call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/specific_heat", C_particle, cstat_particle)
714
call get_option(trim(state(istate_fluid)%option_path)//"/equation_of_state/compressible/stiffened_gas/ratio_specific_heats", gamma, gstat)
717
FLExit("For inter-phase heat transfer, an effective_conductivity needs to be specified for the fluid phase.")
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.")
723
FLExit("For inter-phase heat transfer, ratio_specific_heats needs to be specified for the compressible phase.")
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)
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")
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")
744
! ----- Volume integrals over elements -------------
745
call profiler_tic(internal_energy, "element_loop")
746
element_loop: do ele = 1, element_count(internal_energy)
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
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, &
767
call profiler_toc(internal_energy, "element_loop")
769
call deallocate(nvfrac_fluid)
770
call deallocate(nvfrac_particle)
773
end subroutine assemble_heat_transfer
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, &
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
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
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
810
real, dimension(ele_loc(internal_energy,ele)) :: old_internal_energy_val
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
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|
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
829
call transform_to_physical(x, ele, detwei)
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)
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))
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,:)
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
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))
859
Q = (6.0*k*vfrac_particle_gi*particle_Nu)/(d**2)
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)
868
coefficient_for_matrix = -Q/C_fluid
869
coefficient_for_rhs = Q*(internal_energy_particle_gi/C_particle)
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)
879
if(is_particle_phase) then
880
old_internal_energy_val = ele_val(old_internal_energy_particle, ele)
882
old_internal_energy_val = ele_val(old_internal_energy_fluid, ele)
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
888
! Add to the internal energy equation's RHS
889
call addto(rhs, internal_energy_nodes, rhs_addto)
891
call addto(matrix, internal_energy_nodes, internal_energy_nodes, matrix_addto)
893
end subroutine add_heat_transfer_element
895
end subroutine add_heat_transfer
546
897
end module multiphase_module