!------------------------------------------------------------------------------- ! Code_Saturne version 4.2.1 ! -------------------------- ! This file is part of Code_Saturne, a general-purpose CFD tool. ! ! Copyright (C) 1998-2016 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. !------------------------------------------------------------------------------- module elincl !============================================================================= use paramx use ppthch !============================================================================= !--> DEFINITION DES PARAMETERS ! ========================= ! PERMVI : Mu zero, permeabilite magnetique du vide H/m ! EPSZER : Epsilon zero, permittivite du vide F/m double precision permvi , epszer parameter ( permvi = 1.2566d-6, epszer = 8.854d-12 ) !--> DONNEES EN COMMON POUR LE CHAUFFAGE EFFET JOULE ! =============================================== ! TH, NPOT et NPO sont deja dans ppthch.h ! ----- Fournies par l'utilisateur ! IENTM1 --> indicateur entree matiere premiere ! IELPH1 --> indicateur electrode phase 1 ! IELPH2 --> indicateur electrode phase 2 ! IELPH3 --> indicateur electrode phase 3 ! IELNEU --> indicateur electrode neutre ! ENH --> tabulation enthalpie(temperature) ! USRHO --> - - - - - inverse masse volumique - - - ! SIG --> - - - - - conductivite - - - - ! KAB --> - - - - - coeff absorption - - ! VIS --> - - - - - viscosite - - - - - - ! LCP --> - - - - - Lambda/Cp integer, save :: ientm1(ntypmx),ielph1(ntypmx),ielph2(ntypmx) integer, save :: ielph3(ntypmx),ielneu(ntypmx) ! ENHEL --> tabulation enthalpie (temperature) ! RHOEL --> - - - - - masse volumique - - - ! CPEL --> - - - - - CP - - - ! SIGEL --> - - - - - conductivite elec - - - - ! XLABEL --> - - - - - conductivite thermique - - ! XKABEL --> - - - - - coeff absorption (pour Srad)- - ! VISEL --> - - - - - viscosite dynamique - - - - - - double precision, save :: rhoel (ngazgm,npot), cpel (ngazgm,npot) double precision, save :: sigel (ngazgm,npot), visel (ngazgm,npot) double precision, save :: xlabel(ngazgm,npot), xkabel(ngazgm,npot) ! CL sur les electrodes integer nelemx,nbtrmx parameter (nelemx = 1000 , nbtrmx = 100) integer, save :: nbelec , nbtrf , ntfref integer, save :: ielecc(nelemx),ielect(nelemx),ielecb(nelemx) integer, save :: ibrpr(nbtrmx),ibrsec(nbtrmx) double precision, save :: tenspr(nbtrmx),rnbs(nbtrmx) double precision, save :: zr(nbtrmx) ,zi(nbtrmx) double precision, save :: uroff(nbtrmx) ,uioff(nbtrmx) !--> PARAMETRES POUR LA VERSION ARC ELECTRIQUE ! ======================================== ! IXKABE : valeur lue dans le fichier dp_elec ! = 0 la derniere colonne du fichier est lue mais pas utilisee ! = 1 la derniere colonne du fivhier represente le coefficient ! d'absorption ! = 2 la derniere colonne du fivhier represente le TS radiatif integer, save :: ixkabe ! Grandeurs necessaires au claquage ! NTDCLA : iterration de debut du claquage ! ICLAQ : indicateur pour savoir si on fait actuellement un claquage ! = 0 Pas de claquage ! = 1 Claquage ! XCLAQ ,YCLAQ ZCLAQ : Position de point de claquage integer, save :: ntdcla , iclaq double precision, save :: xclaq , yclaq , zclaq !--> DONNEES SUR LA CORRECTION DES VARIABLES ELECTRIQUES ! EN FONCTION D'UNE INTENSITE DE COURANT DONNEES ! ======================================== ! IELCOR : = 0 pas de correction ! = 1 correction ! COUIMP : intensite de courant impose par l'utilisateur ! pour Arc Electrique ! PUISIM : puissance imposee pour Joule ! DPOT : Delta du potentiel electrique entre l'Anode et la cathode ! (arc et Joule) ! COEJOU : coefficient de correction pour version Joule integer, save :: ielcor double precision, save :: couimp , dpot , puisim , coejou, elcou !--> DONNEES POUR LES ESPECES AYANT UN IMPACT ! SUR LE PROBLEME ELECTRIQUE ! ======================================== ! QESPEL : Charge massique des especes C/kg ! SUSCEP : Susceptibilite (relation champ - mobilite) m2/s/V double precision, save :: qespel(ngazgm), suscep(ngazgm) !============================================================================= end module elincl subroutine uselph & !================ ( nvar , nscal , & mbrom , izfppp , & dt ) !=============================================================================== ! FONCTION : ! -------- ! REMPLISSAGE DES VARIABLES PHYSIQUES POUR LE MODULE ELECTRIQUE ! ----> Effet Joule ! ----> Arc Electrique ! ----> Conduction Ionique ! 1) Masse Volumique ! 2) Viscosite moleculaire ! 3) Chaleur massique Cp ! 4) Lambda/Cp moleculaire ! 4) Diffusivite moleculaire ! ATTENTION : ! ========= ! Il est INTERDIT de modifier la viscosite turbulente VISCT ici ! ======== ! (une routine specifique est dediee a cela : usvist) ! Pour le module electrique, toutes les proprietes physiques sont ! supposees variables et contenues dans le tableau PROPCE ! (meme si elles sont physiquement constantes) ! Remarques : ! --------- ! Cette routine est appelee au debut de chaque pas de temps ! Ainsi, AU PREMIER PAS DE TEMPS (calcul non suite), les seules ! grandeurs initialisees avant appel sont celles donnees ! - dans usipsu : ! . la masse volumique (initialisee a RO0) ! . la viscosite (initialisee a VISCL0) ! - dans usiniv/useliv : ! . les variables de calcul (initialisees a 0 par defaut ! ou a l ! ou a la valeur donnee dans usiniv) ! On peut donner ici les lois de variation aux cellules ! - de la masse volumique ROM kg/m3 ! (et eventuellememt aux faces de bord ROMB kg/m3) ! - de la viscosite moleculaire VISCL kg/(m s) ! - de la chaleur specifique associee CP J/(kg degres) ! - des "diffusivites" associees aux scalaires VISCLS kg/(m s) ! On dispose des types de faces de bord au pas de temps ! precedent (sauf au premier pas de temps, ou les tableaux ! ITYPFB et ITRIFB n'ont pas ete renseignes) ! Il est conseille de ne garder dans ce sous programme que ! le strict necessaire. ! Cells identification ! ==================== ! Cells may be identified using the 'getcel' subroutine. ! The syntax of this subroutine is described in the ! 'cs_user_boundary_conditions' subroutine, ! but a more thorough description can be found in the user guide. !------------------------------------------------------------------------------- ! Arguments !__________________.____._____.________________________________________________. ! name !type!mode ! role ! !__________________!____!_____!________________________________________________! ! nvar ! i ! <-- ! total number of variables ! ! nscal ! i ! <-- ! total number of scalars ! ! mbrom ! te ! <-- ! indicateur de remplissage de romb ! ! izfppp ! te ! <-- ! numero de zone de la face de bord ! ! (nfabor) ! ! ! pour le module phys. part. ! ! dt(ncelet) ! ra ! <-- ! time step (per cell) ! !__________________!____!_____!________________________________________________! ! Type: i (integer), r (real), s (string), a (array), l (logical), ! and composite types (ex: ra real array) ! mode: <-- input, --> output, <-> modifies data, --- work array !=============================================================================== !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use optcal use cstphy use entsor use ppppar use ppthch use ppincl use elincl use field use mesh use cs_c_bindings !=============================================================================== implicit none ! Arguments integer nvar , nscal integer mbrom integer izfppp(nfabor) double precision dt(ncelet) ! Local variables integer iel, ifcvsl,ifcsig, ilelt, nlelt, mode, potr_id integer, allocatable, dimension(:) :: lstelt double precision TungsDens001, CopperDens001, Wk1, Wk2 double precision tp double precision xkr , xbr double precision rom0 , temp0 , dilar , aa , bb , cc double precision srrom1, rhonp1 double precision, allocatable, dimension(:) :: temp double precision, dimension(:), pointer :: cvar_hm, cpro_crom double precision, dimension(:), pointer :: cpro_rom double precision, dimension(:), pointer :: cpro_viscl, cpro_vscalt, cpro_sig double precision, dimension(:), pointer :: cpro_vpoti, cpro_vpotr double precision, dimension(:), pointer :: cpro_cp, cpro_temp double precision, dimension(:), pointer :: cvar_scalt double precision, dimension(:), pointer :: chm integer ipass data ipass /0/ save ipass !=============================================================================== !=============================================================================== ! 0 - INITIALISATIONS A CONSERVER !=============================================================================== ! --- Initialisation memoire ipass = ipass + 1 ! Message au premier passage pour indiquer que l'utilisateur a ! rapatrie le sous-programme. if (ipass.eq.1) then write(nfecra,1000) endif !=============================================================================== ! 2 - ARC ELECTRIQUE !=============================================================================== if (ntcabs.eq.0 .or. ntcabs.eq.1) print*,'uselph start' if ( ippmod(ielarc).ge.1 ) then !Tungsten density * 0.01 (file H-T-Cp_W): TungsDens001 = 19.3D3*1D-2 Wk1=4.28D-6 Wk2=5.8D-10 !dV = 8%, V=V0*(1+(4.28T+0.00058T^2)*10^-6)^3, Ro=Ro0*V0/V !By Peter Hidnert and W. T. Sweeney THERMAL EXPANSION OF TUNGSTEN !Works below 3600K. "Thermal Expansion of Tungsten in the Range 1500-3600 K by a transient interferometric technique", A. P. Miiller 2 and A. Cezairliyan !Copper density * 0.01: CopperDens001 = 8.92D1 allocate(lstelt(ncelet)) ! temporary array for cells selection allocate(temp(ncelet)) call field_get_val_s(ivarfl(isca(ihm)), cvar_hm) call field_get_val_s(icrom, cpro_crom) if (icp.ge.0) call field_get_val_s(icp, cpro_cp) call field_get_key_int(ivarfl(isca(iscalt)), kivisl, ifcvsl) if (ifcvsl.ge.0) then call field_get_val_s(ifcvsl, cpro_vscalt) else cpro_vscalt => NULL() endif !call field_get_key_int (ivarfl(isca(elec_pot_r)), kivisl, ifcsig) call field_get_id('elec_pot_r',potr_id) call field_get_key_int(potr_id, kivisl, ifcsig) print*,'potr_id=',potr_id print*,'ifcsig= ',ifcsig if (ifcsig.ge.0) then call field_get_val_s(ifcsig, cpro_sig) else cpro_sig => NULL() endif write(*,1018)'uselph cathd t=',ntcabs call getcel('domaine_cathode',nlelt,lstelt) do ilelt = 1, nlelt iel = lstelt(ilelt) !Tungsten temperature from cathode enthalpy (file H-T-Cp_W): mode = 1 call usthht (mode,cvar_hm(iel),temp(iel)) !Tungsten density * 0.01 (file H-T-Cp_W): !Ro=Ro0*V0/V !cpro_crom(iel) = TungsDens001/(1+Wk1*temp(iel)+Wk2*temp(iel)**2)**3 cpro_crom(iel) = TungsDens001 !Tungsten specific heat in J/(kg degres) from cathode enthalpy (file H-T-Cp_W): mode = 0 if (icp.ge.0) call usthht (mode,cvar_hm(iel),cpro_cp(iel)) !cpro_cp(iel)=300.0 !Firstly Lambda in kg/(m s) is calculated - Tungsten thermal conductivity from temperature (file prop_W) mode = 4 if(ifcvsl.GT.0) call usthht (mode,temp(iel),cpro_vscalt(iel)) !We divide by CP to have Lambda/CP. We assume that Cp is prefixed. if(icp.le.0) then !If Cp is uniform, we use cp0, cp0 is the reference specific heat !cpro_vscalt(iel) = cpro_vscalt(iel)/cp0 else !If Cp is non uniform, we use the CP calculation above. !cpro_vscalt(iel) = cpro_vscalt(iel)/cpro_cp(iel)!IPCCP IPCVSL endif !Tungsten electrical conductivity in S/m from temperature (file prop_W) mode = 3 if (ifcsig.ge.0) call usthht (mode,temp(iel),cpro_sig(iel))!IPCSIG if (ilelt.eq.3) then write(nfecra,1014)temp(iel) write(nfecra,1015)cpro_sig(iel) write(nfecra,1016)cpro_vscalt(iel) write(nfecra,1017)cpro_cp(iel) endif enddo call getcel('anode',nlelt,lstelt) write(*,1018)'uselph anode t=',ntcabs do ilelt = 1, nlelt iel = lstelt(ilelt) mode = 2 !Copper temperature from anode enthalpy (file H-T_Cu2): call usthht (mode,cvar_hm(iel),temp(iel)) !Copper density * 0.01: cpro_crom(iel) = CopperDens001 !Specific heat J/(kg degres) if (temp(iel).LT.1.37D3) then !if (icp.ge.0) cpro_cp(iel) = 3.7D2+(temp(iel)/ 1.0D1 ) else cpro_cp(iel) = 5.07D2 endif !Lambda/Cp in kg/(m s) if(ifcvsl.GT.0) then !Firstly Lambda is calculated cpro_vscalt(iel) = 4.7D2-( temp(iel) / 1.0D1 ) if(cpro_vscalt(iel).LT.0) then print*,'error anode' write(nfecra,1010)cpro_vscalt(iel) write(nfecra,1011)XYZCEN(1,iel),XYZCEN(2,iel),XYZCEN(3,iel) endif !We divide by Cp to have Lambda/CP. We assume that Cp is prefixed. if(icp.le.0) then !If Cp is uniform, we use cp0 !cpro_vscalt(iel) = cpro_vscalt(iel)/cp0 else !If Cp is non uniform, we use the CP calculation above. !cpro_vscalt(iel) = cpro_vscalt(iel)/cpro_cp(iel)!IPCCP IPCVSL endif !Electrical conductivity in S/m from anode temperature if (ifcsig.ge.0) cpro_sig(iel) = 4.5D7*exp(-temp(iel)/1.0D3) endif enddo print*,'uselph finit' call getcel('chute_cathodique',nlelt,lstelt) do ilelt = 1, nlelt iel = lstelt(ilelt) cpro_sig(iel) = 8.0D3 enddo call getcel('chute_anodique',nlelt,lstelt) do ilelt = 1, nlelt iel = lstelt(ilelt) cpro_sig(iel) = 8.0D3 enddo deallocate(lstelt) deallocate(temp) endif !=============================================================================== ! 3 - CONDUCTION IONIQUE !=============================================================================== ! CETTE OPTION N'EST PAS ACTIVABLE !-------- ! Formats !-------- 1000 format(/, & ' Electrical module: user intervention for ',/, & ' the calculation of physical properties.',/) 1010 format(/,'erreur = ',E14.5) 1011 format(/,'POSX = ',E14.5,/, 'POSY = ',E14.5,/, 'POSY = ',E14.5 ) 1012 format(/,'temp = ',E14.5) 1013 format(/,'enthalpy = ',E14.5) 1014 format(/,'temp= ',E14.5) 1015 format(/,'cpro_sig= ',E14.5) 1016 format(/,'cpro_vscalt=',E14.5) 1017 format(/,'cpro_cp= ',E14.5) 1018 format(A,I6) !---- ! End !---- return end subroutine uselph