Моделирование теплопроводности в тулиевой фольге крыловидной формы для источника ионов

Полный текст

1. Introduction Current-heated metallic foils are commonly used as evaporation elements in ion sources [1]. Their operating temperature is set by the competition between Joule heating, conduction to the cooled contacts, radiation, and evaporation. Because the evaporation rate is strongly temperature dependent, even a moderate error in the predicted temperature field can lead to a large error in the estimated material flux. We consider a thulium foil with rectangular side wings. The wings enlarge the heated area but also disturb the current paths, particularly at the abrupt changes in width. The resulting problem is nonlinear: density, heat capacity, thermal conductivity, and electrical resistivity all vary with temperature, while radiation and evaporation introduce additional temperature-dependent losses. An earlier formulation ranked candidate geometries by the area lying within a specified temperature interval [2]. For an evaporation source, the surface flux is a more direct performance measure. We therefore compute the stationary temperature field and evaluate the local evaporation rate from a vapor-pressure correlation and the Hertz-Knudsen formula [3-5]. © 2026 A. S. Ayriyan, E. A. Ayryan, E. E. Donets, M. P. Veneva This work is licensed under a Creative Commons “Attribution-NonCommercial 4.0 International” license. 2L = 20.0 mm x calculation domain: x ≥ 0, z ≥ 0 2Wmax = 2.5 mm O z 2Wmin = 1.0 mm 2l = 14.0 mm Thickness: h = 0.10 mm Figure 1. Geometry of the wing-shaped thulium foil in the (x, z) plane. The plate is symmetric about both the x- and z-axes: the wide central region of full width 2Wmax = 2.5 mm and length 2l = 14.0 mm narrows to the terminal contacts of width 2Wmin = 1.0 mm, giving a total length 2L = 20.0 mm and a uniform thickness h = 0.10 mm. Due to the symmetry, the computations are restricted to one quarter of the plate (x ≥ 0, z ≥ 0; shaded region) only, with the origin O at the centre The optimization itself requires many solutions for different wing lengths and currents. This paper addresses the direct problem that underlies those repeated evaluations. We discretize one quarter of the symmetric foil on a masked Cartesian grid and use an alternating-direction implicit (ADI) scheme. Each half-step produces independent tridiagonal systems along rows or columns, which we solve in batches on a CUDA-capable GPU. The calculation returns the stationary field, the maximum foil temperature, and the total evaporation flux; these quantities will serve as the objective and constraint data in the subsequent optimization. 2. Mathematical model and direct heat-conduction problem The model couples electric conduction to transient heat conduction in a thin thulium foil of prescribed geometry. Radiation and evaporation enter the thermal balance as losses from the two broad faces of the foil. 1. Rectangular wing-shaped foil We use the one-parameter foil geometry introduced in Ref. [2]. The foil has total length 2L and minimum width 2Wmin. Two symmetric rectangular wings increase the width to 2Wmax over a half- length l, with 0 ≤ l ≤ L. The limits l = 0 and l = L give uniform rectangles of widths 2Wmin and 2Wmax, respectively. The local half-width of the foil is defined as Wl(z) = {Wmax, |z| ≤ l, . Wmin, l < |z| ≤ L The corresponding two-dimensional domain is l = {(x, z) ∶ -L ≤ z ≤ L, -Wl(z) ≤ x ≤ Wl(z)} . All four geometrical parameters are fixed for a single direct calculation (see Figure 1). 2. Heat-conduction problem For a prescribed shape Wl(z) and current I, the temperature satisfies t m(T)cp(T) T = ∇ ⋅ ( (T)∇T) + e(T)|j(x, z; I, l)|2 - Qloss(T). (1) Here m, cp, , and e denote the mass density, specific heat capacity, thermal conductivity, and electrical resistivity of thulium, and j is the in-plane current-density field. The foil thickness h is small relative to its in-plane dimensions, so the temperature is taken to be uniform through the thickness. We nevertheless retain both in-plane coordinates. A one-dimensional approximation, j(z) = I/[2hWl(z)], cannot resolve current crowding at the steps in foil width; the current field is therefore obtained from a separate two-dimensional potential problem. 3. Electric-current distribution The current distribution is computed once for each geometry, before thermal integration begins. This auxiliary problem uses constant conductivity because only the geometry of the current paths is required. The normalized potential (x, z) satisfies ∇ ⋅ ( 0∇ ) = 0, 0 = 1, on the same quarter-domain mask as the heat problem. Because the equation is homogeneous, 0 does not affect the potential shape. Setting 0 = 1 simply fixes an auxiliary normalization; the field is later rescaled to I in Eq. (2). Temperature dependence enters only through e(T) in the Joule source. In the quarter-domain calculation, z = 0 is the mid-plane of the full foil and is assigned = 0.5; the right terminal at z = L is assigned = 0 for x ≤ Wmin. All remaining free boundaries are insulating. This boundary-value problem can be effectively solved by using the over-relaxation method [6, 7]. Numerically, fluxes through faces without an active neighboring node are omitted. The unscaled current-density field is j0 = - 0∇ . The total current carried by the computed quarter-domain field through the terminal section is evaluated as Wmin 0 I(1/4) = h ∫ |(j0)z (x, L)| dx. 0 The physical field used in the Joule source is then obtained by scaling j0 to the prescribed total current I: The Joule source is then ( , ; , ) = I/2 I j x z I l (1/4) 0 j0(x, z). (2) qJ(x, z, T) = e(T)|j(x, z; I, l)|2. The loss term combines radiation with the energy carried by evaporated atoms. Surface fluxes from the two broad faces are converted to a volumetric rate by dividing by the foil thickness: h Qloss(T) = 2 [ rad(T) + evap(T)] . The radiative heat flux from one surface [8] is written as env rad(T) = SB (T4 - T4 ) , where is the effective hemispherical emissivity, SB is the Stefan-Boltzmann constant, and Tenv is the effective radiation temperature of the chamber. The heat flux carried away by evaporation from one surface is estimated as N evap(T) = (T) Hsub . A Here (T) is the evaporation flux, Hsub is the molar enthalpy of sublimation, and NA is Avogadro’s constant. We use a constant reference value for Hsub. Table 1 Saturated vapor pressure and estimated evaporation flux of metallic thulium. The flux was calculated using the Hertz-Knudsen relation with the evaporation coefficient taken equal to unity. The vapor-pressure relation is recommended for 298 K to 1400 K; values below 298 K are shown only as a formal extrapolation. T (K) Psat (Pa) (atom cm-2 s-1) 100 1.88 × 10-111 3.82 × 10-93 200 2.17 × 10-50 3.12 × 10-32 298 2.22 × 10-30 2.61 × 10-12 300 4.16 × 10-30 4.87 × 10-12 400 5.30 × 10-20 5.37 × 10-2 500 5.84 × 10-14 5.30 × 104 600 6.04 × 10-10 5.00 × 108 678 1.21 × 10-7 9.42 × 1010 700 4.35 × 10-7 3.33 × 1011 778 2.25 × 10-5 1.63 × 1013 800 5.94 × 10-5 4.26 × 1013 900 2.69 × 10-3 1.82 × 1015 1000 5.61 × 10-2 3.59 × 1016 At steady state, Eq. (1) reduces to ∇ ⋅ ( (T)∇T) + e(T)|j(x, z; I, l)|2 - Qloss(T) = 0. The in-plane free boundaries are thermally insulated. Radiation and evaporation act through Qloss, not through those edge conditions. 4. Evaporation flux from the temperature field Once the stationary temperature is known, the thulium flux follows directly from the local surface temperature. For solid thulium, the saturated vapor pressure is estimated using the CRC vapor-pressure relation [9] log ( Psat ) = 13.888 - 12270 - 0.9564 log T, 298 K ≤ T ≤ 1400 K. (3) 10 Pa T 10 We calculate the equilibrium flux from the Hertz-Knudsen relation (T) = Psat(T) , √2 mTmkBT where kB is Boltzmann’s constant and mTm = 168.934 21 u is the atomic mass of thulium. We set the evaporation coefficient to unity, so (T) is an equilibrium upper estimate. Table 1 shows the strong temperature dependence. Values below 298 K are formal extrapolations of Eq. (3). For a prescribed wing length l and electric current I, the total evaporation flux from one large surface of the foil is defined as (l; I) = ∫ (Tl,I(x, z)) dx dz. l Figure 2. Radiative and evaporative surface heat fluxes of the thulium foil versus temperature. Solid line: radiation, env rad = SB(T4 - T4 ) with = 0.10; dashed line: evaporation, evap = (T) Hsub/NA, with (T) from the Hertz-Knudsen relation. Below 1000 K radiation exceeds evaporation by more than three orders of magnitude and dominates the surface energy balance This definition refers to one face of the foil. Since the numerical calculation is performed only in one quarter of the domain, the full one-side evaporation flux is reconstructed as The total from both faces is 2 . (l; I) = 4 ∫ l (1/4) (Tl,I(x, z)) dx dz. (4) The number of atoms that ultimately reach the ionization region will be lower than this equilibrium estimate. It also depends on source geometry, transport conductance, background pressure, and deposition on surrounding surfaces. 5. Relative importance of radiation and evaporation env Figure 2 compares the two surface-loss channels. Below 1000 K the radiative flux rad = SB(T4 - T4 ) exceeds the evaporative flux evap = (T) Hsub/NA by at least three orders of magnitude and is therefore the principal loss term. The approximate balance qJ ≃ SBT4 gives Tmax ∝ -1/4. In the numerical model, Tmax changes from 985 K at = 0.05 to 700 K at = 0.20, compared with 830 K at the baseline value = 0.10. This relatively weak fourth-root dependence supports the use of a single effective emissivity. Although evaporation contributes little to cooling below 1000 K, we retain it because (T) is the quantity needed in the subsequent optimization over wing length and operating current. In the regime studied here it uses for diagnostic rather than for cooling. 6. Direct problem For each prescribed pair (l, I), we compute the transient field T(x, z, t) and its stationary limit Tl,I(x, z). The prescribed wing length is aligned with the computational grid: l = kl z, kl = 0, 1, … , Nz, so that the wing edge coincides with a grid line. The current remains fixed during each run. For the computed stationary temperature field, the following output quantities are evaluated: Tmax(l, I) = max (x,z)∈ l Tl,I(x, z), and the total evaporation rate of the surface (4). Optimization over l and selection of a current subject to a maximum- temperature constraint are left to an outer procedure and are not considered here. 3. Thermophysical coefficients used in the numerical calculations The model requires m(T), (T), cp(T), and e(T) over the operating range of the foil. Earlier calculations used room-temperature constants [2]; here we construct temperature-dependent functions from experimental and reference data over approximately 60 K to 800 K. 1. Mass density We reconstructed the density from the linear-expansion data in the TPRC compilation [10], which for thulium draws on the measurements of Singh, Khanduri and Tiang between 90 K and 300 K [11]. Taking m(293 K) = 9.33 g cm-3, as the reference density, the density at temperature T was calculated as m(T) = m(293 K) , [1 + L(T)/L0]3 where L(T)/L0 is measured relative to room temperature. We assume isotropic expansion, as no crystallographic texture is specified for the polycrystalline foil. The resulting correction is small but is retained in the calculations (Table 2). Mass density of metallic thulium reconstructed from thermal expansion data [10, 11]. Table 2 T (K) L/L0 (%) m (g cm-3) 60 -0.163 9.376 78 -0.153 9.373 100 -0.140 9.369 200 -0.083 9.353 293 0.000 9.330 400 0.142 9.290 500 0.266 9.256 600 0.386 9.223 700 0.505 9.190 778 0.599 9.164 800 0.625 9.157 2. Thermal conductivity Single-crystal thulium is anisotropic because of its hexagonal close-packed structure. Ho, Powell and Liley report recommended conductivities parallel and perpendicular to the c axis, together with a polycrystalline value [12]. We use the latter for the foil. The source values, originally in W cm-1 K-1, were converted using 1 W cm-1 K-1 = 100 W m-1 K-1. No additional Voigt-Reuss-Hill averaging was applied. Because the recommended data end near room temperature, the code uses the following smooth continuation rather than a piecewise-constant extension: (T) = 0 + A eff , (5) Tp e(Teff) where 0 = -0.821 188 775 W m-1 K-1, A = 1.07573654 × 10-7, p = 0.829022226. The effective temperature is limited by Teff = min(T, Tmax), Tmax = p 0 , b(1 - p) where 0 and b are defined in the resistivity model below. The cap prevents an unphysical decrease beyond the maximum of the fitted curve. This continuation is a modelling assumption introduced to keep the nonlinear coefficient smooth. 3. Specific heat capacity Below 360 K, the heat capacity is taken from the measurements of Jennings, Hill and Spedding [13], which cover 15 K to 360 K. The original data are given as molar heat capacity Cp,m(T) in J mol-1 K-1. They were converted to mass specific heat by c (T) = Cp,m(T) , M = 168.94 g mol-1. p MTm Tm The difference between cp and cV is small relative to other model uncertainties, so Eq. (1) uses cp(T). Above 360 K, we extend the table using the trend reported in the TPRC compilation [14]. This continuation is less direct than the low-temperature measurements and is used only to cover the foil’s operating range. A constant value cp = 160 J kg-1 K-1 is retained for comparisons with the earlier constant-coefficient model. 4. Electrical resistivity Electrical resistivity directly controls the Joule source: qJ(x, z, T) = e(T)|j(x, z; I, l)|2. For 60 K to 300 K we use the polycrystalline measurements of Colvin, Legvold and Spedding [15]; these data are more appropriate for a foil of unspecified texture than single-crystal values. Above 300 K, we extrapolate from the approximately linear trend in the single-crystal measurements of Edwards and Legvold [16], which provide a(T), b(T), and c(T) up to room temperature. The single-crystal data distinguish the basal-plane and c-axis directions. To obtain a scalar polycrystalline estimate, the conductivity was averaged as (T) i(T) = 1 , i Table 3 Recommended thermal conductivity of metallic thulium from Ho, Powell and Liley [12]. Original values are given in W cm-1 K-1 and converted here to W m-1 K-1. The column kpoly was used in the numerical calculations. T (K) k∥c (W m-1 K-1) k⟂c (W m-1 K-1) kpoly (W m-1 K-1) 50 14.2 9.18 10.2 58 17.0 10.4 11.0 60 17.9 10.4 11.2 70 21.4 10.9 11.8 80 25.4 10.7 12.7 90 24.8 10.5 13.1 100 24.6 11.1 13.5 123.2 24.4 11.7 14.4 150 24.3 12.6 15.3 173.2 24.8 13.0 15.8 200 23.2 13.4 16.2 223.2 22.1 13.8 16.4 250 20.2 13.9 16.7 273.2 16.8 14.0 16.8 298.2 13.8 14.0 16.9 300 12.3 14.1 16.9 400 - - 16.9 500 - - 16.9 600 - - 16.9 700 - - 16.9 778 - - 16.9 800 - - 16.9 and then inverted: eff(T) = 2 a/b(T) + c(T) , 3 (T) eff(T) = 1 . eff The values appear in Table 5; the code represents them by the smooth linear fit e(T) = 0 + bT, with 0 = 1.662 774 21 × 10-7 Ω m, b = 1.730 384 89 × 10-9 Ω m K-1. The same differentiable fit is used in the Joule source and in Eq. (5). Table 4 Specific heat capacity of metallic thulium used in the numerical calculations. Values below 360 K are based on the measurements of Jennings, Hill and Spedding [13]; the high-temperature continuation is based on reference data [14]. T (K) Cp,m (J mol-1 K-1) cp (J g-1 K-1) cp (J kg-1 K-1) 60 24.22 0.143 143 70 24.51 0.145 145 75 24.67 0.146 146 78 24.77 0.147 147 80 24.83 0.147 147 90 25.12 0.149 149 100 25.25 0.150 150 150 25.91 0.153 153 200 26.52 0.157 157 250 26.75 0.158 158 273.15 26.88 0.159 159 298.15 27.02 0.160 160 300 27.03 0.160 160 360 27.37 0.162 162 400 - 0.160-0.165 160-165 500 - 0.165 165 600 - 0.171 171 700 - 0.177 177 778 - 0.180 180 800 - 0.182 182 4. ADI finite-difference scheme on the quarter domain Symmetry about the longitudinal and transverse axes reduces the calculation to the quarter-domain l (1/4) = {(x, z) ∶ 0 ≤ z ≤ L, 0 ≤ x ≤ Wl(z)} . Zero normal heat flux is imposed at x = 0 and z = 0. At the terminal, Ti,Nz = T0, 0 ≤ i ≤ imin, with T0 = 78 K. The other in-plane boundaries are insulated. A uniform Cartesian grid is introduced in the enclosing rectangle [0, Wmax] × [0, L]: xi = i x, i = 0, … , Nx, zk = k z, k = 0, … , Nz. The grid is chosen so that Wmin = imin x, Wmax = imax x, l = kl z. Table 5 Electrical resistivity of metallic thulium used in the numerical calculations. Values up to 300 K are based on polycrystalline data [15]; values above 300 K are model extrapolations based on single-crystal data [16]. T (K) e (µΩ cm) e (Ω m) Origin 60 26.5 2.65 × 10-7 polycrystalline data 70 28.0 2.80 × 10-7 polycrystalline data 78 29.5 2.95 × 10-7 polycrystalline data 80 30.0 3.00 × 10-7 polycrystalline data 100 34.0 3.40 × 10-7 polycrystalline data 120 37.5 3.75 × 10-7 polycrystalline data 160 44.5 4.45 × 10-7 polycrystalline data 200 52.0 5.20 × 10-7 polycrystalline data 240 59.0 5.90 × 10-7 polycrystalline data 280 65.5 6.55 × 10-7 polycrystalline data 300 69.0 6.90 × 10-7 polycrystalline data / joining point 400 86.1 8.61 × 10-7 extrapolation 500 103.2 1.032 × 10-6 extrapolation 600 120.4 1.204 × 10-6 extrapolation 700 137.6 1.376 × 10-6 extrapolation 778 151.0 1.510 × 10-6 extrapolation 800 154.7 1.547 × 10-6 extrapolation We normally set imax = Nx. Active rows contain i = 0, … , iend(k), where iend(k) = {imax, 0 ≤ k ≤ kl, , imin, kl < k ≤ Nz and the active part of a column is k = 0, … , kend(i), where kend(i) = {Nz , 0 ≤ i ≤ imin, . kl, imin < i ≤ imax Before time integration, the current field is stored as j 2 i,k = |j(xi, zk; I, l)|2. At a Picard iterate T˜, the thermal coefficients are Ci,k(T˜) = m(T˜i,k)cp(T˜i,k), i,k(T˜) = (T˜i,k), and i, qeff,i,k(T˜) = j2k e(T˜i,k) - Qloss(T˜i,k). Here Qloss(Tn ) = 2 [ SB ((Tn )4 - T4 ) + (Tn ) Hsub ] . i,k h i,k env i,k NA We discretize diffusion in conservative flux form. At an active node (i, k), the contributions through the vertical faces are F x i+1/2,k n = { i+1/2,k 0 (Ti+1,k - Ti,k) , (i + 1, k) ∈ (1/4), l , (1/4) , (i + 1, k) ∉ l and F x i-1/2,k n = { i-1/2,k 0 (Ti,k - Ti-1,k) , (i - 1, k) ∈ (1/4), l . (1/4) , (i - 1, k) ∉ l Similarly, the fluxes through the horizontal faces are F z i,k+1/2 n = { i,k+1/2 0 (Ti,k+1 - Ti,k) , (i, k + 1) ∈ (1/4), l , (1/4) , (i, k + 1) ∉ l and F z i,k-1/2 n = { i,k-1/2 0 (Ti,k - Ti,k-1) , (i, k - 1) ∈ (1/4), l . (1/4) , (i, k - 1) ∉ l When both neighboring nodes are active, face conductivities are arithmetic means: (Tn ) + (Tn ) n i+1/2,k = i,k i+1,k , 2 (Tn ) + (Tn ) n i,k+1/2 = i,k i,k+1 . 2 If the neighboring node lies outside the foil, the corresponding flux is zero; no ghost temperature is introduced. The corner (imin, kl) requires no special condition. The discrete operators are then written as Fx - Fx x L(n)T i,k = i+1/2,k i-1/2,k , x2 and Fz - Fz z L(n)T i,k = i,k+1/2 i,k-1/2 . z2 Time integration uses a Peaceman-Rachford ADI scheme [17-19], with Picard iteration at each half-step. The first half-step is implicit in x and explicit in z: Tn+1/2,(s+1) - Tn Cn+1/2,(s) i,k i,k = Ln+1/2,(s)Tn+1/2,(s+1) + LnTn + qn+1/2,(s). (6) i,k t/2 x i,k z i,k eff,i,k For each fixed k, Eq. (6) gives a tridiagonal system along the active row i = 0, … , iend(k): ax Tn+1/2,(s+1) + bx Tn+1/2,(s+1) + cx Tn+1/2,(s+1) = rx . (7) i,k For ordinary active nodes, i-1,k i,k i,k i,k i+1,k i,k a x i,k n+1/2,(s) , c = - i-1/2,k x2 x i,k n+1/2,(s) = - i+1/2,k x2 , b n+1/2,(s) n+1/2,(s) n+1/2,(s) and x i,k 2Ci,k = t 2Cn+1/2,(s) i-1/2,k + i+1/2,k + x2 , rx = i,k Tn + LnTn + qn+1/2,(s). (8) i,k t i,k z i,k eff,i,k If a neighboring node in the x direction is absent, the corresponding coefficient ax or cx is set to zero. In particular, at x 0 the symmetry condition gives ax i,k i,k c x iend(k),k = = 0. 0,k = 0, while at the right end of an active row The row k = Nz is not advanced by the x-sweep. It belongs to the current terminal for 0 ≤ i ≤ imin and is fixed as i,Nz Tn+1/2 = T0, 0 ≤ i ≤ imin. The second half-step is implicit in z. The accepted intermediate field provides the explicit x contribution, while Picard iteration updates the implicit coefficients: Tn+1,(s+1) - Tn+1/2 Cn+1,(s) i,k i,k = Ln+1,(s)Tn+1,(s+1) + Ln+1/2Tn+1/2 + qn+1,(s). (9) i,k t/2 z i,k x i,k eff,i,k For each fixed i, Eq. (9) gives a tridiagonal system along the active column k = 0, … , kend(i): az Tn+1,(s+1) + bz Tn+1,(s+1) + cz Tn+1,(s+1) = rz . (10) For ordinary active nodes, i,k i,k-1 i,k i,k i,k i,k+1 i,k n+1,(s) az = - i,k-1/2 , cz n+1,(s) = - i,k+1/2 , i,k z2 n+1,(s) i,k z2 n+1,(s) n+1,(s) b and z i,k 2Ci,k = t + i,k-1/2 + i,k+1/2 , z2 2Cn+1,(s) rz = i,k Tn+1/2 + Ln+1/2Tn+1/2 + qn+1,(s). i,k t i,k x i,k eff,i,k At z = 0, the symmetry condition gives az = 0. For columns inside the side wing, i > imin, the active segment ends at k k i,0 cz 0. For columns inside the = l, and the upper face is insulated; therefore i,kl = central part, 0 ≤ i ≤ imin, the active segment extends to k = Nz, where the Dirichlet condition is imposed by replacing the last equation with i,Nz Tn+1 = T0, 0 ≤ i ≤ imin. Equivalently, a z i,Nz = 0, b z i,Nz = 1, c z i,Nz = 0, r z i,Nz = T0. The Picard iteration in each half-step is stopped when both an absolute and a relative update criterion are satisfied, max |T(s+1) - T(s)| < nl,abs, l (i,k)∈ (1/4) | i,k i,k | max |T(s+1) - T(s)| (i,k)∈ (1/4) | i,k i,k | i,k l < nl,rel. (11) max (1, max l (i,k)∈ (1/4) |T(s+1)|) After each full step, stationarity is tested with the maximum temperature rate, | Tn+1 - Tn | Rn = max | i,k i,k | . (12) T (1/4) | t | (i,k)∈ l | | T The loop ends when Rn < rate or tn ≥ tmax. The default tmax = DBL_MAX - 10 t effectively disables the time limit, so production runs normally stop on the rate criterion. The elapsed physical time then estimates the foil’s approach time to steady operation. 5. GPU implementation of the ADI solver We implemented the ADI scheme in CUDA C for a single GPU. The row systems in Eq. (7) and the column systems in Eq. (10) are independent within their respective half-steps. We compare several ways of solving these systems and dispatching the associated kernels; all produce the same discrete solution. Figure 3 summarizes the implementation. 1. Overall organization i,k All field arrays, including the two temperature layers, squared current density j2 , and the temperature- dependent coefficients Ci,k, i,k and qeff,i,k, are stored on the device in row-major order over the enclosing rectangle, IDX(i, k) = k(Nx + 1) + i. The auxiliary electric-current problem is linear and i,k is solved once on the host before the time loop; only the resulting field j2 is copied to the device. The transient integration then runs entirely on the GPU. Host-device transfers during integration are limited to the maximum temperature change and maximum temperature, and occur only when a convergence test is evaluated. Each time step is carried out by a small set of kernels: a coefficient kernel that evaluates Ci,k, i,k and qeff,i,k from the current temperature estimate with one thread per node; the two sweep kernels that assemble and solve the tridiagonal systems of the first and second half-steps; and reduction kernels that return rates for the Picard test (11) and for the stationary-rate test (12). The two temperature layers are kept in two device buffers that are exchanged by pointer swap after each accepted step, so that no full-field copy is performed inside the time loop. 2. Mapping of the tridiagonal systems onto the device Each tridiagonal system is assigned to one CUDA block, and the threads of that block cooperate on that single system rather than each thread owning a complete system. This keeps the three diagonals and the right-hand side of a given system in a contiguous span of device memory, so that the parallel assembly and scatter phases issue coalesced accesses, and it allows the block to use syncthreads() to separate the assembly, the solve and the scatter phases. The block dimension is fixed at 256 threads. Because the quarter domain is L-shaped, the systems do not all have the same length. The rows split into long rows in the wide central part (0 ≤ k ≤ kl, length imax + 1) and short rows in the side wing (kl < k < Nz, length imin + 1); the columns split analogously into long columns in the central part (0 ≤ i ≤ imin, length Nz + 1) and short columns in the side wing (imin < i ≤ imax, length kl + 1). Each group is handled by a separate kernel launch, so that a launch processes a batch of equal-length systems and the storage for that batch is sized to the exact system length rather than to the largest system in the domain. For the grid used here (Nx = 50, Nz = 400, imax = 50, imin = 20, kl = 280) the four batches contain 281, 119, 21 and 30 systems of length 51, 21, 401 and 281, respectively. For a batch of Nsys systems of length nmax, four device arrays of size Nsys nmax hold the sub-, main- and super-diagonals and the right-hand side; system s occupies the contiguous span starting at s nmax, and its solution overwrites the right-hand side in place. 3. Assembly, solve and scatter within a block Within a block the work proceeds in three barrier-separated phases. In the assembly phase the threads stride over the active nodes of the system, count = tid, tid + blockDim, …, and each thread writes the coefficients a, b, c and the right-hand side r of its nodes according to Eqs. (7)-(8) for the x-sweep and the analogous relations for the z-sweep. The face conductivities are evaluated by arithmetic averaging of the neighbouring nodal values, and a neighbour outside the quarter domain contributes a zero flux, so that the corresponding off-diagonal coefficient vanishes. The symmetry conditions at x = 0 and z = 0 enter as a vanishing lower diagonal at count = 0, and the end of each active segment as a vanishing upper diagonal at count = n - 1; the Dirichlet terminal node (i, Nz) with i ≤ imin is assembled directly as az = 0, bz = 1, cz = 0, rz = T0, i.e. the last equation of the affected columns is replaced rather than eliminated. After a barrier the solve phase reduces each system to its solution, and after a second barrier the scatter phase writes the computed line back into the two-dimensional temperature array, again with the threads striding over the nodes of the system. The same kernel structure is used for both half-steps; only the index running along the system (the row index i for the x-sweep, the column index k for the z-sweep) and the direction of the implicit and explicit operators differ. The assembly and scatter phases are identical across all variants considered below; the variants differ only in the algorithm that implements the solve phase. 4. Tridiagonal-solve algorithms We compared four solvers that differ in how elimination is distributed across the block [20]. Their solutions agree to round-off. 1. Serial Thomas sweep (baseline). Thread 0 executes the classical Thomas algorithm [21]: forward elimination followed by backward substitution. A pivot smaller than 10-300 marks the solution as NaN. The guard was never triggered because the transient term 2Ci,k/ t makes these systems strongly diagonally dominant. This implementation matches the serial CPU sweep bit-for-bit, but the device is under-utilized during this phase because the other block threads remain idle. 2. Two-sided Thomas algorithm. This variant applies Thomas elimination simultaneously from both ends of the system [22, 23]. The two elimination fronts advance towards a central index, after which substitution proceeds outwards. Two threads per block carry out the sweeps, halving the length of the inherently sequential dependency chain. 3. Austin-Berndt-Moulton partition method. Production runs use the memory-efficient partitioning scheme of Austin, Berndt, and Moulton (ABM) [24]. A system of N rows is divided into NP contiguous partitions, each of length Nk ≈ N/NP. Each thread eliminates the interior unknowns of one partition and expresses them in terms of two interface unknowns. These interface unknowns form a reduced tridiagonal system of order 2NP. After this system is solved by thread 0 using the serial Thomas algorithm, the threads recover the interior unknowns of their respective partitions in parallel. The ABM method stores only the boundary-coupling coefficients rather than the full fill-in vectors, requiring approximately 6 kB of shared memory in our kernels. The partition method is applied to systems with N > 100. Smaller systems containing 21 and 51 unknowns, are solved using two-sided elimination. 4. Vendor batched solver (cuSPARSE). For an external comparison, the batched solve was delegated to the NVIDIA cuSPARSE library [25] routine cusparseDgtsv2StridedBatch, which solves a batch of equal-length tridiagonal systems stored in the strided layout already used by the assembly phase. This variant reuses the same assembly and scatter kernels and replaces only the solve, so it measures the cost of a production-quality library solver on the identical data layout. 5. Kernel launch strategy We also compared two launch strategies. In plain mode every kernel of a step - the coefficient kernel, the four assembly/solve batches of each sweep, and the reductions - is launched individually through the ordinary runtime API. Because a single step of the present problem issues of the order of a few tens of very short kernels (each system is at most 401 unknowns long, so a kernel runs only for microseconds), the fixed CPU-side cost of launching each kernel becomes comparable to the kernel’s own execution time, and over the ∼ 8 × 104 steps of a full run this launch overhead accumulates into a noticeable fraction of the wall time. Read parameters; allocate device buffers Solve Φ (SOR) and precompute j2 [CPU] i,k init kernel: [GPU] set T 0 = Tenv, t = 0 Time loop, step n: clip ∆t = tmax - t if overshoot compute coefficients kernel: C, λ, qeff at all nodes First half-step (implicit in x) x Picard s = 1, 2, . . . no Picard converged? yes adi x lines 400 CUDA blocks adi z lines 51 CUDA blocks one CUDA block = one tridiagonal system assembly all threads, coalesced Second half-step (implicit in z) z Picard s = 1, 2, . . . no Picard converged? syncthreads() tridiagonal solve baseline: serial (thread 0); production: ABM partition (all threads) yes swap T n ↔ T n+1; t ← t + ∆t; reduce RT = max |∆T |/∆t syncthreads() scatter all threads → T ∗ no RT < εrate yes or t ≥ tmax? cudaMemcpy T → host; write T (x, z); STOP Figure 3. Structure of the GPU implementation of the ADI solver. The host solves the electric problem once; the transient loop runs entirely on the device. Each ADI half-step launches one kernel per batch of equal-length tridiagonal systems (the x-sweep is shown exploded on the right: 281 long rows of length 51 and 119 short rows of length 21; the z-sweep is organized analogously into 21 long and 30 short columns). Every tridiagonal system is mapped to a single CUDA block, whose threads cooperate in three barrier-separated phases - parallel assembly, the tridiagonal solve, and parallel scatter of the solution. The solve phase is interchangeable: the baseline runs a serial Thomas sweep in thread 0, whereas the production configuration distributes it across the threads of the block through the ABM partition method (Section 5) 1. Electric potential (x, z) for the full left-to-right problem. The potential decreases monotonically from the left contact, where = 1, to the right contact, where = 0. The central value ≃ 0.5 corresponds to the symmetry plane used in the quarter-domain calculation 2. Distribution of the electric-current-density magnitude |j(x, z)| in the full plate domain. The current density increases in the narrow end regions because of the reduced cross-sectional width and exhibits local concentrations near the width-transition corners Figure 4. Electric properties of the full plate geometry: (a) electric potential and (b) magnitude of the electric current density The second mode uses CUDA task graphs [26]. The kernel-call flow within one time step is represented as a task graph and submitted as a single GPU operation. We construct the graph implicitly by capturing the stream workflow with cudaStreamBeginCapture and cudaStreamEndCapture, instantiate it with cudaGraphInstantiate, and execute it with cudaGraphLaunch. This approach is particularly effective for short kernels, whose execution time can be comparable to their launch latency, because it reduces the gaps between successive kernel launches [27]. Since t, the grid, and the batch structure remain fixed, the same graph is used throughout a run. A shortened final step, if needed to reach tmax exactly, is launched separately. 6. Time loop and stopping The Picard iterations of each half-step reuse the coefficient and sweep kernels: the coefficients are re-evaluated from the latest temperature estimate, the tridiagonal systems are re-assembled and re- solved, and the maximum update is reduced on the device. To limit host-device synchronization, the update norm is copied back and tested only on every third iteration and on the last allowed iteration. After a full step the maximum temperature change is reduced on the device and divided by t to form the stationary rate Rn of Eq. (12). The time loop terminates when Rn < rate or when the accumulated T T physical time reaches the prescribed bound tmax; when the next step overshoot tmax, the step size is reduced to t = tmax - tn for that final step in order to land exactly on tmax. 6. Numerical results We solved the coupled problem for a total current I = 2.5 A, with the cold contacts held at T0 = 78 K, the surrounding medium at Tenv = 293 K, and a surface emissivity = 0.10. The grid used Nx = 50, Nz = 400 (grid_scale = 1), with the wing transition at kl = 280 and imin = 20. The time step was t = 1 × 10-3 s, and each Picard update used at most 20 iterations per half-step with tolerances nl,abs = 10-6 and nl,rel = 10-10. The potential decreases monotonically between the contacts (Figure 4a). As expected, the current density is largest in the narrow end regions and shows additional crowding near the re-entrant corners (Figure 4b). This nonuniform |j|2 field is used in the Joule source. Starting from T0 = 78 K, the calculation reached the criterion RT = maxi,k |Tn+1 - Tn |/ t < rate i,k i,k with rate = 1 × 10-2 K s-1 after 77 397 steps, or 77.397 s of simulated time. The maximum temperature was 829.803 K. Figures 5a and 5b show the stationary field. The cold contacts keep the narrow ends near 78 K, while the maximum lies in the broad central region. Changes in the field near the width transitions reflect the redistribution of current there. 1. Two-dimensional temperature map of the stationary regime in the full symmetric thulium plate domain. The narrow end regions are maintained close to the prescribed contact temperature T0 = 78 K, whereas the central wide part reaches the highest temperature 2. Three-dimensional representation of the stationary temperature field T(x, z) in the full symmetric thulium plate domain. Local distortions near the width-transition zones reflect current redistribution and the associated spatial variation of Joule heating Figure 5. Stationary temperature distribution in the full symmetric thulium plate domain: a two-dimensional temperature map and b three-dimensional representation of the temperature field. The temperature decreases toward the narrow cold-contact regions, while the maximum occurs in the central part of the wider section All tridiagonal solvers produced the same converged field and maximum temperature, allowing their execution times to be compared directly. 7. Performance of the tridiagonal-solver strategies We benchmarked the batched solvers on one NVIDIA V100 against a serial CPU implementation. Every variant used Nx = 50, Nz = 400, kl = 280, and t = 1 × 10-3 s, stopped at step 77 397, and produced a bit-for-bit identical Tmax = 829.803 K. Wall-clock differences therefore arise from the solver and launch strategy rather than from different stopping histories. The tested configurations comprise the one-thread Thomas baseline, the two-sided Thomas algorithm, the ABM partition solver, and cusparseDgtsv2StridedBatch. Each in-house solver was measured with ordinary kernel launches and with CUDA task graphs. We also tested device intrinsics in the partition sweep. Each configuration was run five times; Table 6 reports the mean and standard deviation; Figure 6 and 7 show the same data graphically. The serial CPU run took 876 s. ABM partitioning achieved 21-24× speedup (Figure 7), roughly twice that of the serial and two-sided Thomas algorithms (11-12×), because the elimination work is shared across the block. CUDA task graphs reduced runtime by about ten percent over the 77 397 time steps. Explicit use of FMA intrinsics produced no speedup because nvcc enables contraction of eligible multiply-add operations by default (-fmad=true). The fastest configuration, ABM partitioning with a CUDA task graph, completed in 36.3 s (24.1×; Figure 6), about ten percent faster than cuSPARSE at 40.2 s (21.8×). The GPU computations were performed on an NVIDIA Tesla V100-SXM2 (16 GB) accelerator of the “Govorun” supercomputer at the HybriLIT heterogeneous computing platform [28], using code built with the CUDA 11.2 toolkit and GCC 9.5.0 as the host compiler (-arch=sm_70 -O3). For comparison, the CPU version was run on the same platform, on Intel Xeon Platinum 8268 processors (24 cores, 2.9 GHz, 192 GB DDR4) and built with the same compiler. 8. Discussion The calculated temperature field reflects the combined effects of the nonuniform current distribution and heat removal through the cold contacts. The maximum temperature occurs in the broad central Figure 6. Mean wall-clock time of the tridiagonal-solver configurations on an NVIDIA V100 (logarithmic scale; lower is better; error bars show the standard deviation over five runs). The dashed line shows the cuSPARSE timing. Colours distinguish the tridiagonal algorithms. All configurations converge to the same Tmax = 829.803 K at step 77 397 Figure 7. Speedup of the GPU configurations relative to the single-thread CPU run (linear scale; higher is better). The dashed line shows the cuSPARSE result. The ABM partition configurations reach 21-24× speedup; the two configurations that use CUDA task graphs outperform the vendor routine Table 6 Wall-clock time (mean ± standard deviation over five runs) and speedup relative to the single-thread CPU run for the tridiagonal-solver strategies on an NVIDIA V100. All configurations converge at step 77 397 to Tmax = 829.803 K. Configuration Wall time (s) Speedup ABM partition, + task graph 36.3 ± 1.2 24.1× ABM partition, + intrinsics + task graph 36.9 ± 0.8 23.7× cuSPARSE 40.2 ± 3.3 21.8× ABM partition, plain 40.7 ± 1.0 21.5× ABM partition, + intrinsics 41.5 ± 1.7 21.1× Two-sided Thomas, + task graph 72.0 ± 0.6 12.2× Serial Thomas, + task graph 72.4 ± 1.1 12.1× Serial Thomas, plain 81.1 ± 2.8 10.8× Two-sided Thomas, plain 81.3 ± 2.1 10.8× CPU serial 876.0 1.0× region, whereas the narrow end regions remain close to the prescribed contact temperature. The local changes in the temperature field near the width transitions are associated with current crowding at the re-entrant corners. This behaviour confirms the need for a two-dimensional electric model, since the one-dimensional approximation cannot resolve these local variations in the current density. The present results should be considered in the context of the model assumptions. The electric- current distribution is calculated using constant conductivity and is subsequently scaled to the prescribed total current, while the temperature dependence of the electrical resistivity is included in the Joule heat source. In addition, a constant effective emissivity is used, and some thermophysical properties are extrapolated beyond the available experimental ranges. The numerical results therefore characterize the selected foil geometry and operating conditions rather than providing a universal optimum. The ADI discretization produces batches of independent tridiagonal systems. Mapping each system to a CUDA block makes it possible to change the solution algorithm without modifying the assembly procedure or data layout. For the grid considered here, ABM partitioning approximately doubled the performance relative to the serial and two-sided Thomas algorithms. CUDA task graphs provided a further improvement of about ten percent, whereas device intrinsics changed the execution time by less than the run-to-run variation. The fastest configuration, ABM partitioning combined with a CUDA task graph, completed the calculation in 36.3 s and achieved a 24.1× speedup over the single-thread CPU implementation. It was also approximately ten percent faster than cuSPARSE. Compared with the earlier formulation [2], which ranked candidate geometries according to the area within a prescribed temperature interval, the present model provides the evaporation flux as a more direct measure of ion-source performance. The computational efficiency of the GPU implementation makes the direct solver suitable for repeated calculations in an optimization procedure in which the foil geometry and operating current are varied. Future development of the model will couple the electric field to the evolving temperature-dependent conductivity. 9. Conclusions A two-dimensional electric-thermal model of a wing-shaped thulium foil was developed. The model accounts for temperature-dependent material properties, nonuniform Joule heating, radiation, and evaporation. The numerical results show that the proposed approach can determine the stationary temperature field and provide the quantities needed to evaluate the foil performance. All tested tridiagonal solvers produced the same converged solution, while the GPU implementation substantially reduced the computation time. The resulting solver can therefore be used for repeated calculations when optimizing the foil geometry and operating current.

Об авторах

А. С. Айриян

Объединённый институт ядерных исследований; Национальная научная лаборатория им. А. Алиханяна

Email: ayriyan@jinr.ru
ORCID iD: 0000-0002-5464-4392

Candidate of Physical and Mathematical Sciences, Head of the Nonlinear Systems Modeling Methods Sector, Scientific Department of Computational Physics, Meshcheryakov Laboratory of Information Technologies (MLIT), Joint Institute for Nuclear Research (JINR); Research Scientist at the Division of Computational Physics and Information Technologies, A. Alikhanyan National Science Laboratory (Yerevan)

ул. Жолио-Кюки, д. 6, г. Дубна, 141980, Российская Федерация; ул. Братьев Алиханян, д. 2, г. Ереван, 0036, Армения

Э. А. Айрян

Объединённый институт ядерных исследований

Email: ayrjan@jinr.ru
ORCID iD: 0000-0001-6478-4041

Candidate of Physical and Mathematical Sciences, Assistant to the Director of the Meshcheryakov Laboratory of Information Technologies (MLIT), Joint Institute for Nuclear Research (JINR)

ул. Жолио-Кюки, д. 6, г. Дубна, 141980, Российская Федерация

Е. Е. Донец

Объединённый институт ядерных исследований

Email: edonets@jinr.ru
ул. Жолио-Кюки, д. 6, г. Дубна, 141980, Российская Федерация

М. П. Венева

Центр вычислительных наук RIKEN

Автор, ответственный за переписку.
Email: milena.p.veneva@gmail.com
ORCID iD: 0000-0002-6421-4716

Student Trainee at the RIKEN Center for Computational Science

Минатодзима-минами-мати 7-1-26, район Тюо, г. Кобе, префектура Хёго, 650-0047, Япония

Список литературы

Дополнительные файлы