Add EHD force to m_coupling

src/definitions.make

@@ -26,6 +26,8 @@ m_chemistry.o: m_table_data.mod
 m_chemistry.o: m_transport_data.mod
 m_chemistry.o: m_types.mod
 m_chemistry.o: m_units_constants.mod
+m_coupling.o: m_chemistry.mod
+m_coupling.o: m_field.mod
 m_coupling.o: m_gas.mod
 m_coupling.o: m_streamer.mod
 m_coupling.o: m_types.mod
@@ -158,3 +160,4 @@ streamer.o: m_types.mod
 streamer.o: m_units_constants.mod
 streamer.o: m_user_methods.mod
 streamer.o: m_user.mod
+

src/m_coupling.f90

@@ -6,6 +6,8 @@ module m_coupling
   use m_units_constants
   use m_gas
   use m_streamer
+  use m_chemistry
+  use m_field
 
   implicit none
   private
@@ -27,7 +29,8 @@ subroutine add_heating_box(box, dt_vec)
     type(box_t), intent(inout) :: box
     real(dp), intent(in)       :: dt_vec(:)
     integer                    :: IJK, nc
-    real(dp)                   :: J_dot_E
+    real(dp)                   :: J_dot_E, ehd_force(NDIM)
+    real(dp)                   :: E_vector(DTIMES(1:box%n_cell), NDIM)
 
     nc = box%n_cell
     do KJI_DO(1, nc)
@@ -48,8 +51,20 @@ subroutine add_heating_box(box, dt_vec)
 #endif
 
        box%cc(IJK, gas_vars(i_e)) = box%cc(IJK, gas_vars(i_e)) + &
-           gas_heating_efficiency *  J_dot_E * UC_elec_charge * dt_vec(1)
+            gas_heating_efficiency *  J_dot_E * UC_elec_charge * dt_vec(1)
     end do; CLOSE_DO
+
+    E_vector = field_get_E_vector(box)
+
+    ! EHD force
+    do KJI_DO(1, nc)
+       ehd_force = UC_elem_charge * sum(box%cc(IJK, charged_species_itree) * &
+            charged_species_charge) * E_vector(IJK, :)
+
+       box%cc(IJK, gas_vars(i_mom)) = box%cc(IJK, gas_vars(i_mom)) + &
+            gas_EHD_factor * ehd_force * dt_vec(1)
+    end do; CLOSE_DO
+
   end subroutine add_heating_box
 
   !> Update gas number density in the fluid model

src/m_field.f90

@@ -77,6 +77,7 @@ module m_field
   public :: field_bc_homogeneous
   public :: field_from_potential
   public :: field_compute_energy
+  public :: field_get_E_vector
 
 contains
 
@@ -617,4 +618,29 @@ real(dp) function reduce_sum(a, b)
     reduce_sum = a + b
   end function reduce_sum
 
+  function field_get_E_vector(box) result(E_vector)
+    type(box_t), intent(in) :: box
+    real(dp)                :: E_vector(DTIMES(1:box%n_cell), NDIM)
+    integer                 :: nc
+
+    nc = box%n_cell
+
+#if NDIM == 1
+    E_vector(DTIMES(:), 1) = 0.5_dp * (box%fc(1:nc, 1, electric_fld) + &
+         box%fc(2:nc+1, 1, electric_fld))
+#elif NDIM == 2
+    E_vector(DTIMES(:), 1) = 0.5_dp * (box%fc(1:nc, 1:nc, 1, electric_fld) + &
+         box%fc(2:nc+1, 1:nc, 1, electric_fld))
+    E_vector(DTIMES(:), 2) = 0.5_dp * (box%fc(1:nc, 1:nc, 2, electric_fld) + &
+         box%fc(1:nc, 2:nc+1, 2, electric_fld))
+#elif NDIM == 3
+    E_vector(DTIMES(:), 1) = 0.5_dp * (box%fc(1:nc, 1:nc, 1:nc, 1, electric_fld) + &
+         box%fc(2:nc+1, 1:nc, 1:nc, 1, electric_fld))
+    E_vector(DTIMES(:), 2) = 0.5_dp * (box%fc(1:nc, 1:nc, 1:nc, 2, electric_fld) + &
+         box%fc(1:nc, 2:nc+1, 1:nc, 2, electric_fld))
+    E_vector(DTIMES(:), 3) = 0.5_dp * (box%fc(1:nc, 1:nc, 1:nc, 3, electric_fld) + &
+         box%fc(1:nc, 1:nc, 2:nc+1, 3, electric_fld))
+#endif
+  end function field_get_E_vector
+
 end module m_field

src/m_gas.f90

@@ -50,6 +50,9 @@ module m_gas
   ! Joule heating efficiency
   real(dp), public, protected :: gas_heating_efficiency  = 1.0_dp
 
+  ! Factor for the EHD force term (should be 1 by default)
+  real(dp), public, protected :: gas_EHD_factor = 1.0_dp
+
   ! Ratio of heat capacities (polytropic index)
   real(dp), public, protected :: gas_euler_gamma = 1.4_dp
 
@@ -147,6 +150,8 @@ subroutine gas_initialize(tree, cfg)
          "Gas mean molecular weight (kg), for gas dynamics")
     call CFG_add_get(cfg, "gas%heating_efficiency", gas_heating_efficiency, &
          "Joule heating efficiency (between 0.0 and 1.0)")
+    call CFG_add_get(cfg, "gas%EHD_factor", gas_EHD_factor, &
+         "Factor for the EHD force term (should be 1 by default)")
 
     ! Ideal gas law
     gas_number_density = 1e5_dp * gas_pressure / &

src/m_output.f90

@@ -7,6 +7,7 @@ module m_output
   use m_gas
   use m_lookup_table
   use m_units_constants
+  use m_field
 
   implicit none
   private
@@ -452,7 +453,7 @@ subroutine add_variables(box, new_vars, n_var)
                N_inv * box%cc(DTIMES(:), i_electron) * UC_elem_charge
        case ("Je_1")
           ! Add electron current density  (e mu n_e E_vector)
-          E_vector = get_E_vector(box)
+          E_vector = field_get_E_vector(box)
           sigma = LT_get_col(td_tbl, td_mobility, &
                Td(DTIMES(1:nc))) * N_inv * box%cc(DTIMES(1:nc), i_electron) * &
                UC_elem_charge
@@ -963,31 +964,6 @@ subroutine set_power_density_box(box)
     end do; CLOSE_DO
   end subroutine set_power_density_box
 
-  function get_E_vector(box) result(E_vector)
-    type(box_t), intent(in) :: box
-    real(dp)                :: E_vector(DTIMES(1:box%n_cell), NDIM)
-    integer                 :: nc
-
-    nc = box%n_cell
-
-#if NDIM == 1
-    E_vector(DTIMES(:), 1) = 0.5_dp * (box%fc(1:nc, 1, electric_fld) + &
-         box%fc(2:nc+1, 1, electric_fld))
-#elif NDIM == 2
-    E_vector(DTIMES(:), 1) = 0.5_dp * (box%fc(1:nc, 1:nc, 1, electric_fld) + &
-         box%fc(2:nc+1, 1:nc, 1, electric_fld))
-    E_vector(DTIMES(:), 2) = 0.5_dp * (box%fc(1:nc, 1:nc, 2, electric_fld) + &
-         box%fc(1:nc, 2:nc+1, 2, electric_fld))
-#elif NDIM == 3
-    E_vector(DTIMES(:), 1) = 0.5_dp * (box%fc(1:nc, 1:nc, 1:nc, 1, electric_fld) + &
-         box%fc(2:nc+1, 1:nc, 1:nc, 1, electric_fld))
-    E_vector(DTIMES(:), 2) = 0.5_dp * (box%fc(1:nc, 1:nc, 1:nc, 2, electric_fld) + &
-         box%fc(1:nc, 2:nc+1, 1:nc, 2, electric_fld))
-    E_vector(DTIMES(:), 3) = 0.5_dp * (box%fc(1:nc, 1:nc, 1:nc, 3, electric_fld) + &
-         box%fc(1:nc, 1:nc, 2:nc+1, 3, electric_fld))
-#endif
-  end function get_E_vector
-
   subroutine set_gas_primitives_box(box)
     type(box_t), intent(inout) :: box
     integer                    :: IJK, nc, idim