!=============================================================================== ! User source terms definition. ! ! 1) Momentum equation (coupled solver) ! 2) Species transport ! 3) Turbulence (k-epsilon, k-omega, Rij-epsilon, v2-f, Spalart-Allmaras) !=============================================================================== !------------------------------------------------------------------------------- ! Code_Saturne version 5.0.3 ! -------------------------- ! This file is part of Code_Saturne, a general-purpose CFD tool. ! ! Copyright (C) 1998-2017 EDF S.A. ! ! This program is free software; you can redistribute it and/or modify it under ! the terms of the GNU General Public License as published by the Free Software ! Foundation; either version 2 of the License, or (at your option) any later ! version. ! ! This program is distributed in the hope that it will be useful, but WITHOUT ! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS ! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more ! details. ! ! You should have received a copy of the GNU General Public License along with ! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin ! Street, Fifth Floor, Boston, MA 02110-1301, USA. !------------------------------------------------------------------------------- !=============================================================================== ! Purpose: ! ------- !> \file cs_user_source_terms.f90 !> !> \brief Additional right-hand side source terms !> !> See \subpage cs_user_source_terms and \subpage cs_user_source_terms-scalar_in_a_channel !> for examples. !> !> \brief Additional right-hand side source terms for velocity components equation !> (Navier-Stokes) !> !> \section ustsnv_use Usage !> !> The additional source term is decomposed into an explicit part (\c crvexp) and !> an implicit part (\c crvimp) that must be provided here. !> The resulting equation solved by the code for a velocity is: !> \f[ !> \rho \norm{\vol{\celli}} \DP{\vect{u}} + .... !> = \tens{crvimp} \cdot \vect{u} + \vect{crvexp} !> \f] !> !> Note that \c crvexp and \c crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in kg.m/s2 !> - crvimp is expressed in kg/s !> !> The \c crvexp and \c crvimp arrays are already initialized to 0 !> before entering the !> the routine. It is not needed to do it in the routine (waste of CPU time). !> !> \remark The additional force on \f$ x_i \f$ direction is given by !> \c crvexp(i, iel) + vel(j, iel)* crvimp(j, i). !> !> For stability reasons, Code_Saturne will not add -crvimp directly to the !> diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp term is !> treated implicitely only if it strengthens the diagonal of the matrix. !> However, when using the second-order in time scheme, this limitation cannot !> be done anymore and -crvimp is added directly. The user should therefore test !> the negativity of crvimp by himself. !> !> When using the second-order in time scheme, one should supply: !> - crvexp at time n !> - crvimp at time n+1/2 !> !> The selection of cells where to apply the source terms is based on a !> \ref getcel command. For more info on the syntax of the \ref getcel command, !> refer to the user manual or to the comments on the similar command !> \ref getfbr in the routine \ref cs_user_boundary_conditions. ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! Arguments !______________________________________________________________________________. ! mode name role ! !______________________________________________________________________________! !> \param[in] nvar total number of variables !> \param[in] nscal total number of scalars !> \param[in] ncepdp number of cells with head loss terms !> \param[in] ncesmp number of cells with mass source terms !> \param[in] ivar index number of the current variable !> \param[in] icepdc index number of cells with head loss terms !> \param[in] icetsm index number of cells with mass source terms !> \param[in] itypsm type of mass source term for each variable !> (see \ref cs_user_mass_source_terms) !> \param[in] dt time step (per cell) !> \param[in] ckupdc head loss coefficient !> \param[in] smacel value associated to each variable in the mass !> source terms or mass rate (see !> \ref cs_user_mass_source_terms) !> \param[out] crvexp explicit part of the source term !> \param[out] crvimp implicit part of the source term !_______________________________________________________________________________ subroutine ustsnv & ( nvar , nscal , ncepdp , ncesmp , & ivar , & icepdc , icetsm , itypsm , & dt , & ckupdc , smacel , & crvexp , crvimp ) !=============================================================================== !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use entsor use optcal use cstphy use parall use period use mesh use field use cs_c_bindings !=============================================================================== implicit none ! Arguments integer nvar , nscal integer ncepdp , ncesmp integer ivar integer icepdc(ncepdp) integer icetsm(ncesmp), itypsm(ncesmp,nvar) double precision dt(ncelet) double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar) double precision crvexp(3,ncelet), crvimp(3,3,ncelet) ! Local variables character(len=80) :: chaine,name integer, allocatable, dimension(:) :: lstelt integer iel, ilelt, nlelt DOUBLE PRECISION tho, uio type(var_cal_opt) :: vcopt !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! Allocate a temporary array for cells selection !write(*,'(A16,I6)')'ustsnv start, t=',ntcabs allocate(lstelt(ncelet)) ! temporary array for cells selection ! --- Get variable calculation options call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) if (vcopt%iwarni.ge.1) then call field_get_label(ivarfl(ivar), chaine) write(nfecra,1000) chaine(1:8) endif tho = 1.d-30 uio = 0.d0 !call field_get_name(ivarfl(ivar),name) !call field_get_label(ivarfl(ivar), chaine) !if (ntcabs.eq.1) write(*,1001) ntcabs,ivar,name,chaine call getcel('domaine_cathode',nlelt,lstelt) do ilelt = 1, nlelt iel = lstelt(ilelt) crvimp(1, 1, iel) = (-1.0)/ tho crvimp(2, 2, iel) = (-1.0)/ tho crvimp(3, 3, iel) = (-1.0)/ tho crvexp(1, iel) = uio/tho crvexp(2, iel) = uio/tho crvexp(3, iel) = uio/tho enddo call getcel('anode',nlelt,lstelt) do ilelt = 1, nlelt iel = lstelt(ilelt) crvimp(1, 1, iel) = (-1.0)/ tho crvimp(2, 2, iel) = (-1.0)/ tho crvimp(3, 3, iel) = (-1.0)/ tho crvexp(1, iel) = uio/tho crvexp(2, iel) = uio/tho crvexp(3, iel) = uio/tho enddo !-------- ! Formats !-------- 1000 format(' User source terms for variable ',A8,/) 1001 format('ustsnv t=',I3,' Name of ivarfl(',I3,')=',A22,' label=',A22) !---- ! End !---- ! Deallocate the temporary array deallocate(lstelt) return end subroutine ustsnv !=============================================================================== !=============================================================================== ! Purpose: ! ------- !> User subroutine. !> \brief Additional right-hand side source terms for scalar equations (user !> scalars and specific physics scalars). !> !> !> Usage !> ----- !> The routine is called for each scalar, user or specific physisc. It is !> therefore necessary to test the value of the scalar number iscal to separate !> the treatments of the different scalars (if (iscal.eq.p) then ....). !> !> The additional source term is decomposed into an explicit part (crvexp) and !> an implicit part (crvimp) that must be provided here. !> The resulting equation solved by the code for a scalar f is: !> !> \f[ \rho*volume*\frac{df}{dt} + .... = crvimp*f + crvexp \f] !> !> !> Note that crvexp and crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in kg.[scal]/s, where [scal] is the unit of the scalar !> - crvimp is expressed in kg/s !> !> !> The crvexp and crvimp arrays are already initialized to 0 before entering the !> the routine. It is not needed to do it in the routine (waste of CPU time). !> !> For stability reasons, Code_Saturne will not add -crvimp directly to the !> diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp term is !> treated implicitely only if it strengthens the diagonal of the matrix. !> However, when using the second-order in time scheme, this limitation cannot !> be done anymore and -crvimp is added directly. The user should therefore test !> the negativity of crvimp by himself. !> !> When using the second-order in time scheme, one should supply: !> - crvexp at time n !> - crvimp at time n+1/2 !> !> !> The selection of cells where to apply the source terms is based on a getcel !> command. For more info on the syntax of the getcel command, refer to the !> user manual or to the comments on the similar command \ref getfbr in the routine !> \ref cs_user_boundary_conditions. !> WARNING: If scalar is the temperature, the resulting equation !> solved by the code is: !> !> rho*Cp*volume*dT/dt + .... = crvimp*T + crvexp !> !> !> Note that crvexp and crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in W !> - crvimp is expressed in W/K !> !> !> STEP SOURCE TERMS !>=================== !> In case of a complex, non-linear source term, say F(f), for scalar f, the !> easiest method is to implement the source term explicitely. !> !> df/dt = .... + F(f(n)) !> where f(n) is the value of f at time tn, the beginning of the time step. !> !> This yields : !> crvexp = volume*F(f(n)) !> crvimp = 0 !> !> However, if the source term is potentially steep, this fully explicit !> method will probably generate instabilities. It is therefore wiser to !> partially implicit the term by writing: !> !> df/dt = .... + dF/df*f(n+1) - dF/df*f(n) + F(f(n)) !> !> This yields: !> crvexp = volume*( F(f(n)) - dF/df*f(n) ) !> crvimp = volume*dF/df !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! Arguments !______________________________________________________________________________. ! mode name role !______________________________________________________________________________! !> \param[in] nvar total number of variables !> \param[in] nscal total number of scalars !> \param[in] ncepdp number of cells with head loss terms !> \param[in] ncssmp number of cells with mass source terms !> \param[in] iscal index number of the current scalar !> \param[in] icepdc index number of cells with head loss terms !> \param[in] icetsm index number of cells with mass source terms !> \param[in] itypsm type of mass source term for each variable !> \param[in] (see \ref cs_user_mass_source_terms) !> \param[in] dt time step (per cell) !> \param[in] ckupdc head loss coefficient !> \param[in] smacel value associated to each variable in the mass !> \param[in] source terms or mass rate (see \ref cs_user_mass_source_terms) !> \param[out] crvexp explicit part of the source term !> \param[out] crvimp implicit part of the source term !______________________________________________________________________________! subroutine ustssc & ( nvar , nscal , ncepdp , ncesmp , & iscal , & icepdc , icetsm , itypsm , & dt , & ckupdc , smacel , & crvexp , crvimp ) !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use entsor use optcal use cstphy use parall use period use mesh use field use ppincl !=============================================================================== implicit none ! Arguments integer nvar , nscal integer ncepdp , ncesmp integer iscal integer icepdc(ncepdp) integer icetsm(ncesmp), itypsm(ncesmp,nvar) double precision dt(ncelet) double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar) double precision crvexp(ncelet), crvimp(ncelet) ! Local variables character(len=80) :: label integer, allocatable, dimension(:) :: lstelt integer ivar, iel, nlelt, timetestvar integer ilelt, ifac integer ii, jj, iels, ielp,ielhf, MODE integer ifcp,ifcvsl,ifctemp,ifcur!,ifcur2,ifcur3 !pointers integer ifdomtyp,ifnbfrfl,ifintfc,ifnb_htfl !pointers integer nlelt_cath,nlelt_anode,nlelt_fluide,nlelt_chand,nlelt_chcat,source_counter integer, allocatable, dimension(:) :: lstelt_cath,lstelt_anode integer, allocatable, dimension(:) :: lstelt_fluide,lstelt_chand,lstelt_chcat double precision LambdaIfac double precision PROPp, PROPs, ElecSource, DiffSource double precision Q, J,Jp,Js, Jold, kb, eee, Ua, phimat double precision QEM, QI, JEM, JI, h, Me, Uk, phii double precision ARichLaw, LRichLaw, kin_enrg_coef double precision pii, TP , TW, TWmax , TE, cor_tw_coef double precision LamFac, surface, SUMM double precision Ponds, distloc, pondloc, surfloc,temp,source_sum double precision rayo double precision ZZ1 , ZZ2 , ZZ3 , ZZT double precision ZZ4 , ZZ5 , ZZ6 , ZZ7 , ZZ8, ZZ10, ZZ9 double precision ZZA1, ZZA2, ZZA3, ZZAT, TP2, ZZA13, ZZA14 double precision ZZA4, ZZA5, ZZA6, ZZA7, ZZA8, ZZA9, ZZA10 !dist and pond are symbols from module 'mesh' !double precision, dimension(:), pointer :: cpro_viscl, cpro_rom, cpro_sig, cpro_crom !double precision, dimension(:), pointer :: cpro_vpoti, cpro_vpotr double precision, dimension(:), pointer :: cvar_hm,cpro_vscalt,cpro_cp,cpro_temp double precision, dimension(:), pointer :: cpro_domtyp,cpro_nbfrfl,cpro_intfc,cpro_nb_htfl double precision, dimension(:,:), pointer :: cpro_cur !,cpro_cur2,cpro_cur3 character(len=80) :: name !integer cellsnum,chechinit1,chechinit2 !parameter (cellsnum = 93375)!ncelet) !integer inb_from_fl(cellsnum),intrf_face(cellsnum),domtypeloc(cellsnum) !common /geomforsourcei/ chechinit1,chechinit2,inb_from_fl,intrf_face,domtypeloc !=============================================================================== !=============================================================================== ! ENTHALPY SOURCE CALCULATION AT THE INTERFACE OF ELECTRODES !=============================================================================== !if (cellsnum.lt.ncelet) print*,'USTSSC ERROR! Wrong number of cells in geometrical parameters' !POINTERS call field_get_val_s(ivarfl(isca(iscalt)), cvar_hm) !ENTHALPY. Obtaining the field by the pointer. Actually id of enthalpy field is 2. call field_get_id_try('thermal_diffusivity', ifcvsl) !LAMBDA. Dynamic difusivity of enthalpy, it is the ratio between thermal conductivity and specific heat in the enthalpy equation. !call field_get_key_int(ivarfl(isca(iscalt)), kivisl, ifcvsl) !if (ntcabs-ntpabs.eq.1) then ! write(*,'(A,I4,A,I4)')' iscalt=',iscalt,' ifcvsl',ifcvsl ! call field_get_name(ivarfl(isca(iscalt)),name) ! write(*,'(A,A)')'USTSSC name of iscalt=',name ! call field_get_name(ifcvsl,name) ! write(*,'(A,A)')'USTSSC name of ifcvsl=',name !end if if (ifcvsl.ge.0) call field_get_val_s(ifcvsl, cpro_vscalt) if (ifcvsl.lt.0) print*,'Error: cs_user_source_terms.f90 No id of Dynamic difusivity of enthalpy' call field_get_id_try('temperature', ifctemp) !TEMPERATURE. Obtaining pointer to the field. Actually it is 14 if (ifctemp.ge.0) call field_get_val_s(ifctemp, cpro_temp) if (ifctemp.lt.0) print*,'Error: cs_user_source_terms.f90 No id of temperature' call field_get_id_try('specific_heat', ifcp) !SPECIFIC HEAT. Obtaining pointer to Cp. icp turns out to be current_re_1 if (ifcp.ge.0) call field_get_val_s(ifcp, cpro_cp) !obtaining Cp field if (ifcp.lt.0) print*,'Error: cs_user_source_terms.f90 No id of specific heat' call field_get_id_try('domain_type', ifdomtyp) !Domain type. if (ifdomtyp.ge.0) call field_get_val_s(ifdomtyp, cpro_domtyp) if (ifdomtyp.lt.0) print*,'Error: cs_user_source_terms.f90 No id of domain_type' call field_get_id_try('inb_from_fl', ifnbfrfl) !Neighboring cell from fluid for the cell of metal at the interface if (ifnbfrfl.ge.0) call field_get_val_s(ifnbfrfl, cpro_nbfrfl) if (ifnbfrfl.lt.0) print*,'Error: cs_user_source_terms.f90 No id of inb_from_fl' call field_get_id_try('intrf_face', ifintfc) !Interface face between metal and fluid cells if (ifintfc.ge.0) call field_get_val_s(ifintfc, cpro_intfc) if (ifintfc.lt.0) print*,'Error: cs_user_source_terms.f90 No id of intrf_face' call field_get_id_try('inb_heatflux', ifnb_htfl) !Interface face between metal and fluid cells if (ifnb_htfl.ge.0) call field_get_val_s(ifnb_htfl, cpro_nb_htfl) if (ifnb_htfl.lt.0) print*,'Error: cs_user_source_terms.f90 No id of inb_heatflux' call field_get_id_try('current_re', ifcur) !jx = sig * Ex imaginary electric current density [A/m^2] x-component if (ifcur.ge.0) call field_get_val_v(ifcur, cpro_cur) if (ifcur.lt.0) print*,'Error: cs_user_source_terms.f90 No id of current density' !call field_get_id_try('current_re_2', ifcur2) !jy - current density y-component !ifcur2 = ifcur1+1 !if (ifcur2.ge.0) call field_get_val_s(ifcur2, cpro_cur2) !if (ifcur2.lt.0) print*,'Error: cs_user_source_terms.f90 No id of current density y-component' !call field_get_id_try('current_re_3', ifcur3) !jz - current density z-component !ifcur3 = ifcur2+1 !if (ifcur3.ge.0) call field_get_val_s(ifcur3, cpro_cur3) !if (ifcur3.lt.0) print*,'Error: cs_user_source_terms.f90 No id of current density z-component' !CONST DATA kb = 1.38065d-23 !boltzmann constant eee = 1.60218d-19 !electron charge pii = 3.14159265359d0!π number Me = 9.10938188d-31 !electron rest mass h = 6.62607d-34 !planck constant Uk = 1.0d1 !cathode voltage drop Ua = 4.0d0 !anode voltage drop phimat= 4.2d0 !work function Фc = 4.2 eV phii = 1.36d1 !ionization energy ARichLaw = eee*4d0*pii*kb*kb*Me/(h*h*h) !eee*4*pii*kb*kb*Me/(h*h*h)=1.2*10^(-19-23-23-31+34*3)=1.2*10^6 LRichLaw=-1d0*eee*phimat/kb kin_enrg_coef=5d0*kb/(2d0*eee) if(iscal.eq.iscalt .and. ifdomtyp.ge.0 .and. ifnbfrfl.ge.0 .and. ifintfc.ge.0) then ! iscalt - number of the thermal variable, enthalpy for electric arcs !STEP 1. Initialization of inb_from_fl and intrf_face with 0 for entire domain. !ntcabs - current time step !ntpabs - number of the last time step in the previous calculation !ntcabs-ntpabs - time step number in current calculation if (ntcabs-ntpabs.eq.1) then ! .or. chechinit1.ne.666001 write(*,'(A,I5)')'USTSSC Initialization of interface arrays with 0 t=',ntcabs !,' chechinit1=',chechinit1 !chechinit1=666001 ! Allocate temporary arrays for cells selection allocate(lstelt_cath(ncelet)) allocate(lstelt_anode(ncelet)) allocate(lstelt_fluide(ncelet)) allocate(lstelt_chand(ncelet)) allocate(lstelt_chcat(ncelet)) allocate(lstelt(ncelet)) !temporary array for cells selection !Initialization of inb_from_fl and intrf_face with 0 for entire domain. call getcel('all[]',nlelt,lstelt) write(*,'(A,I6,A,I6,A,I6,A,I6)')'all[] nlelt=',nlelt,' ncel=',ncel,' ncelet=',ncelet,' nfac=',nfac do ilelt = 1,nlelt iel = lstelt(ilelt) cpro_nbfrfl(iel)=0.0 cpro_intfc(iel)=0.0 cpro_nb_htfl(iel)=0.0 end do !Initialization of cells types by groups do iel = 1, ncel !domtypeloc(iel) = 0 cpro_domtyp(iel)= 0 end do call getcel('domaine_cathode',nlelt_cath,lstelt_cath) do ilelt = 1,nlelt_cath iel = lstelt_cath(ilelt) !domtypeloc(iel) = 1 cpro_domtyp(iel)= 1 enddo call getcel('anode',nlelt_anode,lstelt_anode) do ilelt = 1,nlelt_anode iel = lstelt_anode(ilelt) !domtypeloc(iel) = 2 cpro_domtyp(iel)= 2 enddo call getcel('domaine_fluide',nlelt_fluide,lstelt_fluide) do ilelt = 1,nlelt_fluide iel = lstelt_fluide(ilelt) !domtypeloc(iel) = 0 cpro_domtyp(iel)= 0 enddo call getcel('chute_anodique',nlelt_chand,lstelt_chand) do ilelt = 1,nlelt_chand iel = lstelt_chand(ilelt) !domtypeloc(iel) = 3 cpro_domtyp(iel)= 3 enddo call getcel('chute_cathodique',nlelt_chcat,lstelt_chcat) do ilelt = 1,nlelt_chcat iel = lstelt_chcat(ilelt) !domtypeloc(iel) = 4 cpro_domtyp(iel)= 4 enddo write(*,'(I3,A)')irangp,' MESH INFORMATION:' write(*,'(I3,A,I6)')irangp,' domaine_cathode Ncells=',nlelt_cath write(*,'(I3,A,I6)')irangp,' anode Ncells=',nlelt_anode write(*,'(I3,A,I6)')irangp,' domaine_fluide Ncells=',nlelt_fluide write(*,'(I3,A,I6)')irangp,' chute_anodique Ncells=',nlelt_chand write(*,'(I3,A,I6)')irangp,' chute_cathodiqu Ncells=',nlelt_chcat ! Deallocate the temporary array deallocate(lstelt_cath) deallocate(lstelt_anode) deallocate(lstelt_fluide) deallocate(lstelt_chand) deallocate(lstelt_chcat) deallocate(lstelt) end if !STEP 2. Search for interface faces. if (ntcabs-ntpabs.eq.2) then !.or. (chechinit1.eq.666001 .and. chechinit2.ne.666002) !chechinit2=666002 !if time step in current calculation is 1, the next one will be 2 anyway. That's why we don't need chechinit2 from now. write(*,'(A,I5)')'USTSSC Calculation of interface arrays. t=',ntcabs do ifac = 1, nfac ii = ifacel(1,ifac) !IFACEL - Index-numbers (id est IEL) of the two (only) neighboring cells for each internal face. jj = ifacel(2,ifac) !positive cpro_nbfrfl element is the number of neighboring cell from fluid !negative cpro_nbfrfl element is the number of neighboring cell from the electrode if (cpro_domtyp(ii).eq.1 .and. cpro_domtyp(jj).eq.4) then !domaine_cathode|chute_cathodique cpro_nbfrfl(ii)=jj cpro_nbfrfl(jj)=-ii cpro_intfc(ii)=ifac else if (cpro_domtyp(ii).eq.4 .and. cpro_domtyp(jj).eq.1) then !chute_cathodique|domaine_cathode cpro_nbfrfl(ii)=-jj cpro_nbfrfl(jj)=ii cpro_intfc(jj)=ifac else if (cpro_domtyp(ii).eq.2 .and. cpro_domtyp(jj).eq.3) then !anode|chute_anodique cpro_nbfrfl(ii)=jj cpro_nbfrfl(jj)=-ii cpro_intfc(ii)=ifac else if (cpro_domtyp(ii).eq.3 .and. cpro_domtyp(jj).eq.2) then !chute_anodique|anode cpro_nbfrfl(ii)=-jj cpro_nbfrfl(jj)=ii cpro_intfc(jj)=ifac end if end do do ifac = 1, nfac ii = ifacel(1,ifac) !IFACEL - Index-numbers (id est IEL) of the two (only) neighboring cells for each internal face. jj = ifacel(2,ifac) if (cpro_domtyp(ii).eq.3 .and. cpro_domtyp(jj).eq.3) then !chute_anodique|chute_anodique !if (abs(xyzcen(1,ii)-xyzcen(1,jj)).lt.1d-8 .and. abs(xyzcen(2,ii)-xyzcen(2,jj)).lt.1d-8) then if (cpro_nbfrfl(ii).lt.0 .and. cpro_nbfrfl(jj).eq.0) then ! ii is close to the anode ielhf=-cpro_nbfrfl(ii) !number of anode element cpro_nb_htfl(ielhf)=jj !number of fluid element else if (cpro_nbfrfl(jj).lt.0 .and. cpro_nbfrfl(ii).eq.0) then ! jj is close to the anode ielhf=-cpro_nbfrfl(jj) !number of anode element cpro_nb_htfl(ielhf)=ii !number of fluid element end if !end if end if end do end if !STEP N. Calculation of interface heat source. if (ntcabs-ntpabs.ge.2) then allocate(lstelt_cath(ncelet)) allocate(lstelt_anode(ncelet)) call field_get_label(ivarfl(isca(iscal)),label) !write(*,'(A7,A12,A11,I6)')'USTSSC ',label,' Source, t=',ntcabs ZZA3 = 0d0 ZZA4 = 0d0 ZZA5 = 0d0 ZZA6 = 0d0 ZZA7 = 0d0 ZZA13= 0d0 ZZA14= 0d0 ZZAT = 0d0 ZZ3 = 0d0 ZZ4 = 0d0 ZZ5 = 0d0 ZZ6 = 0d0 ZZ7 = 0d0 ZZ8 = 0d0 ZZ9 = 0d0 ZZ10 = 0d0 ZZT = 0d0 TWmax=0d0 !CATHODE call getcel('domaine_cathode',nlelt_cath,lstelt_cath) source_counter=0 source_sum=0d0 do ilelt = 1,nlelt_cath iel = lstelt_cath(ilelt) if (cpro_nbfrfl(iel).gt.0) then iels = iel !cathode cell (IFACEL 1, includes IF) ielp = cpro_nbfrfl(iel) !fluid cell (IFACEL 2, includes FJ) ifac = cpro_intfc(iel) !number of feca between cathode cell and fluid cell if (iels.gt.ncel) print*,'ERROR USTSSC iels>ncel ',iels if (ielp.gt.ncel) print*,'ERROR USTSSC ielp>ncel ',ielp if (iels.lt.1 .or. ielp.lt.1) print*,'ERROR USTSSC iel<1' if (ifac.gt.nfac) print*,'ERROR USTSSC ifac>nfac' if (ifac.lt.1) print*,'ERROR USTSSC ifac<1' ponds = pond(ifac) !(FJ·n)/(IJ·n) - share of fluid in assembly of two cells !pond(iface) - (FJ· n)/(IJ·n) for every internal face. F - center of face ifac. Ideal value is 0.5. surfloc = surfan(ifac)!SURFN = RA(ISRFAN-1+ifac) the norm of the surface of an internal face ifac distloc = dist(ifac) !The scalar product between the vectors IJ and n for every internal face. I and J are respectively the centers of the first and the second neighboring cell. The vector n is the UNIT VECTOR normal to the face and oriented from the first to the second cell. MODE = 4 !Cathode Tungsten thermal conductivity from temperature (file prop_W) CALL USTHHT (MODE,cpro_temp(iels),PROPs) ! PROPs = λ of cathode cell Just in case polynomial approx: PROPs=3.29E-12*TE**4 - 3.07E-08*TE**3 + 1.07E-04*TE**2 - 1.71E-01*TE+2.11E+02 PROPp = cpro_vscalt(ielp)*cpro_cp(ielp) !λ/Cp * Cp = λ of fluid cell TP = cpro_temp(ielp) !temperature of fluid cell TE = cpro_temp(iels) !temperature of cathode cell !Wall temperature calculation with respect to thermal conductivity !cor_tw_coef=ponds/(1d0-ponds)*PROPs/PROPp !TW = (TP/(1d0+cor_tw_coef) + TE/(1d0+1d0/cor_tw_coef)) !if (TW.ge.3.68d3) TW=3.68d3 !Clipping of Twall temperature. We don't have data for temperature above 3680K. TW=TE !fluid at the cathode tip is very cold therefore it's temperature doesn't have any physical meaning if(TW.gt.TWmax) TWmax=TW Jp=dabs((cpro_cur(1,ielp)*surfac(1,ifac)+cpro_cur(2,ielp)*surfac(2,ifac)+cpro_cur(3,ielp)*surfac(3,ifac))/surfan(ifac)) !current from fluid cells. Abs is necessary. Vectors j and n can be opposite. Js=dabs((cpro_cur(1,iels)*surfac(1,ifac)+cpro_cur(2,iels)*surfac(2,ifac)+cpro_cur(3,iels)*surfac(3,ifac))/surfan(ifac)) !current from cathode cells. Abs is necessary. Vectors j and n can be opposite. J=(Jp+Js)*0.5d0 !norm of current density component normal to the face ifac !LRichLaw=-1d0*eee*phimat/kb JEM = ARichLaw * TW * TW* exp(LRichLaw/TW)!The electron current density Jem is calculated from the Richardson-Dushman law !There is no Schottky effect. 4.2 eV is experimental data. We should take into account Schottky effect. JI = J - JEM !kin_enrg_coef=5d0*kb/(2d0*eee) QI = JI*(kin_enrg_coef*TW + Uk + phii ) !Heat source from ions current QEM = -1d0*JEM*(kin_enrg_coef*TW + phimat ) !Negative heat source from electron emission. Electrons carry out enthalpy. ElecSource = (QI + QEM)*surfan(ifac) !Source term from electron and ion currents. No Qed, Qray, Qevap. Qcond is acccounted in common part. !DiffSource=0d0 !if (TP.le.4d3 .and. ntcabs-ntpabs.ge.10) then ! DiffSource = -2d0*PROPs*PROPp*(TP-TE)/dist(ifac)/(PROPp*(1d0-ponds)+PROPs*ponds)*surfan(ifac) ! !write(*,'(A15,I5,A7,F10.4,A12,E15.5)')'Cathode cell N=',iels,' Twall=',TW,' DiffSource=',DiffSource !end if crvimp(iels) = 0d0 crvexp(iels) = crvexp(iels) + ElecSource! - DiffSource !at the beginning of time step crvexp(iels)=0 source_counter=source_counter+1 source_sum=source_sum+crvexp(iels) !CRVIMP(ielp) = 0d0 !crvexp(ielp) = crvexp(ielp) + 1d2*DiffSource !- ElecSource !if (ntcabs-ntpabs.eq.10000 .or. ntcabs-ntpabs.eq.40) write(*,'(A18,I5,A2,E15.5)')'cath fluid crvexp(',ielp,')=',crvexp(ielp) !WHERE IS OPPOSITE HEAT SOURCE FOR FLUID?! if (ntcabs-ntpabs.eq.99998) write(*,'(A16,I5,A2,E15.5)')' cathode crvexp(',iels,')=',crvexp(iels) LamFac = (PROPs*PROPp)/((PROPs*(1-Ponds))+(PROPp*Ponds)) !λcath*λfl/(λcath*γcath + λfl*γfl) ZZ2 = LamFac*(TP - TE)* surfan(ifac)/dist(ifac) ZZT = ZZT + ElecSource + ZZ2 !Total power transferred to the cathode = ',E15.7) ZZ3 = ZZ3 + ZZ2 !Convection at the cathode = ',E15.7) ZZ4 = ZZ4 + ElecSource !Heat characteristic of the polarity at the cathode = ',E15.7) Q = JI*kin_enrg_coef*TW* surfan(ifac) ZZ5 = ZZ5 + Q !Q = JI*TW*(5k/2e)*Surf = ',E15.7) Q = JI * Uk*surfan(ifac) ZZ6 = ZZ6 + Q !Q = JI*Uk = ',E15.7) Q = JI * phimat*surfan(ifac) ZZ7 = ZZ7 + Q !Q = JI*phimat = ',E15.7) QEM = -1d0*JEM*(kin_enrg_coef*TW + Uk + phimat)* surfan(ifac) ZZ8 = ZZ8 + QEM !QEM = -JEM(TW*(5k/2e) + Uk + phimat)*Surf = ',E15.7) Q = -1d0* JEM * phimat* surfan(ifac) ZZ9= ZZ9 + Q !QEM = -JEM*phimat = ',E15.7) Q = J* surfan(ifac) ZZ10 = ZZ10 + Q !QEM = -JEM*Uk = ',E15.7) end if end do !write(*,'(A12,I7,A,I6,A,E15.7)')label,ntcabs,' CathodeSource Counter=',source_counter,' Heat transferred to cathode=',source_sum if (irangp.ge.0) then call parsom (source_sum) call parcpt (source_counter) endif write(nfecra,'(A,I6,A,E15.7)')' Cathode Source Counter=',source_counter,' Heat transferred to cathode=',source_sum ! ANODE call getcel('anode',nlelt_anode,lstelt_anode) source_counter=0 source_sum=0d0 do ilelt = 1,nlelt_anode iel = lstelt_anode(ilelt) if (cpro_nbfrfl(iel).gt.0) then iels = iel !anode cell (IFACEL 1, includes IF) ielp = cpro_nbfrfl(iel) !fluid cell (IFACEL 2, includes FJ) ifac = cpro_intfc(iel) !number of feca between cathode cell and fluid cell if (iels.gt.ncel) print*,'ERROR USTSSC iels>ncel',iels if (ielp.gt.ncel) print*,'ERROR USTSSC ielp>ncel',ielp if (iels.lt.1 .or. ielp.lt.1) print*,'ERROR USTSSC iel<1' if (ifac.gt.nfac) print*,'ERROR USTSSC ifac>nfac' if (ifac.lt.1) print*,'ERROR USTSSC ifac<1' surfloc = surfan(ifac)!SURFN = RA(ISRFAN-1+ifac) the norm of the surface of an internal face ifac distloc = dist(ifac) !The scalar product between the vectors IJ and n for every internal face. I and J are respectively the centers of the first and the second neighboring cell. The vector n is the UNIT VECTOR normal to the face and oriented from the first to the second cell. ponds = pond(ifac) !(FJ·n)/(IJ·n) - share of fluid in assembly of two cells !pond(iface) - (FJ· n)/(IJ·n) for every internal face. F - center of face ifac. Ideal value is 0.5. !if (ntcabs-ntpabs.eq.2) write(*,'(A,I6,A,I6)')'anode cell N=',iels,' fluid cell N=',ielp PROPs = 4.7d2- ( cpro_temp(iels) / 1d1 ) !anode λ - thermal conductivity from temperature PROPp = cpro_vscalt(ielp)*cpro_cp(ielp) !λ/Cp * Cp = λ of fluid cell TP = cpro_temp(ielp) !temperature of fluid cell TE = cpro_temp(iels) !temperature of anode cell !Wall temperature calculation with respect to thermal conductivity cor_tw_coef=ponds/(1d0-ponds)*PROPs/PROPp TW = (TP/(1d0+cor_tw_coef) + TE/(1d0+1d0/cor_tw_coef)) Jp=dabs((cpro_cur(1,ielp)*surfac(1,ifac)+cpro_cur(2,ielp)*surfac(2,ifac)+cpro_cur(3,ielp)*surfac(3,ifac))/surfan(ifac)) !current from fluid cells. Abs is necessary. Vectors j and n can be opposite. Js=dabs((cpro_cur(1,iels)*surfac(1,ifac)+cpro_cur(2,iels)*surfac(2,ifac)+cpro_cur(3,iels)*surfac(3,ifac))/surfan(ifac)) !current from anode cells. Abs is necessary. Vectors j and n can be opposite. J=(Jp+Js)*0.5d0 !norm of current density component normal to the face ifac !rayo = sqrt(cdgfac(1,ifac)*cdgfac(1,ifac)+cdgfac(2,ifac)*cdgfac(2,ifac)) !if (ntcabs-ntpabs.eq.20) print*,J,surfan(ifac),rayo,xyzcen(1,iels),xyzcen(2,iels) !**************************************************************** !The heat flux provided by the recombination of the ions at the surface is negligible !The heat flux related to the electronic current is generally the majority term ElecSource = J*(kin_enrg_coef*TW + Ua + phimat) !Why phimat is the same for anode ? crvimp(iels) = 0d0 crvexp(iels) = crvexp(iels) + ElecSource * surfan(ifac) !at the beginning of time step crvexp(iels)=0 source_sum=source_sum+crvexp(iels) source_counter=source_counter+1 if (ntcabs-ntpabs.eq.99999) write(*,'(A16,I5,A2,E15.5)')' anode crvexp(',iels,')=',crvexp(iels) !WHERE IS OPPOSITE HEAT SOURCE FOR FLUID?! !write(*,'(I3,A16,I5,A3,E12.3)')ntcabs,' anode crvexp(',iels,')= ',crvexp(iels) !IDK what is this !SUMM =(PROPs*PROPp)/((PROPs*(1-Ponds))+(PROPp*Ponds)) !λcath*λfl/(λcath*γcath + λfl*γfl) IDK what is this !Calculation of convective flux taking into account temperature of the second layer of fluid cells. The first layer is too cold and uselss. ielhf=cpro_nb_htfl(iels) !number of fluid element one layer distant from the anode !TP=0.5d0*(cpro_temp(ielhf)+cpro_temp(ielp)) !PROPp = 0.5d0*(cpro_vscalt(ielhf)*cpro_cp(ielhf)+cpro_vscalt(ielp)*cpro_cp(ielp)) !LamFac = (PROPs*PROPp)/((PROPs*(1-Ponds))+(PROPp*Ponds)) !λcath*λfl/(λcath*γcath + λfl*γfl) ZZA1 =ElecSource * surfan(ifac) distloc = 0.004d0 !Thickness of the cathode !ZZA2 = -2d0*LamFac*(TE - TP) * surfan(ifac)/dist(ifac) ZZA2 = PROPs*(TE - 300d0) * surfan(ifac)/distloc !Thermal flux inside the anode from top to bottom ZZAT = ZZAT + ZZA2 + ZZA1 !Total power transferred to the anode = ',E15.7) if (TE.GT.450d0) then ZZA3 = ZZA3 + ZZA2 !Convection at the anode = ',E15.7) end if ZZA4 = ZZA4 + ZZA1 !Heat characteristic of the polarity at the anode = ',E15.7) Q = J*kin_enrg_coef*TW* surfan(ifac) ZZA5 = ZZA5 + Q !Q = JI*TW*(5k/2e)*Surf = ',E15.7) Q = J * Ua* surfan(ifac) ZZA6 = ZZA6 + Q !Q = JI*Ua = ',E15.7) Q = J * phimat* surfan(ifac) ZZA7 = ZZA7 + Q !Q = JI*phimat = ',E15.7) if (TE.GT.400d0) then ZZA13 = ZZA13 + ZZA2 !Partial convection at the anode 2 = ',E15.7) end if ZZA14 = ZZA14 + ZZA2 !Total convection at the anode = ',E15.7) end if end do !display output at every time step lead to memory leak and for big time step number the case crashes. !write(*,'(A12,I7,A,I6,A,E15.7)')label,ntcabs,' AnodeSource Counter=',source_counter,' Heat transferred to anode =',source_sum if (irangp.ge.0) then call parsom (source_sum) call parcpt (source_counter) endif write(nfecra,'(A,I6,A,E15.7)')' Anode Source Counter=',source_counter,' Heat transferred to anode =',source_sum deallocate(lstelt_cath) deallocate(lstelt_anode) !SUMMATION IN CASE OF PARALLELISM if (irangp.ge.0) then call parsom (ZZT) call parsom (ZZ3) call parsom (ZZ4) call parsom (ZZ5) call parsom (ZZ6) call parsom (ZZ7) call parsom (ZZ8) call parsom (ZZ9) call parsom (ZZ10) call parsom (ZZAT) call parsom (ZZA3) call parsom (ZZA4) call parsom (ZZA5) call parsom (ZZA6) call parsom (ZZA7) call parsom (ZZA13) call parsom (ZZA14) call parmax (TWmax) !compute of global maximum end if !PRINT TO LISTING write(nfecra,1008) write(nfecra,1022)ZZT !Total power transferred to the cathode = ',E15.7) write(nfecra,1023)ZZ3 !Convection at the cathode = ',E15.7) write(nfecra,1024)ZZ4 !Heat characteristic of the polarity at the cathode = ',E15.7) write(nfecra,1025)ZZ5 !Q = JI*TW*(5k/2e)*Surf = ',E15.7) write(nfecra,1026)ZZ6 !Q = JI*Uk = ',E15.7) write(nfecra,1027)ZZ7 !Q = JI*phimat = ',E15.7) write(nfecra,1028)ZZ8 !QEM = -JEM(TW*(5k/2e) + Uk + phimat)*Surf = ',E15.7) write(nfecra,1029)ZZ9 !QEM = -JEM*phimat = ',E15.7) write(nfecra,1030)ZZ10 !Total current at the cathode',E15.7 write(nfecra,1004)TWmax !Cathode max Twall = ',E15.7 write(nfecra,1008) write(nfecra,1008) write(nfecra,1012)ZZAT !Total power transferred to the anode = ',E15.7) write(nfecra,1013)ZZA3 !Convection at the anode = ',E15.7) write(nfecra,1014)ZZA4 !Heat characteristic of the polarity at the anode = ',E15.7) write(nfecra,1015)ZZA5 !Q = JI*TW*(5k/2e)*Surf = ',E15.7) write(nfecra,1016)ZZA6 !Q = JI*Ua = ',E15.7) write(nfecra,1017)ZZA7 !Q = JI*phimat = ',E15.7) write(nfecra,1031)ZZA13 !Partial convection at the anode 2 = ',E15.7) write(nfecra,1032)ZZA14 !Total convection at the anode = ',E15.7) write(nfecra,1008) end if endif !-------- ! Formats !-------- 1101 format('ustssc t=',I7,' iscal=',I4,' ncepdp=',I7,' ncesmp=',I7,' nvar=',I7) 1102 format('ustssc t=',I3,' Name of ivarfl(isca(',I3,')) = ',A22,' label=',A22) 1103 format('CellN=',I5,' Tcath=',F11.3,' Tfluid=',F11.3,' cor_tw_coef=',E13.5,' Twall=',F11.5) 1104 format('LambdaCathode=',F15.7,' CpCathode=',F15.7,' LambdaFluid=',F15.7,' CpFluid=',F15.7) 1000 format(' User source terms for variable ',A8,/) 1001 format(' Calculated Current = ',E15.7) 1012 format(' Total power transferred to the anode = ',E15.7) 1013 format(' Convection at the anode = ',E15.7) 1014 format(' Heat characteristic of the polarity at the anode = ',E15.7) 1015 format(' Q = JI*TW*(5k/2e)*Surf = ',E15.7) 1016 format(' Q = JI*Ua = ',E15.7) 1017 format(' Q = JI*phimat = ',E15.7) 1022 format(' Total power transferred to the cathode = ',E15.7) 1023 format(' Convection at the cathode = ',E15.7) 1024 format(' Heat characteristic of the polarity at the cathode = ',E15.7) 1025 format(' Q = JI*TW*(5k/2e)*Surf = ',E15.7) 1026 format(' Q = JI*Uk = ',E15.7) 1027 format(' Q = JI*phimat = ',E15.7) 1028 format(' QEM = -JEM(TW*(5k/2e) + Uk + phimat)*Surf = ',E15.7) 1029 format(' QEM = -JEM*phimat = ',E15.7) 1030 format(' Total current at the cathode',E15.7) 1003 format(' TP = ',E15.7) 1004 format(' Cathode max Twall = ',E15.7) 1005 format(' ZZ1 = ',E15.7) 1006 format(' ZZ2 = ',E15.7) ! 1007 FORMAT('ustssc I AM HERE') 1008 format('******* USTSSC *******') 1031 format('Partial convection at the anode 2 = ',E15.7) 1032 format('Total convection at the anode = ',E15.7) !---- ! End !---- return end subroutine ustssc