!------------------------------------------------------------------------------- ! 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. !------------------------------------------------------------------------------- !=============================================================================== ! Function: ! --------- !> \file cs_user_boundary_conditions.f90 !> !> \brief User subroutine which fills boundary conditions arrays !> (\c icodcl, \c rcodcl) for unknown variables. !> !> See \subpage cs_user_boundary_conditions_examples for examples. !> !> \section cs_user_boundary_conditions_intro Introduction !> !> Here one defines boundary conditions on a per-face basis. !> !> Boundary faces may be selected using the \ref getfbr subroutine. !> !> \code getfbr(string, nelts, lstelt) \endcode !> - string is a user-supplied character string containing selection criteria; !> - nelts is set by the subroutine. It is an integer value corresponding to !> the number of boundary faces verifying the selection criteria; !> - lstelt is set by the subroutine. It is an integer array of size nelts !> containing the list of boundary faces verifying the selection criteria. !> !> string may contain: !> - references to colors (ex.: 1, 8, 26, ...) !> - references to groups (ex.: inlet, group1, ...) !> - geometric criteria (ex. x < 0.1, y >= 0.25, ...) !> !> These criteria may be combined using logical operators (\c and,\c or) and !> parentheses. !> !> \par Example !> \code 1 and (group2 or group3) and y < 1 \endcode !> will select boundary faces !> of color 1, belonging to groups 'group2' or 'group3' and with face center !> coordinate y less than 1. !> !> Operators priority, from highest to lowest: !> '( )' > 'not' > 'and' > 'or' > 'xor' !> !> Similarly, interior faces and cells can be identified using the \ref getfac !> and \ref getcel subroutines (respectively). Their syntax are identical to !> \ref getfbr syntax. !> !> For a more thorough description of the criteria syntax, see the user guide. !> !> !> \section bc_types Boundary condition types !> !> Boundary conditions may be assigned in two ways. !> !> !> \subsection std_bcs For "standard" boundary conditions: !> !> One defines a code in the \c itypfb !> array (of dimensions number of boundary faces). !> This code will then be used by a non-user subroutine to assign the !> following conditions. !> The available codes are: !> - \c ientre: Inlet !> - \c isolib: Free outlet !> - \c isymet: Symmetry !> - \c iparoi: Wall (smooth) !> - \c iparug: Rough wall !> !> These integers are defined elsewhere (in paramx.f90 module). !> Their value is greater than or equal to 1 and less than or equal to !> ntypmx (value fixed in paramx.h) !> !> In addition, some values must be defined: !> - Inlet (more precisely, inlet/outlet with prescribed flow, as the flow !> may be prescribed as an outflow): !> - Dirichlet conditions on variables other than pressure are mandatory !> if the flow is incoming, optional if the flow is outgoing (the code !> assigns zero flux if no Dirichlet is specified); thus, !> at face \c ifac, for the variable \c ivar: \c rcodcl(ifac, ivar, 1) !> !> !> - Smooth wall: (= impermeable solid, with smooth friction) !> - Velocity value for sliding wall if applicable: !> - \c rcodcl(ifac, iu, 1) = fluid velocity in the x direction !> - \c rcodcl(ifac, iv, 1) = fluid velocity in the y direction !> - \c rcodcl(ifac, iw, 1) = fluid velocity in the z direction !> - Specific code and prescribed temperature value at wall if applicable: !> - \c icodcl(ifac, ivar) = 5 !> - \c rcodcl(ifac, ivar, 1) = prescribed temperature !> - Specific code and prescribed flux value at wall if applicable: !> - \c icodcl(ifac, ivar) = 3 !> - \c rcodcl(ifac, ivar, 3) = prescribed flux !> . !> Note that the default condition for scalars (other than k and epsilon) !> is homogeneous Neumann. !> !> !> - Rough wall: (= impermeable solid, with rough friction) !> - Velocity value for sliding wall if applicable: !> - \c rcodcl(ifac, iu, 1) = fluid velocity in the x direction !> - \c rcodcl(ifac, iv, 1) = fluid velocity in the y direction !> - \c rcodcl(ifac, iw, 1) = fluid velocity in the z direction !> - Value of the dynamic roughness height to specify in !> - \c rcodcl(ifac, iu, 3) !> - Value of the scalar roughness height (if required) to specify in !> - \c rcodcl(ifac, iv, 3) (values for iw are not used) !> - Specific code and prescribed temperature value at wall if applicable: !> - \c icodcl(ifac, ivar) = 6 !> - \c rcodcl(ifac, ivar, 1) = prescribed temperature !> - Specific code and prescribed flux value at rough wall, if applicable: !> - \c icodcl(ifac, ivar) = 3 !> - \c rcodcl(ifac, ivar, 3) = prescribed flux !> . !> Note that the default condition for scalars (other than k and epsilon) !> is homogeneous Neumann. !> !> - Symmetry (= slip wall): !> - Nothing to specify !> !> - Free outlet (more precisely free inlet/outlet with prescribed pressure) !> - Nothing to prescribe for pressure and velocity. For scalars and !> turbulent values, a Dirichlet value may optionally be specified. !> The behavior is as follows: !> - pressure is always handled as a Dirichlet condition !> - if the mass flow is inflowing: !> one retains the velocity at infinity !> Dirichlet condition for scalars and turbulent values !> (or zero flux if the user has not specified a !> Dirichlet value) !> - if the mass flow is outflowing: !> one prescribes zero flux on the velocity, the scalars, !> and turbulent values !> . !> Note that the pressure will be reset to p0 on the first free outlet !> face found. !> !> !> \subsection nonstd_bcs For "non-standard" conditions: !> !> Other than (inlet, free outlet, wall, symmetry), one defines !> - on one hand, for each face: !> - an admissible \c itypfb value (i.e. greater than or equal to 1 and !> less than or equal to \c ntypmx; see its value in paramx.h). !> The values predefined in paramx.h: !> \c ientre, \c isolib, \c isymet, \c iparoi, \c iparug are in this range, !> and it is preferable not to assign one of these integers to \c itypfb !> randomly or in an inconsiderate manner. To avoid this, one may use !> \c iindef if one wish to avoid checking values in paramx.h. \c iindef !> is an admissible value to which no predefined boundary condition !> is attached. !> Note that the \c itypfb array is reinitialized at each time step to !> the non-admissible value of 0. If one forgets to modify \c itypfb for !> a given face, the code will stop. !> !> - and on the other hand, for each face and each variable: !> - a code !> - \c icodcl(ifac, ivar) !> - three real values !> - \c rcodcl(ifac, ivar, 1) !> - \c rcodcl(ifac, ivar, 2) !> - \c rcodcl(ifac, ivar, 3) !> !> \anchor icodcl \anchor rcodcl !> The value of \c icodcl is taken from the following: !> - 1: Dirichlet (usable for any variable) !> - 3: Neumann (usable for any variable) !> - 4: Symmetry (usable only for the velocity and components of !> the Rij tensor) !> - 5: Smooth wall (usable for any variable except for pressure) !> - 6: Rough wall (usable for any variable except for pressure) !> - 9: Free outlet (usable only for velocity) !> - 13: Dirichlet for the advection operator and !> Neumann for the diffusion operator !> !> The values of the 3 \c rcodcl components are: !> - \c rcodcl(ifac, ivar, 1): !> - Dirichlet for the variable if \c icodcl(ifac, ivar) = 1 or 13 !> - Wall value (sliding velocity, temp) if \c icodcl(ifac, ivar) = 5 !> . !> The dimension of \c rcodcl(ifac, ivar, 1) is that of the !> resolved variable, for instance: !> - U (velocity in m/s), !> - T (temperature in degrees) !> - H (enthalpy in J/kg) !> - F (passive scalar in -) !> - \c rcodcl(ifac, ivar, 2): !> "exterior" exchange coefficient (between the prescribed value !> and the value at the domain boundary) !> rinfin = infinite by default !> - For velocities U, in kg/(m2 s): !> \c rcodcl(ifac, ivar, 2) = (viscl+visct) / d !> - For the pressure P, in s/m: !> \c rcodcl(ifac, ivar, 2) = dt / d !> - For temperatures T, in Watt/(m2 degres): !> \c rcodcl(ifac, ivar, 2) = Cp*(viscls+visct/turb_schmidt) / d !> - For enthalpies H, in kg /(m2 s): !> \c rcodcl(ifac, ivar, 2) = (viscls+visct/turb_schmidt) / d !> - For other scalars F in: !> \c rcodcl(ifac, ivar, 2) = (viscls+visct/turb_schmidt) / d !> (d has the dimension of a distance in m) !> !> - \c rcodcl(ifac, ivar, 3) if \c icodcl(ifac, ivar) = 3 or 13: !> Flux density (< 0 if gain, n outwards-facing normal) !> - For velocities U, in kg/(m s2) = J: !> \c rcodcl(ifac, ivar, 3) = -(viscl+visct) * (grad U).n !> - For pressure P, in kg/(m2 s): !> \c rcodcl(ifac, ivar, 3) = -dt * (grad P).n !> - For temperatures T, in Watt/m2: !> \c rcodcl(ifac, ivar, 3) = -Cp*(viscls+visct/turb_schmidt) * (grad T).n !> - For enthalpies H, in Watt/m2: !> \c rcodcl(ifac, ivar, 3) = -(viscls+visct/turb_schmidt) * (grad H).n !> - For other scalars F in: !> \c rcodcl(ifac, ivar, 3) = -(viscls+visct/turb_schmidt) * (grad F).n !> !> - \c rcodcl(ifac, ivar, 3) if \c icodcl(ifac, ivar) = 6: !> Roughness for the rough wall law !> - For velocities U, dynamic roughness !> \c rcodcl(ifac, iu, 3) = roughd !> - For other scalars, thermal roughness !> \c rcodcl(ifac, iv, 3) = rough !> !> !> Note that if the user assigns a value to \c itypfb equal to \c ientre, \c isolib, !> \c isymet, \c iparoi, or \c iparug and does not modify \c icodcl (zero value by !> default), \c itypfb will define the boundary condition type. !> !> To the contrary, if the user prescribes \c icodcl(ifac, ivar) (nonzero), !> the values assigned to \c rcodcl will be used for the considered face !> and variable (if \c rcodcl values are not set, the default values will !> be used for the face and variable, so: !> - \c rcodcl(ifac, ivar, 1) = 0.d0 !> - \c rcodcl(ifac, ivar, 2) = rinfin !> - \c rcodcl(ifac, ivar, 3) = 0.d0) !> !> Especially, one may have for example: !> - set \c itypfb(ifac) = \c iparoi which prescribes default wall !> conditions for all variables at face ifac, !> - and define IN ADDITION for variable ivar on this face specific !> conditions by specifying \c icodcl(ifac, ivar) and the 3 \c rcodcl values. !> !> The user may also assign to \c itypfb a value not equal to \c ientre, \c isolib, !> \c isymet, \c iparoi, \c iparug, \c iindef but greater than or equal to 1 and less !> than or equal to ntypmx (see values in param.h) to distinguish groups !> or colors in other subroutines which are specific to the case and in !> which itypfb is accessible. In this case though it will be necessary !> to prescribe boundary conditions by assigning values to icodcl and to !> the 3 \c rcodcl fields (as the value of \c itypfb will not be predefined in !> the code). !> !> !> \subsection comp_bcs Boundary condition types for compressible flows !> !> For compressible flows, only predefined boundary conditions may !> be assigned among: \c iparoi, \c isymet, \c iesicf, \c isspcf, \c isopcf, \c iephcf, \c ieqhcf !> !> - \c iparoi : standard wall !> - \c isymet : standard symmetry !> !> - \c iesicf, \c isspcf, \c isopcf, \c iephcf, \c ieqhcf : inlet/outlet !> !> For inlets/outlets, we can prescribe !> a value for turbulence and passive scalars in \c rcodcl(.,.,1) !> for the case in which the mass flux is incoming. If this is not !> done, a zero flux condition is applied. !> !> - \c iesicf: prescribed inlet/outlet (for example supersonic inlet) !> the user prescribes the velocity and all thermodynamic variables !> - \c isspcf: supersonic outlet !> the user does not prescribe anything !> - \c isopcf: subsonic outlet with prescribed pressure !> the user presribes the pressure !> - \c iephcf: mixed inlet with prescribed total pressure and enthalpy !> the user prescribes the total pressure and total enthalpy !> - \c ieqhcf: subsonic inlet with prescribed mass and enthalpy flow !> to be implemented !> !> !> \subsection cons_rul Consistency rules !> !> A few consistency rules between \c icodcl codes for variables with !> non-standard boundary conditions: !> !> - Codes for velocity components must be identical !> - Codes for Rij components must be identical !> - If code (velocity or Rij) = 4 !> one must have code (velocity and Rij) = 4 !> - If code (velocity or turbulence) = 5 !> one must have code (velocity and turbulence) = 5 !> - If code (velocity or turbulence) = 6 !> one must have code (velocity and turbulence) = 6 !> - If scalar code (except pressure or fluctuations) = 5 !> one must have velocity code = 5 !> - If scalar code (except pressure or fluctuations) = 6 !> one must have velocity code = 6 !> !> !> \remarks !> - Caution: to prescribe a flux (nonzero) to Rij, the viscosity to take !> into account is viscl even if visct exists !> (visct=rho cmu k2/epsilon) !> - One have the ordering array for boundary faces from the previous time !> step (except for the fist one, where \c itrifb has not been set yet). !> - The array of boundary face types \c itypfb has been reset before !> entering the subroutine. !> !> !> \subsubsection cs_user_bc_cell_id Cell values of some variables !> !> Cell value field ids !> !> - Density: \c irom !> - Dynamic molecular viscosity: \c iviscl !> - Turbulent viscosity: \c ivisct !> - Specific heat: \c icp !> - Diffusivity(lambda): \c field_get_key_int(ivarfl(isca(iscal)), & !> kivisl, ...) !> !> !> \subsubsection fac_id Faces identification !> !> - Density: \c field id \c ibrom !> - Boundary mass flux (for convecting \c ivar): !> field id \c iflmab !> using \c field_get_key_int(ivarfl(ivar), kbmasf, iflmab) !> - For other values: take as an approximation the value in the adjacent cell !> i.e. as above with \c iel = ifabor(ifac). !> !> Please refer to the !> boundary conditions !> section of the theory guide for more informations. !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! Arguments !______________________________________________________________________________. ! mode name role ! !______________________________________________________________________________! !> \param[in] nvar total number of variables !> \param[in] nscal total number of scalars !> \param[out] icodcl boundary condition code: !> - 1 Dirichlet !> - 2 Radiative outlet !> - 3 Neumann !> - 4 sliding and !> \f$ \vect{u} \cdot \vect{n} = 0 \f$ !> - 5 smooth wall and !> \f$ \vect{u} \cdot \vect{n} = 0 \f$ !> - 6 rough wall and !> \f$ \vect{u} \cdot \vect{n} = 0 \f$ !> - 9 free inlet/outlet !> (input mass flux blocked to 0) !> - 13 Dirichlet for the advection operator and !> Neumann for the diffusion operator !> \param[in] itrifb indirection for boundary faces ordering !> \param[in,out] itypfb boundary face types !> \param[out] izfppp boundary face zone number !> \param[in] dt time step (per cell) !> \param[in,out] rcodcl boundary condition values: !> - rcodcl(1) value of the dirichlet !> - rcodcl(2) value of the exterior exchange !> coefficient (infinite if no exchange) !> - rcodcl(3) value flux density !> (negative if gain) in w/m2 or roughness !> in m if icodcl=6 !> -# for the velocity \f$ (\mu+\mu_T) !> \gradt \, \vect{u} \cdot \vect{n} \f$ !> -# for the pressure \f$ \Delta t !> \grad P \cdot \vect{n} \f$ !> -# for a scalar \f$ cp \left( K + !> \dfrac{K_T}{\sigma_T} \right) !> \grad T \cdot \vect{n} \f$ !_______________________________________________________________________________ subroutine cs_f_user_boundary_conditions & ( nvar , nscal , & icodcl , itrifb , itypfb , izfppp , & dt , & rcodcl ) !=============================================================================== !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use optcal use cstphy use cstnum use entsor use parall use period use ihmpre use ppppar use ppthch use coincl use cpincl use ppincl use ppcpfu use atincl use atsoil use ctincl use cs_fuel_incl use mesh use field use turbomachinery use iso_c_binding use cs_c_bindings !=============================================================================== implicit none ! Arguments integer nvar , nscal integer icodcl(nfabor,nvarcl) integer itrifb(nfabor), itypfb(nfabor) integer izfppp(nfabor) double precision dt(ncelet) double precision rcodcl(nfabor,nvarcl,3) ! Local variables !< [loc_var_dec] integer ifac, ii, iel,iproc,icell integer idim,izone,it_pvbc,ilelt, nlelt,int32size,stat!,fbrnum !cellsnum integer ifcpvbc,ifcpvbcr,ifcur!x,ifcury,ifcurz !pointers integer ipotr, ipoti, f_id, ipotva,keyvar integer ndimve parameter (ndimve = 3) !dimensions number integer ncelglob4 integer icelcur(ncelgb) !parameter (fbrnum = 6390)!nfbrgb) !parameter (cellsnum = 100000)!ncelgb - global ncel) integer, allocatable, dimension(:) :: lstelt character(len=80) :: filename,format_string,filenameout !double precision, dimension(:), pointer :: cpro_cury,cpro_curz double precision, dimension(:,:), pointer :: cpro_cur,cpro_pvbc,cpro_pvbcrel double precision elcurr , elcurr1, surface double precision ray , Rarc , jmax , bb, mu0, pii,mu0div4pi,currelcoef double precision rr1, rr2, debit, surftot, inletvel double precision vol(ncelgb),coord(3,ncelgb),elcur(3,ncelgb) double precision elcurcheck,pvrelcoef,rfminrc,elcursum,maxBC,checksum,pvbc_mag !< [loc_var_dec] common /curoutput/ filenameout !=============================================================================== allocate(lstelt(nfabor)) ! temporary array for boundary faces selection !=============================================================================== ! Assign boundary conditions to boundary faces here ! For each subset: ! - use selection criteria to filter boundary faces of a given subset ! - loop on faces from a subset ! - set the boundary condition for each face !=============================================================================== ! --- For boundary faces of entree assign an inlet for all phases and assign a cathode for "electric" variables. !========================================================================================================= !This subroutine doesn't work for gas mixtures !************************************************** ! CONSTANTS * !************************************************** Rarc = 1.5D-3 jmax = 1.27D8 elcurr = 2.0D1 elcurr1 = 0.0D0 bb = 1.96D3 jmax = 6*elcurr / (pi*Rarc*Rarc) elcurr = 0.0D0 mu0 = 1.2566370614d-6 !Vacuum permeability pii = 3.14159265359d0 !π number mu0div4pi=mu0/4d0/pii it_pvbc=50 !PV BC recalculation period pvrelcoef=0.99d0 !Relaxation coefficient for vector potential boundary condition int32size=2147483647 !Inlet Velocity rr2 = 3.0D-3 !external radius rr1 = 1.5D-3 !internal radius debit = 0.11D-3 surftot = pi * ((rr2**2)-(rr1**2)) inletvel = debit/surftot if ( ippmod(ielarc).ge.2 ) then !ONLY FOR ELECTRIC ARCS !********************************************************************************************* ! POINTERS * !********************************************************************************************* call field_get_key_id("variable_id", keyvar) call field_get_id('elec_pot_r', f_id) call field_get_key_int(f_id, keyvar, ipotr) call field_get_id('vec_potential', f_id) call field_get_key_int(f_id, keyvar, ipotva) !ifcpvbc=-1 !ifcpvbcr=-1 call field_get_id_try('PotVecBC', ifcpvbc) !Obtaining pointer to PV boundary condition directly calculated from electric current. if (ifcpvbc.ge.0) call field_get_val_v(ifcpvbc, cpro_pvbc) !obtaining PV boundary condition directly calculated from electric current. if (ifcpvbc.lt.0) print*,'ERROR: No PotVecBC id. cs_user_boundary_conditions.f90' !error message just in case. call field_get_id_try('PVBCRel', ifcpvbcr) !Obtaining pointer to PV boundary condition with relaxation. if (ifcpvbcr.ge.0) call field_get_val_v(ifcpvbcr, cpro_pvbcrel) !obtaining PV boundary condition with relaxation if (ifcpvbcr.lt.0) print*,'ERROR: No PVBCRel id. cs_user_boundary_conditions.f90' !error message just in case. call field_get_id_try('current_re', ifcur) if (ifcur.ge.0) call field_get_val_v(ifcur, cpro_cur) if (ifcur.lt.0) print*,'ERROR: No current_re id. cs_user_boundary_conditions.f90' !error message just in case. !call field_get_id_try('current_re_2', ifcury) !ifcury=ifcurx+1 !if (ifcury.ge.0) call field_get_val_s(ifcury, cpro_cury) !if (ifcury.lt.0) print*,'ERROR: No current_re_2 id. cs_user_boundary_conditions.f90' !error message just in case. !call field_get_id_try('current_re_3', ifcurz) !ifcurz=ifcury+1 !if (ifcurz.ge.0) call field_get_val_s(ifcurz, cpro_curz) !if (ifcurz.lt.0) print*,'ERROR: No current_re_3 id. cs_user_boundary_conditions.f90' !error message just in case. !*************************************************************************************** ! INITIALIZATION OF PV BC WITH ZERO AT FIRST TIME STEP * !*************************************************************************************** if (ntcabs-ntpabs.eq.1 .and. ifcpvbcr.ge.0) then call getfbr('all[]', nlelt, lstelt) write(*,'(A,I5)')'global number of boundary faces nfbrgb=',nfbrgb do ilelt = 1, nlelt ifac = lstelt(ilelt) do idim = 1, ndimve cpro_pvbcrel(ifac,idim)=0d0 end do end do end if !*************************************************************************************** ! ELECTRIC CURRENT OUTPUT FOR PARALLEL CASE FOR BIG MESHES * !*************************************************************************************** if (irangp.ge.0 .and. ncelgb.gt.int32size .and. ifcpvbc.ge.0 .and. ifcpvbcr.ge.0) then !STEP 1. Defining of file name for current thread = "geom_eleccur_+(1..(nrangp-1))+.out" !nrangp - number of processes (=1 if sequental) !irangp - process rank r (0 < r < n_processes) in distributed parallel run if (ntcabs-ntpabs.eq.1) then !write(*,'(A, if (irangp < 10) then format_string = "(A13,I1,A4)" else if(irangp <100) then format_string = "(A13,I2,A4)" else format_string = "(A13,I3,A4)" end if write (filenameout,format_string) "geom_eleccur_", irangp,".out" !inquire(FILE=filenameout, exist=file_exists) open(unit=impusr(10), iostat=stat, file=filenameout, status='new',form='formatted') close(impusr(10)) end if !STEP 2. Cleaning of the file from old data. if (mod(ntcabs-ntpabs+1,it_pvbc).eq.0) then open(unit=impusr(10), iostat=stat, file=filenameout, status='replace',form='formatted') close(impusr(10)) end if !STEP 3. The file writing. if ((mod(ntcabs-ntpabs,it_pvbc).eq.0 .and. ntcabs-ntpabs.ge.it_pvbc) .or. & ((ntcabs-ntpabs).eq.2)) then write(*,'(A20,A,I6)')filenameout,'Output of electric current and geometry. Time step ',ntcabs open(unit=impusr(10), iostat=stat, file=filenameout, status='old',form='formatted') do iel=1,ncel write(impusr(10),'(I6,E18.8,E18.8,E18.8,E18.8,E18.8,E18.8,E18.8)') & iel,volume(iel),xyzcen(1,iel),xyzcen(2,iel),xyzcen(3,iel), & cpro_cur(1,iel),cpro_cur(2,iel),cpro_cur(3,iel) end do close(impusr(10)) !attempt to equalize the progress for every processor checksum=1d0 call parsom(checksum) write(*,'(I3,A,F10.3)')irangp,' checksum=',checksum end if end if !*************************************************************************************** ! BOUNDARY CONDION FOR VECTOR POTENTIAL EVERY it_pvbc time steps * !*************************************************************************************** if (ifcpvbc.ge.0 .and. ifcpvbcr.ge.0) then ! Recalculation starts at time step number 2,1*it_pvbc,2*it_pvbc,.. ! Vector Potential is calculated with the formula similar to the Biot Savart law: ! cpro_pvbc(iface,1) = Ax= mu0/4π *∫∫∫jx*dV/|vector Rface - vector Rcell| ! cpro_pvbc(iface,2) = Ay= mu0/4π *∫∫∫jy*dV/|vector Rface - vector Rcell| ! cpro_pvbc(iface,3) = Az= mu0/4π *∫∫∫jz*dV/|vector Rface - vector Rcell| if ((mod(ntcabs-ntpabs,it_pvbc).eq.0 .and. ntcabs-ntpabs.ge.it_pvbc) .or. & (ntcabs-ntpabs.eq.2)) then !SINGLE PROCESSOR if (irangp.lt.0) then call getfbr('all[]', nlelt, lstelt) write(*,'(I3,A64,I3)')irangp,' Boundary condition for vector potential calculation. Time step=',ntcabs do ilelt = 1, nlelt ifac = lstelt(ilelt) cpro_pvbc(ifac,1)=0d0 cpro_pvbc(ifac,2)=0d0 cpro_pvbc(ifac,3)=0d0 do iel=1,ncel !Calculation of |vector Rface - vector Rcell| rfminrc=sqrt((cdgfbo(1,ifac)-xyzcen(1,iel))**2+ & (cdgfbo(2,ifac)-xyzcen(2,iel))**2+ & (cdgfbo(3,ifac)-xyzcen(3,iel))**2) !Calculation of integral ∫∫∫j*dV/|vector Rface - vector Rcell| cpro_pvbc(ifac,1)=cpro_pvbc(ifac,1)+cpro_cur(1,iel)*volume(iel)/rfminrc cpro_pvbc(ifac,2)=cpro_pvbc(ifac,2)+cpro_cur(2,iel)*volume(iel)/rfminrc cpro_pvbc(ifac,3)=cpro_pvbc(ifac,3)+cpro_cur(3,iel)*volume(iel)/rfminrc end do cpro_pvbc(ifac,1)=cpro_pvbc(ifac,1)*mu0div4pi cpro_pvbc(ifac,2)=cpro_pvbc(ifac,2)*mu0div4pi cpro_pvbc(ifac,3)=cpro_pvbc(ifac,3)*mu0div4pi end do end if !PARALLEL CALCULATION !SMALL MESHES - if ncelgb value fits to __int32 - MPI restriction !integer(kind=4) values range from -2,147,483,648 to 2,147,483,647 if (irangp.ge.0 .and. ncelgb.le.int32size) then ncelglob4 = ncelgb call paragv(ncel, ncelglob4, xyzcen(1,:), coord(1,:)) !MPI paragv(__int32 size of the local array,__int32 size of the global array,array,g_array) call paragv(ncel, ncelglob4, xyzcen(2,:), coord(2,:)) call paragv(ncel, ncelglob4, xyzcen(3,:), coord(3,:)) call paragv(ncel, ncelglob4, cpro_cur(1,:), elcur(1,:)) call paragv(ncel, ncelglob4, cpro_cur(2,:), elcur(2,:)) call paragv(ncel, ncelglob4, cpro_cur(3,:), elcur(3,:)) call paragv(ncel, ncelglob4, volume, vol) call getfbr('all[]', nlelt, lstelt) do ilelt = 1, nlelt ifac = lstelt(ilelt) !Here we don't use loop "do idim = 1,ndimve" to avoid calculation of the same rfminrc three times for every dimension cpro_pvbc(ifac,1)=0d0 cpro_pvbc(ifac,2)=0d0 cpro_pvbc(ifac,3)=0d0 do icell=1,ncelgb !Calculation of |vector Rface - vector Rcell| rfminrc=sqrt((cdgfbo(1,ifac)-coord(1,icell))**2 & +(cdgfbo(2,ifac)-coord(2,icell))**2 & +(cdgfbo(3,ifac)-coord(3,icell))**2) !Calculation of integral j*dV/|vector Rface - vector Rcell| cpro_pvbc(ifac,1)=cpro_pvbc(ifac,1)+elcur(1,icell)*vol(icell)/rfminrc cpro_pvbc(ifac,2)=cpro_pvbc(ifac,2)+elcur(2,icell)*vol(icell)/rfminrc cpro_pvbc(ifac,3)=cpro_pvbc(ifac,3)+elcur(3,icell)*vol(icell)/rfminrc end do cpro_pvbc(ifac,1)=cpro_pvbc(ifac,1)*mu0div4pi cpro_pvbc(ifac,2)=cpro_pvbc(ifac,2)*mu0div4pi cpro_pvbc(ifac,3)=cpro_pvbc(ifac,3)*mu0div4pi !Calculation of Max PV BC for output and debugging. pvbc_mag=sqrt(cpro_pvbc(ifac,1)**2 + cpro_pvbc(ifac,2)**2 + cpro_pvbc(ifac,3)**2) if (pvbc_mag.gt.maxBC) maxBC=pvbc_mag end do write(*,'(I3,I6,A,E12.5)')irangp,ntcabs,' MPI BC for vector potential calculation. Max Value=',maxBC end if !BIG MESHES with ncelgb bigger than __int32 range !BETA VERSION. ALGORITHM IS TESTED ONLY WITH SMALL MESH WITH 93375 CELLS !ncelgb = ncelg2 in majgeo ncelg2 = mesh->n_g_cells = n_g_elts[0]; in cs_preprocess.c from MPI_Allreduce if (irangp.ge.0 .and. ncelgb.gt.int32size) then write(*,'(I3,A,I5)')irangp,' files are read t=',ntcabs iel=1 elcurcheck=0d0 elcursum = 0d0 do iproc=0,nrangp-1 if (iproc < 10) then format_string = "(A13,I1,A4)" else if(iproc <100) then format_string = "(A13,I2,A4)" else format_string = "(A13,I3,A4)" end if write (filename,format_string) "geom_eleccur_", iproc,".out" open(unit=impusr(12), file=filename, status='old',form='formatted') do while (.true.) !icelcur read(impusr(12),*,end=999) icelcur(iel),vol(iel),coord(1,iel),coord(2,iel),coord(3,iel), & elcur(1,iel),elcur(2,iel),elcur(3,iel) if(iproc.eq.irangp) then elcurcheck=elcurcheck+dabs(cpro_cur(1,icelcur(iel))-elcur(1,iel))+ & dabs(cpro_cur(2,icelcur(iel))-elcur(2,iel))+ & dabs(cpro_cur(3,icelcur(iel))-elcur(3,iel)) elcursum = elcursum + dabs(cpro_cur(1,icelcur(iel)))+ & dabs(cpro_cur(2,icelcur(iel)))+ & dabs(cpro_cur(3,icelcur(iel))) !if (mod(iel,10000).eq.0) then !output for debugging ! write(*,'(I3,I6,E17.7,E17.7,E17.7,E17.7,E17.7,E17.7,E17.7)') & ! irangp,iel,vol(iel),coord(1,iel),coord(2,iel), & ! coord(3,iel),elcur(1,iel),elcur(2,iel),elcur(3,iel) ! write(*,'(I3,I6,E17.7,E17.7,E17.7,E17.7,E17.7,E17.7,E17.7)') & ! irangp,icelcur(iel),volume(icelcur(iel)),xyzcen(1,icelcur(iel)),xyzcen(2,icelcur(iel)), & ! xyzcen(3,icelcur(iel)),cpro_cur(1,icelcur(iel)),cpro_cur(2,icelcur(iel)),cpro_cur(3,icelcur(iel)) !end if end if iel=iel+1 end do 999 continue close(impusr(12)) end do if (elcursum.gt.0d0) then write(*,'(I3,A,E15.7,I6,I6)')irangp,' Electric current from file delta/Sum(j)=',elcurcheck/elcursum maxBC=0d0 call getfbr('all[]', nlelt, lstelt) do ilelt = 1, nlelt ifac = lstelt(ilelt) !Here we don't use loop "do idim = 1,ndimve" to avoid calculation of the same rfminrc three times for every dimension cpro_pvbc(ifac,1)=0d0 cpro_pvbc(ifac,2)=0d0 cpro_pvbc(ifac,3)=0d0 do icell=1,ncelgb !Calculation of |vector Rface - vector Rcell| rfminrc=sqrt((cdgfbo(1,ifac)-coord(1,icell))**2 & +(cdgfbo(2,ifac)-coord(2,icell))**2 & +(cdgfbo(3,ifac)-coord(3,icell))**2) !Calculation of integral j*dV/|vector Rface - vector Rcell| cpro_pvbc(ifac,1)=cpro_pvbc(ifac,1)+elcur(1,icell)*vol(icell)/rfminrc cpro_pvbc(ifac,2)=cpro_pvbc(ifac,2)+elcur(2,icell)*vol(icell)/rfminrc cpro_pvbc(ifac,3)=cpro_pvbc(ifac,3)+elcur(3,icell)*vol(icell)/rfminrc end do cpro_pvbc(ifac,1)=cpro_pvbc(ifac,1)*mu0div4pi cpro_pvbc(ifac,2)=cpro_pvbc(ifac,2)*mu0div4pi cpro_pvbc(ifac,3)=cpro_pvbc(ifac,3)*mu0div4pi !Calculation of Max PV BC for output and debugging. pvbc_mag=sqrt(cpro_pvbc(ifac,1)**2 + cpro_pvbc(ifac,2)**2 + cpro_pvbc(ifac,3)**2) if (pvbc_mag.gt.maxBC) maxBC=pvbc_mag end do write(*,'(I3,A48,E12.5,A,I5)')irangp,' BC for vector potential calculation. Max Value=',maxBC,' nlelt=',nlelt else write(*,'(I3,A,E15.7,I6,I6)')irangp,'ERROR: Sum(j)=0d0 Electric current from file delta=',elcurcheck end if end if end if !RELAXATION OF BOUNDARY CONDITION FOR VECTOR POTENTIAL AT EVERY TIME STEP call getfbr('all[]', nlelt, lstelt) do ilelt = 1, nlelt ifac = lstelt(ilelt) do idim = 1, ndimve cpro_pvbcrel(ifac,idim)=cpro_pvbcrel(ifac,idim)*pvrelcoef+(1-pvrelcoef)*cpro_pvbc(ifac,idim) icodcl(ifac,ipotva + idim-1) = 1 rcodcl(ifac,ipotva + idim-1,1) = cpro_pvbcrel(ifac,idim) end do end do !If there is no field for vector potential boundary condition we impose PVNV for bord_lateral and PVNF for all other faces else if (ntcabs-ntpabs.eq.2) print*,"ERROR: No pointer to PV BC or PV BC with relaxation. 2." call getfbr('entree', nlelt, lstelt) do ilelt = 1, nlelt ifac = lstelt(ilelt) do idim = 1, ndimve icodcl(ifac,ipotva + idim-1) = 3 rcodcl(ifac,ipotva + idim-1,3) = 0.d0 end do end do call getfbr('bord_lateral', nlelt, lstelt) do ilelt = 1, nlelt ifac = lstelt(ilelt) do idim = 1, ndimve icodcl(ifac,ipotva + idim-1) = 1 !only here should be A=0 according to Isabelle Choquet 2012 On the choice of electromagnetic model for arcs rcodcl(ifac,ipotva + idim-1,1) = 0.d0 end do end do call getfbr('bord_superieur', nlelt, lstelt) do ilelt = 1, nlelt ifac = lstelt(ilelt) do idim = 1, ndimve icodcl(ifac,ipotva + idim-1) = 3 rcodcl(ifac,ipotva + idim-1,3) = 0.d0 end do end do call getfbr('bord_lateral_nozzle', nlelt, lstelt) do ilelt = 1, nlelt ifac = lstelt(ilelt) do idim = 1, ndimve icodcl(ifac,ipotva + idim-1) = 3 rcodcl(ifac,ipotva + idim-1,3) = 0.d0 end do end do call getfbr('bord_anode', nlelt, lstelt) do ilelt = 1, nlelt ifac = lstelt(ilelt) do idim = 1, ndimve icodcl(ifac,ipotva + idim-1) = 3 rcodcl(ifac,ipotva + idim-1,3) = 0.d0 end do end do call getfbr('base_anode', nlelt, lstelt) do ilelt = 1, nlelt ifac = lstelt(ilelt) do idim = 1, ndimve icodcl(ifac,ipotva + idim-1) = 3 rcodcl(ifac,ipotva + idim-1,3) = 0.d0 end do end do call getfbr('cathode', nlelt, lstelt) do ilelt = 1, nlelt ifac = lstelt(ilelt) do idim = 1, ndimve icodcl(ifac,ipotva + idim-1) = 3 rcodcl(ifac,ipotva + idim-1,3) = 0.d0 end do end do end if end if !*************************************************************************************** ! THERMAL, INLET AND ELECTRIC BOUNDARY CONDITIONS * !*************************************************************************************** call getfbr('entree', nlelt, lstelt) if (ntcabs-ntpabs.eq.1) write(*,'(A,I5)')'BC entree nlelt=',nlelt do ilelt = 1, nlelt ifac = lstelt(ilelt) itypfb(ifac) = ientre izone = 1 ! Zone Number (from 1 to n) ! Zone localization for a given face izfppp(ifac) = izone rcodcl(ifac,iu,1) = 0.0d0 rcodcl(ifac,iv,1) = 0.0d0 rcodcl(ifac,iw,1) = inletvel !Enthalpy ii = iscalt icodcl(ifac,isca(ii)) = 1 rcodcl(ifac,isca(ii),1) = 14.0D3 !corresponds to the temperature of 300 Kelvin for argon at atmospheric pressure (see dp_ELE) !Electric potential if (ippmod(ielarc).ge. 2) then icodcl(ifac,ipotr) = 3 rcodcl(ifac,ipotr,3) = 0.0 end if end do write(nfecra,1011)rcodcl(ifac,iw,1) ! --- For boundary faces of bord_lateral and bord_superieur assign an free outlet for all phases and example of electrode for Joule Effect by direct conduction. !isolib - free outlet face or more precisely free inlet/outlet with forced pressure. DOESN'T WORK IN CS4.2. GIVES UNSTABLE PRESSURE WHICH RESULTS IN UNSTABLE ELECTRIC ARC !ifrent - free outlet, free inlet (based on Bernoulli relationship) face. ! ================================================================================================================================ call getfbr('bord_lateral', nlelt, lstelt) if (ntcabs-ntpabs.eq.1) write(*,'(A,I5)')'BC bord_lateral nlelt=',nlelt do ilelt = 1, nlelt ifac = lstelt(ilelt) itypfb(ifac) = ifrent !isolib is unstable with CS4.2, but was used with CS1.3 and worked okay. izone = 2 ! Zone Number (from 1 to n) izfppp(ifac) = izone ! Zone location for a given face ! --- Handle Scalars ! Enthalpy in J/kg (By default zero flux with ISOLIB). Nothing to do ! Mass fraction of the (n-1) gas mixture components (Zero flux by defaut with ISOLIB). Nothing to do ! Specific model for Joule Effect by direct conduction: ! If you want to make a simulation with an imposed Power PUISIM ! (you want to get PUISIM imposed in useli1 and PUISIM = Amp x Volt) ! you need to impose IELCOR=1 in useli1 ! The boundary conditions will be scaled by COEJOU coefficient ! for example the electrical potential will be multiplied bu COEJOU ! (both real and imaginary part of the electrical potential if needed) ! COEJOU is automatically defined in order that the calculated dissipated power ! by Joule effect (both real and imaginary part if needed) is equal to PUISIM ! At the beginning of the calculation, COEJOU ie equal to 1; ! COEJOU is writing and reading in the result files. ! If you don't want to calculate with by scaling, ! you can impose directly the value. !Electric potential if (ippmod(ielarc).ge. 2) then icodcl(ifac,ipotr) = 3 rcodcl(ifac,ipotr,3) = 0.0 end if end do call getfbr('bord_superieur', nlelt, lstelt) if (ntcabs-ntpabs.eq.1) write(*,'(A,I5)')'BC bord_superieur nlelt=',nlelt do ilelt = 1, nlelt ifac = lstelt(ilelt) itypfb(ifac) = ifrent !isolib is unstable with CS4.2, but was used with CS1.3 and worked okay. izone = 3 izfppp(ifac) = izone !Electric potential if (ippmod(ielarc).ge. 2) then icodcl(ifac,ipotr) = 3 rcodcl(ifac,ipotr,3) = 0.0 end if end do call getfbr('bord_lateral_nozzle', nlelt, lstelt) if (ntcabs-ntpabs.eq.1) write(*,'(A,I5)')'BC bord_lateral_nozzle nlelt=',nlelt do ilelt = 1, nlelt ifac = lstelt(ilelt) itypfb(ifac) = iparoi !smooth solid wall face, impermeable and with friction. izone = 4 izfppp(ifac) = izone !Enthalpy ii = iscalt icodcl(ifac,isca(ii)) = 3 rcodcl(ifac,isca(ii),3) = 0.0D0 !zero flux of enthalpy !Electric potential if (ippmod(ielarc).ge. 2) then icodcl(ifac,ipotr) = 3 rcodcl(ifac,ipotr,3) = 0.0 end if end do call getfbr('bord_anode', nlelt, lstelt) if (ntcabs-ntpabs.eq.1) write(*,'(A,I5)')'BC bord_anode nlelt=',nlelt do ilelt = 1, nlelt ifac = lstelt(ilelt) itypfb(ifac) = iparoi izone = 5 izfppp(ifac) = izone !Enthalpy ii = iscalt icodcl(ifac,isca(ii)) = 1 rcodcl(ifac,isca(ii),1) = 14.0D3 !corresponds to the temperature of 300 Kelvin for copper (see H-T_Cu2) !Electric potential if (ippmod(ielarc).ge. 2) then icodcl(ifac,ipotr) = 1 rcodcl(ifac,ipotr,1) = 0.0 end if end do call getfbr('base_anode', nlelt, lstelt) if (ntcabs-ntpabs.eq.1) write(*,'(A,I5)')'BC base_anode nlelt=',nlelt do ilelt = 1, nlelt ifac = lstelt(ilelt) itypfb(ifac) = iparoi izone = 6 izfppp(ifac) = izone !Enthalpy ii = iscalt icodcl(ifac,isca(ii)) = 3 rcodcl(ifac,isca(ii),3) = 0.d0!14.0D3 corresponds to the temperature of 300 Kelvin for copper (see H-T_Cu2) !Electric potential if (ippmod(ielarc).ge. 2) then icodcl(ifac,ipotr) = 1 rcodcl(ifac,ipotr,1) = 0.0 end if end do !For Argon start with 100Amp is unstable. But for 30Amp it is okay. Therefore we use relaxation from 30Amp to 100Amp. !currelcoef=0.99d0 !Relaxation coefficient for the electric current value. !couimp = couimp*currelcoef+100.0d0*(1-currelcoef) !Imposed current intensity (electric arc) in Amp !if (mod(ntcabs-ntpabs,50).eq.0) write(*,'(I5,A,F10.4)')ntcabs,' Imposed current couimp=',couimp call getfbr('cathode', nlelt, lstelt) if (ntcabs-ntpabs.eq.1) write(*,'(A,I5)')'BC cathode nlelt=',nlelt do ilelt = 1, nlelt ifac = lstelt(ilelt) itypfb(ifac) = iparoi izone = 7 izfppp(ifac) = izone !Enthalpy ii = iscalt icodcl(ifac,isca(ii)) = 3 rcodcl(ifac,isca(ii),3) = 0.0D0 !zero flux of enthalpy !Electric potential if (ippmod(ielarc).ge. 2) then icodcl(ifac,ipotr) = 1 rcodcl(ifac,ipotr,1) = -pot_diff !call field_get_val_v(iprpfl(idjr(1)), djr) !zelcurr = djr(3,iel) !call field_get_val_s(iprpfl(idjr(3)), djr3) iel = ifabor(ifac) surface = surfbn(ifac) elcurr = elcurr + cpro_cur(3,iel) * surface !elcou = elcou + cpro_cur(3,iel) * surfac(3,ifac) ray = sqrt((cdgfbo(1,ifac)**2)+ (cdgfbo(2,ifac)**2)) elcurr1 = elcurr1 + (jmax * exp((-1.0D0)*bb*ray)* surface) !I don't know what is this end if end do if (irangp.ge.0) then call parsom (elcurr) call parsom (elcurr1) end if elcurr = ABS(elcurr) write(nfecra,1009)elcurr write(nfecra,1010)elcurr1 !-------- ! Formats !-------- 1009 format(/,'elcurr = ',E14.5) 1010 format(/,'elcurr1= ',E14.5) 1011 format(/,'Inlet Velocity = ',E14.5) ! 1010 format(/,' Nx = ',E14.5,/, ' Ny = ',E14.5,/, ' Nz = ',E14.5) !---- ! End !---- deallocate(lstelt) ! temporary array for boundary faces selection return end subroutine cs_f_user_boundary_conditions