Assessment of absorbed power density and temperature rise for nonplanar body model under electromagnetic exposure above 6 GHz
AAssessment of absorbed power density and temperature rise for nonplanar body model under electromagnetic exposure above 6 GHz
Yinliang Diao , Essam A. Rashed , and Akimasa Hirata Department of Electrical and Mechanical Engineering, Nagoya Institute of Technology, Nagoya 466-8555, Japan College of Electronic Engineering, South China Agricultural University, Guangzhou 510642, China Department of Mathematics, Faculty of Science, Suez Canal University, Ismailia 41522, Egypt Center of Biomedical Physics and Information Technology, Nagoya Institute of Technology, Nagoya 466-8555, Japan
Abstract
The averaged absorbed power density (APD) and temperature rise in body models with nonplanar surfaces were computed for electromagnetic exposure above 6 GHz. Different calculation schemes for the averaged APD were investigated. Additionally, a novel compensation method for correcting the heat convection rate on the air/skin interface in voxel human models was proposed and validated. The compensation method can be easily incorporated into bioheat calculations and does not require information regarding the normal direction of the boundary voxels, in contrast to a previously proposed method. The APD and temperature rise were evaluated using models of a two-dimensional cylinder and a three-dimensional partial forearm. The heating factor, which was defined as the ratio of the temperature rise to the APD, was calculated using different APD averaging schemes. Our computational results revealed different frequency and curvature dependences. For body models with curvature radii of >30 mm and at frequencies of >20 GHz, the differences in the heating factors among the APD schemes were small.
Keywords: numerical dosimetry, FDTD, conformal thermal model, staircasing error
1. Introduction
With the advancement of technology, there has been a trend of adopting higher frequencies to achieve larger bandwidths and hence higher data transfer rates. Additionally, fifth-generation mobile communications allow the use of frequency bands above those employed for current wireless communications (typically <6 GHz) (Lee et al et al et al et al et al et al et al (2016) evaluated the specific absorption rate (SAR) and temperature rise in human head models up to 30 GHz. Strong correlations between the peak averaged SAR and temperature rise were observed at frequencies up to 3β4 GHz, and weaker correlations were observed at higher frequencies. Dosimetric studies for frequencies above 6 GHz have been conducted for both near- and far-field exposure scenarios (Alekseev et al et al et al et al et al et al et al et al et al et al et al (2020) revealed that the normal component of the IPD correlates well with the skin-temperature rise, regardless of the incident angle. Several studies have been performed on the correlation between the power-density averaging area and the skin surface temperature (Hashimoto et al et al et al et al (2018) reported that the averaged APD over 4 cm exhibited a strong correlation with the skin-temperature rise in the frequency range of 30β300 GHz and was reasonable and conservative at frequencies as low as 10 GHz. Kageyama et al (2019) developed an exposure system and measured the temperature rise in forearm skin exposed to a focused beam generated by a lens antenna at 28 GHz. In the foregoing studies, cubic models with planar surfaces were generally adopted for the assessment of millimeter-wave exposure, in contrast to dosimetric studies or frequencies below ~10 GHz (Hirata et al et al et al et al et al et al et al et al et al (2007) proposed a conformal scheme for correcting the heat flux. However, this method requires accurate knowledge of the normal directions for each boundary voxel. A simple correction method was adopted by Laakso (2009); the heat convection rate for the skin voxel π» was corrected as π»/βπ , where π represents the number of neighboring air voxels. However, the total surface area of a voxel model was still slightly overestimated (Laakso and Hirata 2011). The objective of this study was to reliably assess the averaged APD and temperature rise for body models with nonplanar surfaces at frequencies above 6 GHz. Different calculation schemes for the averaged APD were proposed and evaluated. We proposed a local compensation method for the heat convection at the model boundary in bioheat calculations. The heating factors (ratios of the surface temperature rise to the APD) for the nonplanar body model were evaluated to investigate the discrepancies caused by the different curvature radii and frequencies. . Models and numerical methods Two models were adopted for evaluation of the APD and temperature rise. The first was a 2D cylindrical multilayer model, as shown in Fig. 1 a). This model comprised three types of tissues: skin (thickness of 1.4 mm), fat (4 mm), and muscle. Different outer radii of the cylindrical models (i.e., 20, 30, 40, and 50 mm) were considered. The other model was a three-dimensional (3D) partial forearm model extracted from the XCAT phantom (Segars et al d from the forearm model. For EM calculations, the FDTD method (Taflove and Hagness 2005) was used. A 15-layer convolutional perfectly matched layer (Roden and Gedney 2000) was adopted to truncate the simulation domain. For 2D analysis, both transverse-electric (TE) and transverse-magnetic (TM) plane incident waves were considered as radiation sources. For 3D analysis, a 4 Γ 1 half-wave dipole antenna array was used. The antenna array was located in the front of the forearm, as shown in Fig. 1 c). Different distances ( π ) between the antenna and the forearm, i.e., 5, 10, 15, 20, 30, and 40 mm, were considered. The four dipole elements were fed by delta-gap voltage sources with the same amplitudes and phases. The local SAR was calculated as the average of the electric field components, which were defined on the edges of the target voxel. The simulated frequencies ranged from 6 to 60 GHz for the 2D analyses. The dielectric properties were obtained from the work of Gabriel, Lau, and Gabriel (1996) and are presented in Table I. For he 3D analysis, the working frequency was set as 28 GHz. The forearm model comprised seven tissues, and their dielectric properties at 28 GHz are presented in Table II. TABLE I. Dielectric Properties of the Tissues of the Multilayer Model Tissue Conductivity π [S/m] Relative permittivity π π Freq [GHz]
6 10 20 30 45 60 6 10 20 30 45 60
Skin
Fat
Muscle
TABLE II. Dielectric Properties of the Tissues of the Forearm Model at 28 GHz Tissue Conductivity π [S/m] Relative permittivity π π Skin 25.8 16.6 Fat 1.70 3.70 Muscle 33.6 24.4 Bone (Cortical) 4.94 5.17 Bone (Cancellous) 8.87 7.51 Blood 37.0 23.9 Tendon 23.6 13.9
There are two definitions of the spatial-average APD in ICNIRP-2020 (ICNIRP 2020). The first is as follows: π ab = β¬ ππ₯ππ¦ β« π(π₯, π¦, π§) β SAR(π₯, π¦, π§)ππ§/π΄ π§ max , (1) where π§ = 0 corresponds to the body surface, π§ max encloses most of the power deposition in the body, and π΄ = 4 cm represents the averaging area. The other definition is as follows: π ab = β¬ Re[πΊ] β ππ/π΄ π΄ = β¬ Re[π¬ Γ π― β ] β ππ/π΄ π΄ , (2) where ππ represents the integral variable vector normal to the body surface. In IEEE C95.1 (Bailey et al at frequencies between 6 and 300 GHz. This veraging area generally provides good estimations of the surface temperature rise (Funahashi et al et al The widely used bioheat equation (Pennes 1948) was employed for calculation of the temperature rise inside the human body model. Under the assumption of the steady-state condition, the following equation was used for the temperature rise Ξπ : β β (πβΞπ) β π΅Ξπ + π π£ = 0, (3) where π [W/(mβΒ°C)] represents the thermal conductivity; π΅ [W/(m βΒ°C)] is a coefficient related to the blood perfusion rate; and π π£ [W/m ] represents the power-loss density, which is related to the SAR by π π£ = SAR β π , where π [kg/m ] represents the tissue mass density. The temperature rise Ξπ is defined at the center of each voxel. The Neumann boundary condition in (4) was applied to the boundary of the model: π πΞπππ + π»Ξπ = 0, (4) here π» [W/(m βΒ°C)] represents the heat convection rate from the skin to the surrounding air. In this study, π» was set as 8 W/(m βΒ°C) (Hirata et al For the implementation of the Neumann boundary condition, the heat flux through the boundary of the voxel via convection must equal the heat reaching that voxel from its neighboring voxels via conduction (Bernardi et al
π» β π , where π represents the effective area of the air/skin interface. For a voxel belonging to a planar surface parallel to the voxel surfaces, π = Ξ , whereas in stepped boundaries, π > Ξ . To reduce the uncertainties in bioheat calculations, an accurate estimation of the surface area is required. Mullikin and Verbeek (1993) proposed a simple method for estimating surface area of 3D binary objects. This method assigns surface-area weights to each boundary voxel. The boundary voxels can be classified into five types: π π , π = 1, 2, β¦ ,5 , where π represents the number of voxel faces exposed to the background. Fig. 3 a) shows a local boundary region containing different types of boundary voxels. The total surface area of the model can be estimated as follows: π΄ = (β π€ π β π π5π=1 ) β Ξ , (5) where π π represents the number of voxels of type π π , and π€ π represents the weight for voxel type π π . The following optimized π€ π values were reported: π€ = 0.8940 , π€ = 1.3409 , and π€ = 1.5879 (Mullikin and Verbeek 1993). If the spatial resolution is sufficiently high, there is no obvious deviation of the boundary region from a plane. In such regions, only three types of voxels ( π ) exist. Voxels of type π can occur on sharply curved surfaces. For type π , the weights π€ = 2 , and π€ = 8/3 were presented by Mullikin and Verbeek (1993). Nonetheless, for a sphere model, there are significantly less π voxels than π voxels; therefore, the weights of the π voxels are insignificant. A simple and straightforward compensation method involves adopting a factor π€ π for voxel type π π . Consequently, the heat flux through the boundary voxel becomes π»π€ π Ξ . Thus, the total heat flux through the body surface can be accurately estimated. However, particular attention should be paid to local stepped regions, as shown in Figs. 3 b)βd). In these locally planar regions, only single type of boundary voxel exists. With the foregoing compensation method, the local effective surface area would be underestimated as π€ π < βπ . To resolve this issue, we propose a local correction method using a 3 Γ 3 Γ 3 moving cube, which is centered at the target boundary voxel (outlined in red in Fig. 3). Assuming that π πβ² represents the otal number of voxels of type π π within the moving cube, if there is only one type of boundary voxel within the cube (i.e., π πβ² = π π΄β² , π π΄β² is the total number of boundary voxels), the heat transfer through the central voxel is compensated as π»βπΞ . Thus, the heat convection rate in stepped regions, as shown in Figs. 3 b)βd), is corrected. (a) (b) (c) (d) Fig. 3. Air/skin interface for voxel models. Black voxels indicate boundary voxels; yellow voxels indicate internal tissues. The boundary voxels are labeled according to the voxel type; the center voxels are outlined in red. The region in (a) contains different types of boundary voxels. The regions in (b) to (d) contain one type of boundary voxel. The effective surface areas for the central voxels are (b) Ξ , (c) β2Ξ , and (d) β3Ξ .
3. Validation of compensation method for bioheat calculation
To validate the proposed compensation method for the heat flux at the air/skin interface, we considered a homogeneous sphere model with radius of 60 mm. The heat conduction rate was set as 1.0 W/(mβΒ°C), and a 200-W/m heat source was uniformly distributed inside the sphere. The heat convection rate between the skin and air was set as 8 W/(m βΒ°C). The analytical solution to this problem is Ξπ(π) = β200π /6 + 0.62 Β°C, where π represents the radial distance. The calculation results of the conformal and proposed compensation methods are presented in Table III and Fig. 4. Implementation of the conformal method requires knowledge of the normal directions for each boundary voxel, which were obtained using a 3D Sobel operator in this study. As indicated by Table III, the total surface area was 2.2% larger than the real surface area for the conformal method. This resulted in an underestimation of the temperature rise, with a mean value of 0.4886 Β°C for the sphere surface, compared with the exact solution of 0.5 Β°C. The Ξπ calculated using the proposed method was almost identical to the exact solution, as shown in Table III and Fig. 5. Moreover, compared with the conformal method, the proposed method provided a smaller standard deviation of Ξπ for the boundary voxels. ABLE III. Comparison of Thermal Calculation Results Obtained using Different Compensation Methods. Surface area [mm ] Max Ξπ [Β°C] Mean surface Ξπ [Β°C] SD of surface Ξπ [Β°C] Exact solution 4.523 Γ 10 (+2.2%) 0.6069 (β2.11%) 0.4886 (β2.3%) 0.001178 Proposed method 4.521 Γ 10 (β0.04%) 0.6179 (β0.34%) 0.4997 (β0.06%) 0.001054 (a) (b) Fig. 4 Distributions of the temperature rise on the surface of a homogeneous sphere. The boundary condition was implemented using (a) conformal and (b) proposed methods. Fig. 5 Temperature-rise distributions along the sphere radius for different compensation methods. n ellipsoidal model, as shown in Fig. 6 a), was also adopted for validation. The three principal semi-axes of the ellipsoid were set as π = π =
50, and π =
100 mm. The thermal parameters were identical to those used for the spherical model. Spatial resolutions of 1, 0.5, and 0.25 mm were employed. The ellipsoid was tilted around the y-axis in steps of 15 Β° . The calculated maximum temperature rises inside the ellipsoids for different spatial resolutions and tilt angles are presented in Fig. 6 b). As shown, the maximum temperature rises inside the models and on the surfaces were almost independent of the tilt angle, and the relative standard deviations caused by the tilt angle were smaller than ~0.3%. The discrepancies of the maximum temperature rises inside the models and on the model surfaces from those obtained via the MATLAB Partial Differential Equation Toolbox TM (version R2020a) (MathWorks 2020) using the finite-element method (FEM) were within ~0.8% and ~0.5%, respectively. (a) (b) Fig. 6 Validation of the proposed compensation method using prolate ellipsoidal models. (a) Ellipsoidal model (original position), with an arrow indicating the rotation direction; (b) maximum temperature rises inside the model and on the model surface.
4. APD and temperature rise for nonplanar models
The calculated SAR and temperature rise in the 2D multilayer models for TE and TM incident waves are presented in Figs. 7 a) and b), respectively. As shown, the EM power depositions were distributed superficially above 10 GHz. The peak local SARs for the TM cases were lower than those for the TE cases. This is because for the TE cases, the electric fields were parallel to the axis of the infinitely long cylindrical model. The maximum temperature rises were generally slightly larger for the TM cases than for the TE cases. As shown in Fig. 7 b), the EM wave penetrated into the cylindrical model from different incident angles, resulting in broader SAR and temperature distributions compared with the TE case. This phenomenon was attributed to Brewsterβs effect (Li et al (a) (b) Fig. 7. Distributions of the SAR and temperature rise on the 2D multilayer models with a 20-mm radius for (a) TE and (b) TM incident waves. In Fig. 8, the orange regions indicate the integration volumes for different APD calculation schemes at various frequencies. For the cylindrical model with a 20-mm radius, increasing the depth to enclose all the power transmitted through L1 is impossible at 6 and 10 GHz. Therefore, the lengths of L2 and L3 are set as 20 mm at 6 and 10 GHz for the 20-mm-radius model. This setting is in accordance with the definition of the 10-g spatial-average SAR, where the side length of the equivalent integration cube is approximately 20 mm (Hirata et al π as π»πΉ π . At 6 GHz, π»πΉ and π»πΉ were ~20% smaller than those at higher frequencies; while π»πΉ and π»πΉ were ~13% and ~5% smaller than those at higher frequencies for TE and TM incident waves, respectively. For schemes 1 and 2, the upper bounds of the integration volume were tangent planes to the cylinders, whereas for schemes 3 and 4, the upper bounds were bent along the surface. Therefore, the integration volumes for the latter two schemes were smaller than those for schemes 1 and 2, esulting in larger heating factors. The differences between π»πΉ and π»πΉ and between π»πΉ and π»πΉ decreased with the increasing curvature radius. This is mainly attributed to the reduction in the deviation of the model surface from a plane. π»πΉ and π»πΉ decreased steadily with an increase in the radii, whereas π»πΉ and π»πΉ were relatively stable for all the radii. In general, among the schemes examined, π»πΉ seems to be most stable across all radii and frequencies. For a cylindrical radius of β₯
30 mm, the four schemes provided comparable heating-factor results, with relative standard deviations of < ~5% for frequencies of β₯
20 GHz and < ~10% for frequencies of β₯ A partial human forearm model was adopted for 3D analysis. The temperature rises were calculated using the heat convection compensation method described in Section II. The calculated local SAR and Ξπ distributions in the partial forearm model are shown in Figs. 11 a) and b), respectively. The differences in the heating factors among the different averaging schemes were marginal and smaller than those for the 2D cylindrical models. In the 2D cases, a plane wave was used as the radiation source. Fig. 12 shows the heating factors for different distances between the antenna and the forearm. The largest heating factors were observed for π = π β€
15 mm) and were mainly attributed to the oscillating peak power density in this field region (Colombi et al π =
40 mm, the heating factors were approximately 0.22 Β°Cβm /W for schemes 1 and 2 and 0.23 Β°Cβm /W for schemes 3 and 4. These values were slightly larger than those for 2D multilayer models with 30-mm radii (the curvature radius of the forearm was approximately 30β40 mm). This difference was attributed to the field non-uniformity and the thicker fat tissue of the forearm model. As reported by lekseev et al (2008), a thicker fat layer tends to produce a slightly larger surface-temperature rise in the skin. (a) (b) Fig. 11 Distributions of the (a) local SAR and (b) temperature rise in the partial forearm model. The antennaβforearm distance is 10 mm. The antenna accepted power is normalized to 20 dBm. The dark red box in (a) indicates the integration volume for APD averaging scheme 1. Fig. 12 Heating factors for a 4 Γ 1 dipole antenna array with different antennaβforearm distances.
5. Discussion and concluding remarks
In bioheat calculations, the boundary conditions significantly influence the calculation results. The discretization of the model surface with cubes results in a larger model surface area; hence, the temperature rise is underestimated. This issue is particularly crucial for the assessment of millimeter-wave exposure, owing to the extremely shallow penetration depth. In the study of Laakso (2009), the heat convection rate for the skin voxel, i.e., π» , was corrected as π»/βπ , where π represents the number of exposed faces of the voxel. However, the total surface area of the voxel model was lightly overestimated (Laakso and Hirata 2011). This problem can be solved by introducing a compensation factor π , which is defined as the ratio of the real body surface area (if known) to the calculated one. The compensated heat transfer from one boundary voxel then becomes π»βππΞ , π < 1 . This method is applicable to whole-body exposure scenarios, as the total heat flux through the body surface is corrected. For localized exposure scenarios, however, compensation of the whole-body surface area may still lead to underestimation of the heat flux though the local surface parallel to the grid axes (as shown in Figs. 3 b)βd)), because π < 1 . In this paper, a new conformal compensation method was proposed for reducing the uncertainties in FD bioheat calculations for voxel models with nonplanar surfaces. This method does not require knowledge of the normal directions of the boundary voxels and can be easily incorporated into FD iterative programs with a minimal computational load. Additionally, we confirmed that the proposed method provides an accurate estimation of the temperature rise. The variation in the maximum temperature rise was observed to be smaller than ~1% for different spatial resolutions and different model rotation angles (shown in Fig. 6). Because no clear definition of the calculation scheme for the APD is prescribed in the existing standards, four schemes were proposed, particularly for models with nonplanar surfaces. The relationships between the maximum temperature rise and the averaged APD were determined using different schemes. Compared with previous studies using planar models, comparable heating factors for the APD were observed for curvature radii larger than ~30 mm and frequencies above ~20 GHz. Above ~20 GHz, all four schemes provided stable heating factors, and the heating factors decreased below ~10 GHz. As reported by Hirata, Funahashi, and Kodera (2019), even though there may be optimized physical quantities for estimating the surface-temperature rise below 10 GHz, the APD is a good surrogate, with deviations smaller than ~20%. In general, if the curvature radius is larger than ~30 mm, the uncertainties in the heating factors for different APD schemes are small, with relative standard deviations smaller than ~5% for frequencies above ~20 GHz and within ~10% for frequencies above 6 GHz. In the 3D cases, the relative standard deviations of the heating factors are within ~3.5% for all the APD schemes. In previous studies, slightly larger heating factors were observed for TM-like exposure using 2D infinite-slab models (Samaras and Kuster 2019, Li et al et al Β° C/W β m at frequencies of β₯
10 GHz. For the 3D forearm model, when the antennaβforearm istance was shorter, the EM power was more narrowly distributed, and larger heating factors were observed. In such cases, the coupling effects between the antenna and the body are predominant under exposure to a reactive near field (Funahashi et al et al et al (2017), who reported that the heating factors for beams with diameters larger than the side length of the averaging area are comparable to those for incident plane waves. In general, the calculated heating factors were in good consistence with < /W, which is suggested by ICNIRP-2020 guideline (ICNIRP 2020). In this study, the uncertainty from the incidence angle was not considered. Owing to the rotational symmetry of the multilayer models, the results of our 2D simulation were independent of the incident angle of the exposed plane wave. Nakae et al (2020) also found that the heating factors for the APD are insensitive to the incident angle. In conclusion, the averaged APD and temperature rise in body models with nonplanar surfaces under EM exposure at > Acknowledgement
This work was partly supported by the Ministry of Internal Affairs and Communications, Japan. eferences
Alekseev S I, Radzievsky A A, Logani M K and Ziskin M C 2008 Millimeter wave dosimetry of human skin
Bioelectromagnetics TM -2019 βIEEE standard for safety levels with respect to human exposure to electric, magnetic, and electromagnetic fields, 0 Hz to 300 GHzβ IEEE Access IEEE Trans. Microw. Theory Tech. IEEE Antennas Wirel. Propag. Lett. IEEE Access IEEE Trans. Electromagn. Compat. Phys. Med. Biol. Health Phys.
IEEE Access Phys. Med. Biol. Phys. Med. Biol. IEEE Access IEEE Trans. Biomed. Eng. Ann. Telecommun. IEEE Trans. Electromagn. Compat. Phys. Med. Biol. Health Phys. pp 1β4 Kageyama I, Masuda H, Morimatsu Y, Ishitake T, Sakakibara K, Hikage T and Hirata A 2019 Comparison of temperature elevation between in physical phantom skin and in human skin during local exposure to a 28 GHz millimeter-wave pp 766β9 Kanezaki A, Hirata A, Watanabe S and Shirai H 2010 Parameter variation effects on temperature elevation in a steady-state, one-dimensional thermal model for millimeter wave exposure of one- and three-layer human tissue
Phys. Med. Biol. Phys. Med. Biol. Phys. Med. Biol. IEEE Commun. Mag. Phys. Med. Biol. Phys. Med. Biol. Bioimaging IEEE Access Bioelectromagnetics Phys. Med. Biol. J. Appl. Physiol. Microw. Opt. Technol. Lett. Phys. Med. Biol. N221--N229 Samaras T and Kuster N 2019 Theoretical evaluation of the power transmitted to the body as a function of angle of incidence and polarization at frequencies >6 GHz and its relevance for standardization
Bioelectromagnetics Phys. Med. Biol. Med. Phys. Computational Electrodynamics: The Finite-Difference Time-Domain Method
IEEE Microw. Mag. Phys. Med. Biol. IEEE Trans. Antennas Propag. Bioelectromagnetics173β89