nova página do texto(beta)
Inglês (pdf)
Artigo em XML
Referências do artigo
Tradução automática
Enviar este artigo por email
Citado por SciELO 
Similares em
SciELO 
Permalink
Introduction
Ultrasounds are mechanical waves with frequencies above 20 kHz that interact with media while propagating through them. This interaction produces phenomena such as reflection, refraction, scattering, attenuation, changes in sound velocity, that may modify the acoustic pattern[1]. The ultrasonic beam pattern is strongly dependent on the geometry of the vibrating element, the operating frequency, the transducer manufacturing, among others[2]. Ultrasonic waves energy can be concentrated in a well-located target or region named focus; this pattern modification is called focused ultrasound (FUS). The focalization can be achieved by different ways such as acoustic lenses, spherically concave transducer, electronic or phase-controlled arrays, etc.[3] [4] [5].
Focused ultrasound as being a non-ionizing energy and having the capacity to deliver it in a specific region has been under research for clinical applications for the last decades. FUS interaction with biological tissue can cause thermal effects[6], nonthermal effects such as cavitation[4] [6], nonlinear behavior like both focus distortion and shift[7] [8] and bioeffects[9] [10]. In cancer treatment, the main goal is to achieve a thermal effect due to FUS by increasing the target temperature over 60 °C which results in coagulative necrosis of tissue[4] [10]. Extracorporeal and intracavitary devices based on FUS are mostly designed to treat malignant tumors in liver, kidney, prostate, breast, osteosarcoma, etc.[3] [4] [5] [11]. FUS in cancer treatment acts as adjuvant method to radiotherapy, immunotherapy and chemotherapy[12]. Yet, FUS effects are being studied in brain malignancies like Parkinson disease, to open the blood brain barrier, in the treatment of varicose veins, and cancer treatment, among others[13] [14] [15] [16].
Computational modeling is a powerful tool used to solve multi-physic scenarios and get an approximate behavior of a real problem. The adequate selection of the physics involved in the problem, the geometry definition, physics configuration and boundary conditions settings influence the numerical solution. The finite element method (FEM) allows acoustic propagation modeling; but it can be challenging when solving high-frequency large geometries, i.e., tridimensional spaces[17] [18]. In acoustics, FEM requires both domain and subdomains discretization with a mesh containing frequency-dependent element size; at least 10 elements per wavelength (λ)[19] [20] [21]. In large or multi-domain geometries, this results in an increase of the number of degrees of freedom which implies more computational resources and longer solution time[17] [22]. On the other hand, the boundary element method (BEM) allows faster computation of acoustic radiation patterns in large or infinite homogeneous domains, called exterior problems [23] [24] [25]. BEM modeling requires surface or boundary discretization[21] [23] the problem solution is obtained by using the free-space Green’s function which satisfies the Sommerfeld radiation condition at infinity [21] [26]. However, in multidomain geometries, BEM is limited to obtain the scattering effect over the surface of a subdomain located in the radiation pattern, but not the propagation within it: interior problem Ω[23]. Even though BEM can be used to solve interior problems or the acoustic propagation in a finite domain, FEM solution is more approximated to the real problem[23][27]. Meanwhile, when it refers to obtaining the solution of the exterior problem Ωe, FEM has limitations due to infinite domain discretization[23]; cost-effectiveness comparison between both methods has been analyzed by Harari, I. and Hughes, T.[27]. Nonetheless, both methods present advantages in time-harmonics acoustic model solution whether is the interior problem or the exterior problem. In this sense, Assaad et al.[28] combined FEM solution with BEM to compute the radiation pattern of a circular transducer in a 3D geometry and its 2D axisymmetric approach. In order to do so, the authors modeled the piezoelectric vibration by means of FEM, while the load and propagation in water was modeled by BEM. The electrical impedance modulus and the far-field directivity obtained from the 3D model and the axisymmetric model show an identical pattern. Furthermore, the electric impedance of 3D model of a square shape transducer was compared with measurements that showed good agreement, except in the impedance peaks amplitude.
In the biomedical field, Santiago, A. et al. compared BEM and FEM methods to estimate the elastic modulus from two generated images including a lesion: one without pressure applied and one after a known pressure was applied[29]. By using the same quantity of elements to mesh the geometries, FEM results showed a higher mean error in elastic modulus estimation of the lesion than BEM results. Computational time to obtain the solutions by FEM was around 6 h compared to 60 s by using BEM[29]. Recently, Shen, F. et al.[18] combined BEM and FEM to simulate the transcranial acoustic propagation in 2D and 3D models. For the 2D model, the transducer was a linear array with 512 elements; while for the 3D model, the ultrasonic transducer consisted of 256 source points. Both models were solved at 0.5 MHz. The implemented algorithms corrected the focus shift due to aberration and attenuation caused by the skull bone[18].
In ultrasonic thermal ablation applications, Gélat, P. et al.[30] used BEM to model both acoustic and scattered fields of a high intensity focused ultrasound (HIFU) phased array through human ribs for liver and pancreatic tumors treatment. The transducer consisted of 256 circular elements placed over a spherical bowl shape with a focal distance of 18 cm, each element had a 3 mm radius, and a frequency of 1 MHz. The study was performed considering two possible scenarios of beam focusing at 3 cm depth from the rib cage: intercostally (between ribs 10-11) or trans-costally (across rib 10). The total acoustic pressure or surface velocity of each element was average to attain a focal peak pressure of 10 MPa in free space. For the intercostal case, for the desired focal pressure, the deposited acoustic pressure in the surface of the ribs reached 2.2 MPa, while it reached 1.8 MPa for the transcostal case. In both cases, that amount of acoustic pressure could lead to both bone and tissue heating. Afterwards, Gélat, P. et al.[31] solved the inverse problem by adding constraints to their previous work in order to achieve a desired acoustic pressure on the surface of the rib like geometry. Among the constraints, the element velocity magnitudes were defined in a specific dynamic range. The dimensions of the geometry were modified to work with a frequency of 100 kHz, this way a faster solution could be obtained. The results showed both the velocity magnitude and phase distributions of the 256 elements considering the constraint which could lead to the desired acoustic pressure deposition in the rib like geometry. Furthermore, Van’t Wout, E. et al.[17] implemented a numerical preconditioner to reduce the converge time of the BEM formulation developed in[30] [31] for targeting tissue behind the ribcage. The improved algorithms solved a realistic transcostal simulation of HIFU propagation within 2 h by using a desktop computer with 12 processors and 80 GB RAM [17]. A more complex model for HIFU field radiation through the ribcage based on BEM was reported by Haqshenas, S.R. et al.[22]; which included an idealized subdomain with fat layer properties besides two ribs. The HIFU transducer geometry, the operating frequency, focal distance and ribcage position were the same as reported in[30]; the fat layer was placed between 1 mm to 15 mm from the ribs facing the transducer[22]. The reported results showed that the maximum pressure reached at the focus diminished in 15 % considering transmission and scattered field through the ribs compared with the maximum pressure without scatters. After adding the fat layer before the ribs, the focal pressure showed a decay of 50 %. Moreover, Haqshenas, S.R. et al.[22], implemented an algorithm to couple the exterior problem solution with the interior problem calculation to model the HIFU beam inside a kidney subdomain surrounded by a fat layer. Numerical modeling results showed that the focal pressure decreased 5 % in the absence of the fat layer, and after considering the fat subdomain, the focal pressure within the kidney decreased 15 %. Computational time for solving the two anatomical models (fat-ribs and fat-layer) varied from 82 min to 40 h approximately in a workstation with 32 processors and 512 GB RAM[22].
Motivations of the study
The aim of this work is to compare the acoustic pattern of a 3D focused radiating surface in a homogeneous lossless medium by means of both FEM and BEM towards the design of biomedical ultrasonic extracorporeal applicators intended for thermal ablation therapies. The proposed geometry represents a high computational cost, and its shape differs from the typical geometries used in HIFU applications. Therefore, a simplified medium was considered because our priority was to evaluate the performance of both methods in terms of ultrasonic radiation pattern characteristics and solution time. Besides, its potential clinical application could reach areas where a high focal acoustic energy is required at medium or high penetration depths. Additionally, the focal dimension should be capable of attaining targets smaller than 20 mm with an applicator of, relatively, both reduced dimensions and fabrication costs. Nowadays, to our best knowledge, in the ultrasonic ablation field, BEM use has been restricted to liver ablation through the ribcage with 256 elements placed in a spherical bowl applicator and there are few studies oriented to compare BEM and FEM performances in the biomedical field.
Materials and methods
In order to achieve the goal of this work, four tridimensional time-harmonic acoustic models for a defined radiating surface were developed. Two models consisted of solving the acoustic pattern of the entire radiating surface; while the remaining models considered only half of the geometry, i.e. by defining a symmetrical plane. Moreover, FEM and BEM models were performed for each case: entire geometry and symmetric geometry, as described in Table 1. The symmetric cases represent a simplified version of the entire geometry cases that reduce computational resources. The radiating surface SR proposed in all four cases simulates a mono-element transducer with the shape of a cylindrical section with a radius R = 15 mm, a width b = 6 mm and the length of the chord of the cylindrical boundary a = 21.22 mm (see Figure 1A). The origin coordinates (x, y, z) for all models were set at (x=0, y=0, z=0), as shown in Figure 1A. The models were solved by using the acoustic module of COMSOL Multiphysics 5.4® on a PC with Intel® Xeon® W-2133 CPU @ 3.60 GHz processor and 128 GB RAM.
Table 1 Model identifier based on geometry type and computational method used to obtain the solution.
| Geometry | FEM | BEM |
|---|---|---|
| Entire | Case A | Case C |
| Symmetrical | Case B | Case D |
Figure 1 Geometries designed to solve the FEM and BEM models. A) Proposed geometry of the radiating surface, and Case C: entire radiating surface SR emitting in an infinite void with water acoustic properties to solve the exterior problem, B) Case A: cylindrical domain with entire radiating surface SR emitting in a finite domain: interior problem, C) Case B: symmetrical geometry obtained from cutting the geometry of case A in half used to solve the interior problem, D) Case D: symmetrical geometry obtained from cutting the geometry of case C in half used to solve the exterior problem. The dimensions of the radiating surface are the same in all four models; the geometrical origin was located at coordinates (x=0, y=0, z=0). The propagation medium was set as water for both interior and exterior problems, i.e. FEM and BEM models, respectively.
FEM acoustic radiation modeling
Figures 1B and 1C show the geometries used to model the acoustic beam pattern in the interior problem, where Ωi represents the finite domain, i.e. cases A and B. The time-harmonic wave equation for pressure acoustic propagation used for FEM modeling is represented by Equation (1) as
where p is the acoustic pressure (Pa), ω is the angular frequency (rad/s), and cs is the sound velocity (m/s). The domain consists of a cylinder of 15 mm of radius and a height of 10 mm. The radiating surface center is placed in the middle of the lateral side of the cylinder, i.e. there is a 2 mm distance from SR to each base of the cylinder. The domain acoustic properties were defined as to emulate water which is considered as a lossless medium[32]; the sound velocity was set at 1500 m/s and the density at 1000 kg/m3. The boundary conditions used to solve the model were configured as acoustic impedance, Z (Rayls), in both bases and lateral side of the cylindric domain (boundaries 1, 2 and 3; see Figure 1BB). The acoustic impedance is defined by Equation (2) as
where ρ is the medium density (kg/m3); therefore, Z = 1.5 MRayls which corresponds to water to decrease wave reflection[32]. The operating frequency was set at 1 MHz. Besides, the boundary condition of the SR was set as acoustic pressure with an arbitrary initial value of 500 Pa.
In case B, as the geometry used in the model is half of the geometry in case A, a symmetry boundary condition was configured in the cross-sectional plane of the cylindric domain (boundary 4; see Figure 1C). Boundaries 1, 2 and 3 remained as impedance condition with water acoustic properties. FEM models were solved by using a mesh consisting of tetrahedral elements of quadratic order. A convergence analysis was realized in both cases; therefore, the element size was varied from λ/6 up to λ/12 for case A and from λ/6 up to λ/14 for case B. The resulting meshes consisted of 7, 942, 317 elements and 63, 823, 707 elements, respectively, for case A, and 3, 946, 491 elements and 50, 797, 081 elements, respectively, for case B. For case A, the solution time for λ/12 mesh was 1 day 2 h 32 min 15 s (123.74 GB physical memory/192.68 GB virtual memory); for λ/6 mesh was 14 min 39 s (53.12 GB physical memory/58.34 GB virtual memory). For case B, the solution time for λ/14 was 14 h 19 min 9 s (123.8 GB physical memory/165.96 GB virtual memory); for λ/6 was 5 min 55 s (23.22 GB physical memory/25.5 GB virtual memory). Both cases were tried to be solved with finer meshes, but computational resources were not enough.
BEM acoustic radiation modeling
The geometries depicted in Figures 1A and 1D were used to solve the acoustic beam pattern by means of BEM. The main geometry in the models represents the entire radiating surface or half of it; and the dimensions b, a and R are the same as in cases A and B. In this scenario, there is no need to define a finite domain, because BEM is used to solve the exterior problem Ωe; however, the infinite void (space) that surrounds SR has the acoustic properties of water. The governing Helmholtz equation defined for BEM modeling is given by Equations (3), (4) and (5) as follows
where pt is the total acoustic pressure (Pa), keq is the wave number (rad/m), ρc is the complex density of the medium (kg/m3), cc is the speed of sound (m/s) that could have complex values depending on the problem, pb is the background pressure wave (Pa) and p is the pressure (Pa). For both cases C and D, cc = c, ρc = ρ, and pt = p; after substituting the previous values in eqs. (3)-(5), the governing equation results in eq. (1). Similarly to cases A and B, the boundary condition of SR was set as acoustic pressure with an arbitrary value of 500 Pa. In case D, half of the radiating surface is modeled; therefore, in axis y = 0, an infinite symmetry plane was configured as shown in Figure 1D. For both cases C and D, the condition at infinity was set as an outgoing wave, in order to satisfy the Sommerfeld radiation condition[21][26]. As well as in FEM models, the operating frequency was set at 1 MHz.
BEM models were solved by using a mesh consisting of triangular elements of quadratic order. A convergence analysis was performed in both cases by varying the boundary element size from λ/6 up to λ/20. The resulting meshes for case C consisted of 5, 472 elements and 59, 244 elements, respectively, and 200 elements and 29, 724 elements, respectively, for case D. For case C, the solution time for λ/20 mesh was 1 min 20 s (3.75 GB physical memory/4.41 GB virtual memory), and for λ/6 mesh was 10 s (1.61 GB physical memory/2.2 GB virtual memory). For case D, the solution time for λ/20 mesh was 58 s (2.56 GB physical memory/3.09 GB virtual memory), and for λ/6 mesh was 12 s (1.31 GB physical memory/1.89 GB virtual memory).
BEM model validation
After performing the convergency analysis, BEM results were validated with FEM results by using the models solved with a λ /10 mesh as reference[19] [21]. Then, focus length and width at -3 dB pressure decay, full width at half maximum (FWHM) pressure, and full length at half maximum (FLHM) pressure were calculated. Contours at -3 dB and -6 dB decay and their corresponding areas were obtained in zy-, xy-, and zx- planes for all cases. The focal region produced by high-intensity focused ultrasound (HIFU), typically, has an ellipsoidal shape[4]; however, its shape depends on the radiator geometry and its frequency, among other parameters[2] [33]. Therefore, an elliptical shape ratio Er was calculated in the three planes of each model to obtain the similarity percentage of the contours to an ellipse value; where Er = 100 % is the ideal[32]. Finally, ellipsoidal approximations of focal volumes at -3 dB and -6 dB were estimated by means of the major and minor axis of the contours at zy- and xy- planes.
Results and discussion
The convergence analysis was carried out by obtaining the maximum acoustic pressure Pmax value along zy- plane of each model. Then, the relative error was calculated by taking the solution of the model meshed with an element size of λ/10 as reference[32]; where Pmax started to show stability for each case (see Table 2). For case A, the model was solved for meshes from λ/6 up to λ/12, due to available computational resources, the relative errors obtained were 0.061 % and 0.032%, respectively. For case B, the model was solved for meshes from λ/6 up to λ/14, due to computational resources, the relative errors calculated were 0.021 % and 0.057 %, respectively. For case C, the simulations were done from λ/6 up to λ/20, the relative errors calculated were 0.111 % and 0.064 %, respectively. Finally, for case D, the simulations were done from λ/6 up to λ/20, and the relative errors obtained were 4.141 % and 0.126 %, respectively. From now on, the results presented correspond to the models solved with a λ/10 mesh for all cases.
Table 2 Convergence analysis results for the four cases. Solution from λ/10 mesh models was used as reference for calculating the relative error.
| Case | Mesh size | No. elements | DOF | Solution time | Relative error [%] |
|---|---|---|---|---|---|
| A | λ/6 | 7942317 | 10673982 | 14 min 39 s | 0.062 |
| λ /8 | 18862851 | 25293358 | 44 min 1 s | 0.035 | |
| λ /10 | 36888845 | 49391271 | 9 h 2 min 7 s | - | |
| λ /12 | 63823707 | 85374660 | 1 day 2 h 32 min 15 s | 0.032 | |
| B | λ /6 | 3946491 | 5314575 | 5 min 55 s | 0.022 |
| λ /8 | 9394791 | 12616975 | 14 min 34 s | 0.005 | |
| λ /10 | 18402298 | 24664459 | 40 min 4 s | - | |
| λ /12 | 31828411 | 42619267 | 5 h 37 min 9 s | 0.038 | |
| λ /14 | 50797081 | 67965941 | 14 h 19 min 9 s | 0.057 | |
| C | λ /6 | 5472 | 14191 | 10 s | 0.111 |
| λ /8 | 9660 | 24821 | 14 s | 0.051 | |
| λ /10 | 15216 | 38883 | 21 s | - | |
| λ /12 | 21360 | 54403 | 28 s | 0.013 | |
| λ /14 | 29340 | 74525 | 38 s | 0.034 | |
| λ /16 | 37644 | 95445 | 48 s | 0.047 | |
| λ /18 | 48048 | 121627 | 1 min 2 s | 0.059 | |
| λ /20 | 59244 | 149777 | 1 min 20 s | 0.064 | |
| D | λ /6 | 200 | 577 | 3 s | 4.141 |
| λ /8 | 4844 | 12495 | 11 s | 0.014 | |
| λ /10 | 7540 | 19331 | 16 s | - | |
| λ /12 | 10756 | 27467 | 22 s | 0.077 | |
| λ /14 | 14676 | 37363 | 29 s | 0.050 | |
| λ /16 | 18864 | 47925 | 36 s | 0.118 | |
| λ /18 | 23988 | 60831 | 46 s | 0.138 | |
| λ /20 | 29724 | 75267 | 58 s | 0.127 |
Figure 2A depicts the pressure distribution along z- axis for all cases and it is observed that the acoustic pressure reaches its maximum near 15 mm. The acoustic pressure difference ΔPmax at its maximum value was 41.101 Pa between cases A and C which correspond to the models with the entire radiating surface, which represents an error of 1.48 % taking as reference the case A (FEM complete geometry). The distance difference Δz between Pmax of both cases was 0.451 mm. For the symmetrical models, cases B and D, ΔPmax was 38.649 Pa that represents an error of 1.37 % taking as reference case B (FEM symmetric geometry) and Δz was 0.452 mm. Figures 2B and 2C show the pressure distribution along y- axis and x- axis at the z position of Pmax, where Pmax is centered at x = 0 and y = 0 in accordance with the proposed geometry. Figure 2B depicts that the pressure distribution along y- axis is narrower than in the x- axis, see Figure 2C. This behavior is strongly related to the geometrical shape of the surface radiator [2] [33], presenting an oval shape in xy- plane at z position of Pmax (see Figure 3) with its mayor axis over x- axis and differentiating it form the classical HIFU focal shape.
Figure 2 Pressure distribution (A) along z- axis, (B) along y- axis, and (C) along x- axis. All three pressure distributions were taken at the maximum acoustic pressure value Pmax obtained from zy- plane.
Figure 3 Pressure contours at -3 dB and -6 dB decays for all four cases. A) -3 dB contours in zy- plane, B) -3 dB contours in xy- plane, C) -3 dB contours in zx- plane, D) -6 dB contours in zy- plane, E) -6 dB contours in xy- plane, and F) -6 dB contours in zx- plane.
In this sense, it is worth noting that the focal region dimensions in each plane are also different as shown in Figure 3. In Figures 3A-3C depicts the contours at -3 dB Pmax decay in zy-, xy-, zx- planes for all four cases. In Figure 3A, the contours show an elliptic focal shape, where case A and case B present at slight displacement towards the left (surface radiator location). These displacements are related with the mentioned Δz between cases A and C, and cases B and D; while, in Figures 3B and 3C, all contours are centered in x = 0 and y = 0; i.e. no error was found over x and y positions of Pmax. Figures 3D and 3F show the contours at -6 dB pressure decay in zy-, xy-, zx- planes. From Figure 3D, it can be observed an enlarged region towards the left (SR location) which describes the beam path through focalization. Also, isolated contours can be appreciated for cases A and B (FEM modeling) but not for cases C and D (BEM modeling). In a clinical application, such as FUS thermal ablation, the enlarged area and the appearance of the isolated contours could produce “hot spots” or undesired hot regions[34] at the zone nearby the transducer and its prediction is a fundamental step to avoid patient damage. Figure 4 depicts the pressure distribution along y- axis at the position in z- axis where the isolated contours are presented (z = 10.82 mm). The pressure peaks reach a maximum value of 1456.03 Pa for case A, and 1445.33 Pa for case B, that represent half decay of Pmax. Focal dimensions are usually measured at half decay of P[35] which could directly relate to the size of the produced thermal lesion.
Figure 4 Pressure distribution along y- axis at z = 10.82 mm, location of the two isolated contours.
However, the thermal lesion size depends on the exposure time to ultrasound energy, cross-sectional intensity accumulation, tissue absorption coefficient, non-linear effects, frequency, among others[35]. In Figure 3E it can be observed that the focal regions keep an elliptical fashion as in Figure 3B. In Figure 3F, the contour shows an enlargement along z- axis from 14 to 18 mm and along x- axis which could produce non-uniform thermal lesion, and care must be taken in order to avoid healthy tissue damage.
Table 3 shows the lengths and widths in zy-, xy- and zx- planes calculated at -3 dB and -6 dB (FLHM and FWHM) pressure decays for all four cases. As mentioned before, both solution time and memory resources increase when solving high-frequency large FEM models; therefore, a model simplification must be pursued to reduce computational cost without compromising the approximated solution. Hence, results from cases A and B (FEM) were taken as references to evaluate cases C and D (BEM) respectively. Case B represents the simplification of case A, which could be considered as the more realistic model in this study; as well as Case D represents simplification of case C. Then, the most representative ratio (i.e. equal to one when ideal) in length (lengthBEM-Case (C or D)/lengthRef-Case (A or B)) and width (widthBEM-Case (C or D)/widthRef-Case (A or B)) were found at -3 dB and -6 dB pressure decays. In both cases (C/A and D/B), the worst ratio values were found at zx-plane width (0.92 for C/A and D/B at -3 dB and 0.90 for C/A and D/B at -6 dB) and at xy-plane length at -6 dB (0.94 for C/A and D/B), representing errors up to 9.87 %. All other ratios approximate unity with an error minor to 1.7 %. Additionally, it is important to address that both lengths in the zy- and zx- planes should agree as both distances represent the focus extent in the propagation axis z. In the same way, the widths in both zy- and xy- planes must have close values because they correspond to the focal width along y- axis. In consequence, the length in xy- planes and the width in zx- planes present close values due to its relationship with the dimension in the x- axis (see Figure 3 and Table 3).
Table 3 Focal length and width in zy-, xy- and zx- planes calculated at -3 dB and -6 dB (FLHM and FWHM) pressure decays for the four cases.
| Pressure level |
Dimension | Plane | Case A | Case B | Case C | Case D |
|---|---|---|---|---|---|---|
| -3 dB | Length [mm] | zy | 4.2586 | 4.2296 | 4.2909 | 4.2631 |
| xy | 3.5695 | 3.5732 | 3.5743 | 3.5751 | ||
| zx | 4.2270 | 4.2270 | 4.265 | 4.2639 | ||
| Width [mm] | zy | 0.8914 | 1.0707 | 0.8932 | 1.0730 | |
| xy | 0.8914 | 1.0712 | 0.8927 | 1.0724 | ||
| zx | 3.8963 | 3.8955 | 3.5872 | 3.5862 | ||
| -6 dB | FLHM [mm] | zy | 6.0302 | 6.0302 | 6.0302 | 5.9296 |
| xy | 5.4097 | 5.4156 | 5.0728 | 5.0742 | ||
| zx | 5.9712 | 5.9704 | 6.0343 | 6.0342 | ||
| FWHM [mm] | zy | 1.2121 | 1.2121 | 1.2121 | 1.2121 | |
| xy | 1.2270 | 1.4011 | 1.2140 | 1.3933 | ||
| zx | 5.6680 | 5.6689 | 5.1098 | 5.1093 |
Table 4 shows the calculated areas within the contours at -3 dB and -6 dB pressure decays for all four models. The areas from cases A and C have in all three -3 dB contours and xy- plane in -6 dB contour similar values. This can be due to the geometry type used to model the radiation pattern, i.e., the model includes the entire radiating surface. However, -6 dB contours areas in zy- and zx- planes present differences of 1.62 mm2 and 2.272 mm2, respectively, due to the isolated contours presence in the zy- plane and the wider extension of the focus in zx- plane (see Figures 3D-F). As well as before, areas from cases B and D (half of the geometry) show differences of 1.259 mm2 and 2.095 mm2 in the -6 dB contours at zy- and zx- planes, respectively. In both geometries, entire and half, the larger areas and focus distortions were found at - 6 dB contours FEM solutions in the zy- and zx- planes.
Table 4 Areas within contours at -3 dB and -6 dB decays of P max .
| Areas within contours |
Plane | Case A [mm2] |
Case B [mm2] |
Case C [mm2] |
Case D [mm2] |
|---|---|---|---|---|---|
| -3 dB | zy | 3.175 | 3.856 | 3.197 | 3.896 |
| xy | 2.543 | 3.134 | 2.553 | 3.147 | |
| zx | 12.837 | 12.844 | 12.086 | 12.327 | |
| -6 dB | zy | 11.031 | 11.420 | 9.411 | 10.161 |
| xy | 5.089 | 5.988 | 4.932 | 5.794 | |
| zx | 26.640 | 26.650 | 24.368 | 24.555 |
Table 5 shows the ellipsoidal shape ratio E[32] obtained in the four cases. In all cases, the shape ratio calculated is higher than 92 %, except in -6 dB contours of the zy- planes due to the presence of the isolated contours in cases A and B, and the enlarged area located between 11-13 mm in z- axis (see Figure 3D). This behavior can be avoided by modifying the SR geometrical characteristics. Even though the typical focal region shape is ellipsoidal[33], there are ultrasonic transducers that could produce different focal shapes such as conical[36], annular[37] or even two simultaneous foci[38]. In addition, it is important to remark that all four cases have a similar percentage of shape ratio in their corresponding planes which approximates to the typical focal shape[33] that was taken as reference.
Table 5 Elliptical shape ratio percentage.
| Contour level | Plane | Case A [%] | Case B [%] | Case C [%] | Case D [%] |
|---|---|---|---|---|---|
| -3 dB | zy | 93.76 | 92.10 | 94.10 | 92.36 |
| xy | 98.16 | 95.90 | 98.16 | 95.67 | |
| zx | 99.33 | 99.34 | 99.34 | 99.34 | |
| -6 dB | zy | 52.03 | 49.42 | 61.00 | 56.62 |
| xy | 97.63 | 99.50 | 98.05 | 95.81 | |
| zx | 99.74 | 99.75 | 99.31 | 99.32 |
Table 6 shows the ellipsoidal volume approximations estimated at -3 dB and -6 dB pressure decay levels. For the calculation, half-length of the contour in the zy- plane and both half-length and half-width of the contours in the xy- plane were considered as the ellipsoid semi-axis. By making this assumption, the enlarged regions in the left side of Figure 3D are ignored; therefore, further analysis in the contribution of these areas should be performed. In the -3 dB contour level, the volume difference between cases A and C was 0.0738 mm3, while between cases B and D was 0.0817 mm3. In the -6 dB contour pressure decay, the volume difference was 1.514 mm3 and 1.2354 mm3 between cases A and C, and B and D, respectively.
Table 6 Ellipsoidal volume approximation at -3 dB and -6 dB decays of Pmax.
| Contour level |
Case A [mm3] |
Case B [mm3] |
Case C [mm3] |
Case D [mm3] |
|---|---|---|---|---|
| -3 dB | 7.0951 | 8.4766 | 7.1689 | 8.5583 |
| -6 dB | 20.9575 | 23.5579 | 19.4435 | 22.3225 |
The larger differences between FEM and BEM areas and volumes, e.g. at -6 dB pressure decay contours, may be related to spurious reflections produced by the subdomain cylindrical boundaries and the impedance condition used to limit the interior problem. While BEM solutions account for no-reflections as they were modeled in an infinite or unbounded space. Nevertheless, these differences are reasonable small and validate the use of BEM to model the acoustic pattern of a FUS biomedical applicator prior to solving the interior problem. Additionally, the analyzed data suggest that FEM model needs improving, such as modifying the propagation medium shape and dimensions, and the inclusion of perfectly matched layers (PML). PMLs simulate an infinite domain in the exterior boundaries of the domain and reduces unwanted waves reflections[21] [24].
Conclusions
Four acoustic propagation models based on the same radiating surface geometry are presented in this work. The models consisted of an entire or a symmetric radiating surface propagating in a lossless medium. The results were obtained by means of FEM and BEM, for the interior Ω and exterior Ωe problems solution, respectively. The analyzed data by using the four models solved with a λ/10 mesh show that the average ΔPmax and the average Δz reached at the focal spot are 39.875 Pa and 0.4515 mm, respectively. Low values of errors were computed for Pmax estimated with BEM cases, and taking as references the FEM cases: 1.48 % for entire geometries and 1.37 % for symmetric ones.
Isolated contours related to undesirable hot spots that were identified with the FEM models were not estimated using BEM models. This could be due to the computation method of BEM, which has reached better results on the named “exterior problem” than on the “interior one”. A future work could be oriented to assess and solve these “erroneous estimations” of BEM.
As for dimensional assessment, almost all of the length and width ratios (BEM/FEM) approximate unity with an error minor to 1.7%, except in zx-plane width at -3 dB and -6 dB, and at xy-plane length at -6 dB, presenting errors up to 9.87 %. The maximum areas differences between the four cases in each plane at both contours pressure decays vary from 0.6 mm2 to 2.3 mm2; which addresses that symmetric geometries could be used instead of the entire geometries. The Er was over 92 % in most of the planes, except at -6 dB contours of zy- plane, where the average Er was 54.76 %; still, the four models presented a similar beam pattern formation. These results indicate that BEM can be used to solve the acoustic propagation problem in lossless media in a faster and with less computational resources fashion than FEM. In this sense, BEM represents a powerful tool during the design of the beam pattern of ultrasonic radiators intended for thermal ablation therapies. Nevertheless, it is important to mention that BEM presents limitations when solving the interior problem, and FEM continues to be the best method to solve acoustic wave interactions within media. Further work includes geometrical radiator modification to attain deeper penetration distances in a lossless medium, domain sectioning in subdomains and combining meshes size to solve the interior problem, modeling both ultrasonic propagation in biological tissue emulators and tissue heating due to FUS.
Author contributions
R. M. V. conceptualization, data curation, formal analysis, investigation, methodology, resources, validation, visualization, writing - original draft, writing - review, and editing; I. B. project administration, resources, supervision, validation, writing - review, and editing.
