!------------------------------------------------------------------------------- ! 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. !------------------------------------------------------------------------------- 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 parall use paramx use numvar use optcal use cstphy use entsor use ppppar use ppthch use ppincl 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 ifcp,ifcvsl,ifcsig,ifctemp!,ifvisc !pointers integer iel,ilelt, nlelt, mode !numbers, indexes integer, allocatable, dimension(:) :: lstelt !cells array double precision TungsDens001, CopperDens001, Wk1, Wk2 !material properties double precision TungsDens, CopperDens double precision, allocatable, dimension(:) :: temp double precision, dimension(:), pointer :: cvar_hm, cpro_crom double precision, dimension(:), pointer :: cpro_rom double precision, dimension(:), pointer :: cpro_vscalt, cpro_sig !,cpro_viscl double precision, dimension(:), pointer :: cpro_cp, cpro_temp double precision rayo, center_sigm, quadr_factr,artif_sig !double precision rom0 , temp0 , dilar , aa , bb , cc !double precision srrom1, rhonp1 !double precision coorx,coory,coorz,elcurx,elcury,elcurz character(len=80) :: name 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 !write(*,'(A,I6)')'uselph start, t=',ntcabs !=============================================================================== ! 2 - ARC ELECTRIQUE !=============================================================================== if ( ippmod(ielarc).ge.1 ) then !*************************************************************************************** ! POINTERS * !*************************************************************************************** if (ntcabs-ntpabs.eq.1) print*,'POINTERS TEST:' !call field_get_key_id("variable_id", keyvar) !call field_get_id('enthalpy', f_id) !call field_get_key_int(f_id, keyvar, ipotr) !call field_get_val_s(ivarfl(isca(ihm)), cvar_hm) !call field_get_id('vec_potential', f_id) !call field_get_key_int(f_id, keyvar, ipotva) call field_get_val_s(ivarfl(isca(ihm)), cvar_hm) !ENTHALPY. Obtaining the field by the pointer. Actually id of enthalpy field is 2. call field_get_name(ivarfl(isca(ihm)),name) if (ntcabs-ntpabs.eq.1) write(*,1021)ivarfl(isca(ihm)),name !test of enthalpy pointer ihm and iscalt. They must be equal. call field_get_name(ivarfl(isca(iscalt)),name) if (ntcabs-ntpabs.eq.1) write(*,1022)ivarfl(isca(iscalt)),name !I found ihm and iscalt to be the same 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 (ntcabs-ntpabs.eq.1 .and. ifctemp.eq.0) print*,'No temperature id. cs_user_physical_properties.f90' !error message just in case. call field_get_val_s(icrom, cpro_crom) !DENSITY call field_get_name(icrom,name) if (ntcabs-ntpabs.eq.1) write(*,1023)icrom,name 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 call field_get_name(ifcp,name) if (ntcabs-ntpabs.eq.1) write(*,1024)ifcp,name call field_get_key_int(ivarfl(isca(ihm)), kivisl, ifcvsl) !LAMBDA. Dynamic difusivity of enthalpy, it is the ratio between thermal conductivity and specific heat in the enthalpy equation. if (ifcvsl.ge.0) then !test just in case call field_get_val_s(ifcvsl, cpro_vscalt) call field_get_name(ifcvsl,name) if (ntcabs-ntpabs.eq.1) write(*,1025)ifcvsl,name !test of Lambda pointer else cpro_vscalt => NULL() if (ntcabs-ntpabs.eq.1) print*,'ERROR: No Lambda id. cs_user_physical_properties.f90' endif call field_get_id_try('elec_sigma', ifcsig) !SIGMA. Diffusivity of electric potential. if (ifcsig.ge.0) then call field_get_val_s(ifcsig, cpro_sig) call field_get_name(ifcsig,name) if (ntcabs-ntpabs.eq.1) write(*,1026)ifcsig,name else cpro_sig => NULL() if (ntcabs-ntpabs.eq.1) print*,'ERROR: No Sigma id. cs_user_physical_properties.f90' endif !call field_get_val_s(iviscl, cpro_viscl) !call field_get_id_try('molecular_viscosity', ifvisc) !MOLECULAR VISCOSITY !call field_get_val_s(ifvisc, cpro_viscl) !call field_get_name(ifvisc,name) !if (ntcabs.eq.1) write(*,1027)ifvisc,name !*************************************************************************************** ! PHYSICAL PROPERTIES * !*************************************************************************************** !Tungsten density * 0.01 (file H-T-Cp_W): TungsDens001 = 19.3D3*1D-2 TungsDens = 19.3D3 Wk1=4.28D-6 !coefficients for linear approximation of temperature dependency of tungsten density 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 till 3600 Thermal Expansion of Tungsten in the Range 1500-3600 K by a Transient Interferometric Technique 1 A. P. Miiller 2 and A. Cezairliyan 2 !Copper density * 0.01: CopperDens001 = 8.92D1 CopperDens = 8.94D3 allocate(lstelt(ncelet)) ! temporary array for cells selection allocate(temp(ncelet)) call getcel('domaine_cathode',nlelt,lstelt) do ilelt = 1, nlelt !if (cvar_hm(iel).le.14000.0) cvar_hm(iel)=14000.0 iel = lstelt(ilelt) !cpro_viscl(iel)=1d-1 !Tungsten temperature from cathode enthalpy (file H-T-Cp_W): mode = 1 call usthht (mode,cvar_hm(iel),temp(iel)) if (ifctemp.ge.0) cpro_temp(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!TungsDens! !Tungsten specific heat in J/(kg degres) from cathode enthalpy (file H-T-Cp_W): mode = 0 if (ifcp.ge.0) call usthht (mode,cvar_hm(iel),cpro_cp(iel)) !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)) !W/(m*K) !We divide by CP to have Lambda/CP. We assume that Cp is prefixed. if(ifcp.lt.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,'(A28)')'Example of metal properties:' write(nfecra,1019)cvar_hm(iel),temp(iel),cpro_crom(iel),cpro_cp(iel),cpro_vscalt(iel),cpro_sig(iel) endif enddo call getcel('anode',nlelt,lstelt) do ilelt = 1, nlelt iel = lstelt(ilelt) !cpro_viscl(iel)=1d-1 mode = 2 !Copper temperature from anode enthalpy (file H-T_Cu2): call usthht (mode,cvar_hm(iel),temp(iel)) if (ifctemp.ge.0) cpro_temp(iel) =temp(iel) !Copper density * 0.01: cpro_crom(iel) = CopperDens001!CopperDens! !Specific heat J/(kg degres) if (temp(iel).LT.1.37D3) then !if (ifcp.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 ) !W/(m*K) 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(ifcp.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 if (ilelt.eq.3) then write(nfecra,1020)cvar_hm(iel),temp(iel),cpro_crom(iel),cpro_cp(iel),cpro_vscalt(iel),cpro_sig(iel) endif enddo call getcel('chute_cathodique',nlelt,lstelt) do ilelt = 1, nlelt iel = lstelt(ilelt) rayo = sqrt(xyzcen(1,iel)*xyzcen(1,iel)+xyzcen(2,iel)*xyzcen(2,iel)) !Sigma=center_sigm+quadr_factr*R2 - parabolic distribution center_sigm = 5.073d3 !Electrical conductivity value at the domain central axis quadr_factr = -1.722d10 artif_sig = quadr_factr*rayo*rayo+center_sigm !Old version: cpro_sig(iel) = 8.0D3 if(artif_sig.lt.0.13214d-3) artif_sig = 0.13214d-3 !Clipping. Value = Electrical conductivity of Argon at 300K if(xyzcen(3,iel).le.4d-5) then cpro_sig(iel)=artif_sig*0.2d0 else if (xyzcen(3,iel).gt.4d-5 .and. xyzcen(3,iel).le.1d-4) then if(cpro_sig(iel).lt.artif_sig*0.8d0) cpro_sig(iel)=artif_sig*0.8d0 else if (xyzcen(3,iel).gt.1d-4) then if(cpro_sig(iel).lt.artif_sig) cpro_sig(iel)=artif_sig end if !cpro_sig(iel)=8d3 enddo call getcel('chute_anodique',nlelt,lstelt) do ilelt = 1, nlelt iel = lstelt(ilelt) rayo = sqrt(xyzcen(1,iel)*xyzcen(1,iel)+xyzcen(2,iel)*xyzcen(2,iel)) !Sigma=center_sigm+quadr_factr*R2 - parabolic distribution center_sigm = 4.0d3 !Electrical conductivity value at the domain central axis quadr_factr = -9.9d8 artif_sig = quadr_factr*rayo*rayo+center_sigm !Old version: cpro_sig(iel) = 8.0D3 if(artif_sig.lt.0.13214d-3) artif_sig = 0.13214d-3 !Clipping. Value = Electrical conductivity of Argon at 300K if (cpro_sig(iel).lt.artif_sig) cpro_sig(iel) = artif_sig 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 ) 1019 format('H= ',E14.5,' W T= ',E14.5,' W Ro*0.01= ',E14.5,' W Cp= ',E14.5,' W Lambda= ',E14.5,' W Sigma= ',E14.5 ) 1020 format('H= ',E14.5,' CuT= ',E14.5,' CuRo*0.01= ',E14.5,' CuCp= ',E14.5,' CuLambda= ',E14.5,' CuSigma= ',E14.5 ) 1021 format(' ihm Id=',I4,' F_name=',A22,' Must be Id=2 name=enthalpy') 1022 format(' iscalt Id=',I4,' F_name=',A22,' Must be Id=2 name=enthalpy') 1023 format(' Density Id=',I4,' F_name=',A22,' Must be Id=7 name=density') 1024 format(' Cp Id=',I4,' F_name=',A22,' Must be Id=23 name=specific_heat') 1025 format(' Lambda Id=',I4,' F_name=',A22,' Must be Id=25 name=thermal_diffusivity') 1026 format(' Sigma Id=',I4,' F_name=',A22,' Must be Id=26 name=elec_pot_r_diffusivity') 1027 format('ViscosityId=',I4,' F_name=',A22,' Must be Id=9 name=molecular_viscosity') ! 1019 format(/,' W temp= ',E14.5,/,' W Ro*0.01= ',E14.5,/,' W Cp= ',E14.5,/,' W lambda= ',E14.5,/,' W sigma= ',E14.5 ) ! 1020 format(/,'Cu temp= ',E14.5,/,'Cu Ro*0.01= ',E14.5,/,'Cu Cp= ',E14.5,/,'Cu lambda= ',E14.5,/,'Cu sigma= ',E14.5 ) !---- ! End !---- return end subroutine uselph