Add the dielectric heating term

* Add the dielectric heating term

The implementation add the RaT number which quantifies the heat added to the system according to the existing Rayleigh number.

* Add the dielectric heating term

The implementation add the RaT number which quantifies the heat added to the system according to the existing Rayleigh number.

Author: ygaillard

src/Namelists.f90

@@ -77,7 +77,7 @@ subroutine readNamelists(tscheme)
       &    mpi_transp,l_adv_curl,mpi_packing
 
       namelist/phys_param/                                        &
-      &    ra,rae,raxi,pr,sc,prmag,ek,gamma,epsc0,epscxi0,radratio,Bn,      &
+      &    ra,rae,rat,raxi,pr,sc,prmag,ek,gamma,epsc0,epscxi0,radratio,Bn,      &
       &    ktops,kbots,ktopv,kbotv,ktopb,kbotb,kbotxi,ktopxi,     &
       &    s_top,s_bot,impS,sCMB,xi_top,xi_bot,impXi,xiCMB,       &
       &    nVarCond,con_DecRate,con_RadRatio,con_LambdaMatch,     &
@@ -322,7 +322,9 @@ subroutine readNamelists(tscheme)
       l_AB1    =.false.
       l_bridge_step=.true.
       l_onset  =.false.
-      l_ehd_dep=.true.
+      l_ehd_dep=.false.
+      l_ehd_die=.false.
+
 
       if ( mode == 1 ) then
          !-- Only convection:
@@ -429,6 +431,13 @@ subroutine readNamelists(tscheme)
 
       if ( ra == 0.0_cp .and. rae == 0.0_cp ) l_heat=.false.
 
+      if ( rae > 0.0_cp ) then 
+         l_ehd_dep = .true.
+         if ( rat > 0.0_cp ) then
+            l_ehd_die = .true.
+         end if
+      end if
+
       if ( ek < 0.0_cp ) l_non_rot= .true.
       if ( l_non_rot ) then
          l_corr=.false.
@@ -929,6 +938,7 @@ subroutine writeNamelists(n_out)
       write(n_out,*) "&phys_param"
       write(n_out,'(''  ra              ='',ES14.6,'','')') ra
       write(n_out,'(''  rae             ='',ES14.6,'','')') rae
+      write(n_out,'(''  rat             ='',ES14.6,'','')') rat
       write(n_out,'(''  raxi            ='',ES14.6,'','')') raxi
       write(n_out,'(''  pr              ='',ES14.6,'','')') pr
       write(n_out,'(''  sc              ='',ES14.6,'','')') sc
@@ -1373,6 +1383,8 @@ subroutine defaultNamelists
 
       !----- Namelist phys_param:
       ra          =0.0_cp
+      rae         =0.0_cp
+      rat         =0.0_cp
       raxi        =0.0_cp
       ek          =1.0e-3_cp
       pr          =one

src/get_nl.f90

@@ -23,14 +23,14 @@ module grid_space_arrays_mod
    use radial_functions, only: or2, orho1, beta, otemp1, visc, r, or3, &
        &                       lambda, or4, or1
    use physical_parameters, only: radratio, LFfac, n_r_LCR, prec_angle, ViscHeatFac,    &
-        &                         oek, po, dilution_fac, ra, rae, gamma, opr, OhmLossFac, &
+        &                         oek, po, dilution_fac, ra, rae, rat, gamma, opr, OhmLossFac, &
         &                         epsPhase, phaseDiffFac, penaltyFac, tmelt
    use horizontal_data, only: sinTheta, cosTheta, phi, O_sin_theta_E2, &
        &                      cosn_theta_E2, O_sin_theta
    use parallel_mod, only: get_openmp_blocks
    use constants, only: two, third, one
    use logic, only: l_conv_nl, l_heat_nl, l_mag_nl, l_anel, l_mag_LF, l_adv_curl, &
-       &            l_chemical_conv, l_precession, l_centrifuge, l_phase_field, l_ehd_dep
+       &            l_chemical_conv, l_precession, l_centrifuge, l_phase_field, l_ehd_dep, l_ehd_die
 
    implicit none
 
@@ -456,6 +456,10 @@ subroutine get_nl(this, time, nR, nBc, lRmsCalc)
 
          end if  ! Viscous heating and Ohmic losses ?
 
+         if ( l_ehd_die ) then
+            this%heatTerms(:,nPhi)= &
+            &  opr * rae/rat * radratio**2/(1.0D0-radratio)**4 * or4(nR)
+         end if
       end do
       !$omp end parallel
 

src/get_td.f90

@@ -11,7 +11,7 @@ module nonlinear_lm_mod
    use truncation, only: lm_max, lm_maxMag, m_min
    use logic, only : l_anel, l_conv_nl, l_corr, l_heat, l_anelastic_liquid, &
        &             l_mag_nl, l_mag_kin, l_mag_LF, l_conv, l_mag,          &
-       &             l_chemical_conv, l_single_matrix, l_double_curl
+       &             l_chemical_conv, l_single_matrix, l_double_curl, l_ehd_die
    use radial_functions, only: r, or2, or1, beta, epscProf, or4, temp0, orho1, l_R
    use physical_parameters, only: CorFac, epsc,  n_r_LCR, epscXi
    use blocking, only: lm2l, lm2lmA, lm2lmS
@@ -60,7 +60,7 @@ subroutine initialize(this,lmP_max)
       allocate( this%VxBrLM(lmP_max), this%VxBtLM(lmP_max), this%VxBpLM(lmP_max))
       bytes_allocated = bytes_allocated + 6*lmP_max*SIZEOF_DEF_COMPLEX
 
-      if ( l_anel) then
+      if ( l_anel .or. l_ehd_die ) then
          allocate( this%heatTermsLM(lmP_max) )
          bytes_allocated = bytes_allocated+lmP_max*SIZEOF_DEF_COMPLEX
       end if
@@ -120,7 +120,7 @@ subroutine set_zero(this)
             this%VStLM(lm)=zero
             this%VSpLM(lm)=zero
          end if
-         if ( l_anel ) this%heatTermsLM(lm)=zero
+         if ( l_anel .or. l_ehd_die ) this%heatTermsLM(lm)=zero
          if ( l_chemical_conv ) then
             this%VXitLM(lm)=zero
             this%VXipLM(lm)=zero
@@ -478,15 +478,15 @@ subroutine get_dsdt(this, nR, nBc, dsdt, dVSrLM)
       if ( nBc == 0 ) then
 
          dsdt(1)=epsc*epscProf(nR)!+opr/epsS*divKtemp0(nR)
-         if ( l_anel ) then
+         if ( l_anel .or. l_ehd_die ) then
             if ( l_anelastic_liquid ) then
                dsdt(1)=dsdt(1)+temp0(nR)*this%heatTermsLM(1)
             else
                dsdt(1)=dsdt(1)+this%heatTermsLM(1)
             end if
          end if
 
-         if ( l_anel ) then
+         if ( l_anel .or. l_ehd_die ) then
             if ( l_anelastic_liquid ) then
                !$omp parallel do
                do lm=lm_min,lm_max

src/logic.f90

@@ -73,6 +73,7 @@ module logic
    logical :: l_LCR          ! Switch for zero electrical conductivity beyond r_LCR
    logical :: lVerbose       ! Switch for detailed information about run progress
    logical :: l_ehd_dep      ! Switch for dilectrophoretic force
+   logical :: l_ehd_die      ! Switch for dielectric heating
 
    logical :: l_PressGraph   ! Store pressure in graphic files
 

src/phys_param.f90

@@ -38,6 +38,7 @@ module physical_parameters
    real(cp) :: radratio       ! aspect ratio
    real(cp) :: ra             ! Rayleigh number
    real(cp) :: rae            ! Electric Rayleigh number
+   real(cp) :: rat            ! Thermal heating from dielectric heat
    real(cp) :: raxi           ! Chemical composition-based Rayleigh number
    real(cp) :: sc             ! Schmidt number (i.e. chemical Prandtl number)
    real(cp) :: ek             ! Ekman number

src/rIter.f90

@@ -21,7 +21,7 @@ module rIter_mod
        &            l_precession, l_centrifuge, l_adv_curl,          &
        &            l_double_curl, l_parallel_solve, l_single_matrix,&
        &            l_temperature_diff, l_RMS, l_phase_field,        &
-       &            l_onset, l_DTrMagSpec, l_ehd_dep
+       &            l_onset, l_DTrMagSpec, l_ehd_dep, l_ehd_die
    use radial_data, only: n_r_cmb, n_r_icb, nRstart, nRstop, nRstartMag, &
        &                  nRstopMag
    use radial_functions, only: or2, orho1, l_R
@@ -683,7 +683,7 @@ subroutine transform_to_lm_space(this, nR, lRmsCalc, dVSrLM, dVXirLM, dphidt)
          call spat_to_qst(this%gsa%VSr, this%gsa%VSt, this%gsa%VSp, &
               &           dVSrLM, this%nl_lm%VStLM, this%nl_lm%VSpLM, l_R(nR))
 
-         if ( l_anel ) call scal_to_SH(this%gsa%heatTerms, this%nl_lm%heatTermsLM, &
+         if ( l_anel .or. l_ehd_die ) call scal_to_SH(this%gsa%heatTerms, this%nl_lm%heatTermsLM, &
                             &          l_R(nR))
       end if
       if ( l_chemical_conv ) then