GamBet shielding calculation walkthrough

In a previous article, I discussed the basis of the H*(10) standard for radiation safety. Here, I’ll continue by demonstrating techniques for shielding calculations. I’ll use a published simple benchmark calculation to compare the results to several other Monte Carlo codes and to illustrate some of the advantages of GamBet. For this, I am grateful to the people at RayXpert who have made reports on a number of benchmark calculations available on the Internet.  In particular, I admire their patience and fortitude in learning four other codes besides their own to make the comparisons.

Geometry of the benchmark calculation

Figure 1. Geometry of the benchmark calculation taken from the RayXpert report.

Figure 1 shows the geometry taken from the report H*(10) dose rate evolution with some different shielding. Targets are located 100 cm from a Co60 point source with activity 10^9 Bq[1]. The dose at a depth of 1.0 cm was calculated for a bare target and with three shield configurations. The setup illustrates two points I made in the last article about H*(10) calculations:

  1. The results are not sensitive to the detector geometry. In this case, a 4.0 cm radius sphere of tissue equivalent material is used rather than the ICRU 15.0 cm radius sphere.
  2. The radiation weighting factor (Wr) is taken as unity for the photons and secondary electrons and there is no organ correction for the phantom (Wt = 1.0). In this case, the biological dose (Sv) equals the physical dose (Gy).

I limited my calculations to the bare target and to one behind a 2 cm lead shield. Here are the results reported by RayXpert. The values are given in units of μGy/hour.

             Dosimex-G  Mercurad  Microshield MCNPX RayXpert
No shielding    349        353         347      355   365
2 cm lead       143        152         141      146   149

The targets subtend only a small fraction of the solid angle surrounding the source, so it would be foolish to use an isotropic source. Instead, this is an ideal opportunity to use the GamBet Radiation Source Tool, illustrated in Fig. 2. After the user enters values in the dialog fields, the program creates a standard source file. The idea is to fill a circle of a specified radius in the direction of interest with a large, random distribution of particles. The flux assigned to each particle is calculated to represent the specified source activity. In this case, I created a distribution of 200,000 photons of energy 1.17 MeV and an equal number at 1.33 MeV. The distribution filled a circle of radius 4.0 cm at a distance of 92.0 from the source, ensuring complete irradiation of the detector. The energy deposition was symmetric about the line from the source through the detector, so I further reduced run time by using the option for cylindrical scoring in GamBet. The air region was represented by the standard Penelope dry air model and the tissue-equivalent target was represented by the custom material definition discussed in the previous article:

Material
 Name TissueEquivalent
 Component O 1.0000
 Component C 0.1940
 Component H 2.1032
 Component N 0.0390
 Density 1.00
 Insulator
End
GamBet Radiation Source Tool dialog

Figure 2. GamBet Radioactive Source Tool dialog.

Figure 3 shows the dose distribution calculated by GamBet for the bare detector with photons entering from the left[2]. Even for this simple case, there is quite a lot going on. The dose at the surface is lower than the peak dose in the material. This is the result of the buildup of the secondary-electron distribution. The energy loss-rate for electrons is much higher than that of photons of the same energy. The depth to reach the equilibrium dose is about 0.6 cm, hence the choice of the H*(10) measurement depth as 1.0 cm for γ rays in the MeV range. The dose rate decreases gradually through the material because of γ ray attenuation. Note the difference in the air dose rate at the entrance and exit of the target. The enhanced downstream dose rate results from electrons driven out of the material. The H*(10) dose rate determined by GamBet was 338 μGy/hour, in good absolute agreement with the values reported by RayXpert.

Dose-rate disttribution for a bare target

Figure 3. Dose-rate distribution for a bare target.

Figure 4 shows the dose-rate distribution with 2.0 cm lead shield placed a short distance in front of the target. The shield extends to the outer radius — the portion irradiated by the 4.0 cm radius beam is clearly visible. The dose buildup depth in lead is shorter than the element size, so it is not visible[3]. Also, note the absence of a buildup region in the target. Most of the equilibrium secondary electron distribution from the lead crosses the small air gap and enters the target. The GamBet H*(10) prediction is 146 μGy/hour, again comparable to the values determined by the other codes. Finally, Figure 5 shows dose-rate scans along the diameter of the spherical target for the two cases. The dashed line is the H*(10) point.

Dose-rate distribution for a 2.0 cm lead shield

Figure 4. Dose-rate distribution for a 2.0 cm lead shield.

Dose-rate scans along a diameter of the spherical target

Figure 5. Dose-rate scans along a diameter of the spherical target.

One conclusion of this study is that all radiation-shielding codes give reasonable results. This is not surprising because the physics and mathematics base for Monte Carlo radiation codes is over half a century old. The question within this constraint is: why choose GamBet? Beyond its low price, GamBet has several advantages, some of which are apparent in this study:

  • Zoning in GamBet is performed with conformal finite-element meshes, allowing plots such as those of Figure 3 and 4. Complex processes may enter even the simplest radiation calculations, so detailed pictures of dose distributions are essential for evaluating results.
  • GamBet is fast for two reasons: 1) the code is designed to encourage division of solutions into stages and 2) the code supports efficient parallel processing. For this example, the run time with 1.6 million incident model photons (4 parallel processes with 400,000 photons each) was 7 minutes.
  • The element approach in GamBet is more effective than combinatorial solid geometry to represent complex systems. The advantage is not apparent in this calculation, but is overwhelming when representing something like a human body phantom.
  • Although all Monte Carlo physics engines are about equal for MeV γ rays, Penelope has far more detailed support for low-energy X-rays.
  • GamBet has several advanced features, such as integration of calculated 3D electric and magnetic fields and direct coupling to electron beam and thermal transport programs.

Footnotes

[1] This article reviews the decay scheme of Co60: Modeling radioactive sources with GamBet.

[2] GamBet uses standard units of Gy/s. The conversion factor to μGy/hour is 3.6E9.

[3] The element size in GamBet determines the display resolution and the accuracy for fitting shapes, but has no effect on the physics of the radiation interactions.

[4] Contact us : techinfo@fieldp.com.

[5] Field Precision home page: www.fieldp.com.

H*(10) calculations with GamBet

I recently had an inquiry from a prospective customer. He had investigated several Monte Carlo options and was impressed with one package because the developer emphasized that the software could perform H*(10) calculations. In reality, every Monte Carlo X-ray code has this capability. The only requirement is that the user be aware of the definition. To remove any doubts, I will review the significance and technical background of H*(10) calculations in this article.

In applications to personnel shielding, it is not sufficient simply to know the radiation field at a location. By radiation field, I mean the fluxes of energetic electrons and photons along with their energy spectra. The critical issue is how do the particles interact with tissue and how does the dose (deposited energy divided by mass) build up inside the organism. A further complication is that the same dose may have different biological effects, depending on the radiation and tissue type. The standard unit of physical dose is the gray (Gy), equal to 1 joule/kg of deposited energy. The unit of biological dose is the sievert (Sv) which also has units of J/kg. The difference is that the biological dose includes weighting factors that depend on the radiation type (Wr) and the tissue type (Wt) to indicate the relative potential for biological harm. Dose values are related by

D(Sv) = D(Gy)*Wr*Wt

The radiation weighting factor is Wr = 1.0 for photons and electrons, so we need not worry about it in GamBet calculations. The tissue weighting factor Wt is difficult to estimate, and therefore there is considerable disagreement about the values. For whole-body irradiation, the convention is to take Wt = 1.0. The implication is that for most shielding calculations at energies of interest for medical applications, sieverts are equivalent to grays.

The critical issue is the dose buildup through interactions with tissue. To standardize measurements and calculations, the ICRU (International Commission on Radiation Units and Measurements) defined the Directional Dose Equivalent, H*(d). The quantity applies to radiation moving predominantly in one direction, such as X-rays outside a radiation shield. A different quantity, Hp (the personnel dose equivalent), applies to approximately isotropic radiation fields, like the time that Mr. Spock was trapped in the Propulsion Chamber.

The ideal calculation or measurement defined by the ICRU uses a 30-cm-diameter sphere of tissue-equivalent plastic with a density of 1 g/cm^3 and a mass composition of 76.2%  oxygen,  11.1% carbon, 10.1% hydrogen and 2.6% nitrogen. The quantity H*(d) is the biological dose at a depth d below the surface of the sphere in the direction of the radiation (Figure 1). The most common choice of d is 10.0 mm, hence the term H*(10). The reasoning is as follows. When a flux of gamma rays enters a material, the dose increases moving away from the surface because of the generation of secondary electrons (e.g., Compton electrons). The electrons deposit energy more rapidly than photons. The secondary electron density grows with distance until the production rate equals the absorption rate. At this equilibrium point, the dose reaches a maximum. For gamma rays in the 1 MeV range (typical of radioactive sources) in a material with density 1 gm/cm^3 (e.g., water), the equilibrium depth is a 10.0 mm or less. The following article describes benchmark calculations that illustrate this effect.

Geometry to define the Directional Dose Equivalent

Figure 1. Geometry to define the Directional Dose Equivalent.

To conclude, here are a couple practical suggestions for H*(10) calculations in GamBet. First, in defining custom materials, GamBet follows the Penelope convention using the stochiometric composition (the relative number of atoms in an equivalent molecule) rather than mass fractions. To make the conversion, assume that F0, FC, FH and FN are the stochiometric fractions for oxygen, carbon, hydrogen and nitrogen. Because the quantities are relative, we take F0 = 1.00. Using the atomic weights and mass fraction of 0.762 for oxygen, we can determine the other quantities from these equations:

FC*12.011/15.994 = 0.111/0.762
FH*1.00794/15.994 = 0.101/0.762
FN*14.00674/15.994 = 0.026/0.762

Here’s the result as it would appear in the Composition section of a GamBet script:

Material
  Name TissueEquivalent
  Component O 1.0000
  Component C 0.1940
  Component H 2.1032
  Component N 0.0390
  Density 1.00
  Insulator
End

Finally, note it is not necessary to use a 30 cm sphere of the tissue equivalent material in calculations. Any block of material will give about the same answer, as long as the object has depth and width greater than d.

The next article describes a walkthrough example that gives some good insights into the physics of the H*(10) calculation. The example illustrates several GamBet techniques and new features and provides an opportunity to make comparisons with several other radiation-shielding codes.

Footnotes

[1] Contact us : techinfo@fieldp.com.

[2] Field Precision home page: www.fieldp.com.

Variance reduction techniques in GamBet

The GamBet software suite calculates the interactions of energetic photons, electrons and positrons with matter. In a previous sequence of three articles (Monte Carlo methods versus moment equations: Part A, Part B and Part C), I discussed how Monte Carlo calculations approximate the behavior of large collections of particles. The essence of the Monte Carlo method is to follow the interactions of a finite set of model particles with the assumption that they represent the average behavior of the large set. The advantage is that there is a close connection with the interaction physics, so that a Monte Carlo calculation is often called a simulation. The disadvantage is that it is a statistical process, so the accuracy of predictions is limited by random variations set by the number of model particles or events, N. The accuracy in estimating an average quantity is usually expressed as <x> ± σ, where <x> is the ideal average and σ is the standard deviation set by the finite number of samples. The standard deviation equals the square root of the variance, the average of the squared differences from the mean value. The standard deviation scales as 1/√N. Improving the accuracy of estimates by a factor of 10 requires 100 times the number of model particles. The implication is that high accuracy leads to long run times.

The term variance reduction[1] applies to methods to optimize Monte Carlo calculations to gain an edge on the 1/√N limit. Variance reduction is not a formalized mathematical method, but rather a set of common-sense fixes that can go a long way toward reducing run time. The essence is finding ways not to waste time on particles that will not contribute to critical results. The success and validity of variance reduction depends strongly on user judgements. GamBet uses the Penelope package for interaction physics. Penelope has built-in features for variance reduction that are implemented in GamBet with the following commands:

ENHANCE NReg NSplit [ElecP, PhotP, PosiP, ElecS, PhotS, PosiS]

GamBet is unique among Monte Carlo codes in its use of a finite-element conformal meshes to represent the solution volume. Users can divide the space into a number of regions to represent different materials or sections of an object. The Enhance command improves statistics in a critical region (such as a detector). Particles entering the region are split into NSplit particles with statistical weight 1/NSplit. With the optional string parameters, the operation may be limited to specific particle types (primary or secondary electrons, photons and positrons).

REDUCE NReg NKill [ElecP, PhotP, PosiP, ElecS, PhotS, PosiS]

The Reduce command is the inverse of Enhance. The number of particle entering a specified non-critical region are reduced while increasing their statistical weight.

FORCE [ELEC,PHOT,POSI] [HELAS,...]  Factor [NReg]

This command instructs the program to increase the probability of low-probability interactions like bremsstrahlung emission. The statistical weight of reaction products is decreased to preserve the correct energy balance between reactions.

Beyond the Penelope techniques, GamBet is structured to help achieve short run times.  The code was designed to encourage the division of calculations into manageable segments. For example, an initial segment could address a radioactive source with shielding and collimation, while a second segment could address the interaction of forward-directed radiation with tissue. The segments are connected by an escape file which records the set of model particles that reach the boundary of the first segment. The key to variance reduction in GamBet is filtering and transforming escape files for optimal performance in a following segment. The escape distribution can be modified within a GamBet run with the following command:

ESCAPEFILTER Condition01 Condition02 ...

The conditions are strings like X>0.15, T<5.0E6,…. Particles must meet the combined conditions to be included in the escape file. Conditions may apply to spatial locations, kinetic energy and particle type. The idea is to limit particles to those that will play a role in the following segment and to limit the size of the escape file. Recently we expanded the EscapeFilter conditions to include particle direction: Ux>0.1,Uz>0.95,Ur<0.25,… Here, the quantities are the components of a unit vector pointing along the direction of the velocity. One application is to limit particles to those that are aimed toward a detector or target.

GenDist as a component in a variance reduction strategy

Figure 1. GenDist as a component in a variance reduction strategy.

The GamBet package includes GenDist, a powerful utility to create or to modify large particle distributions. GenDist can act as an additional stage between calculation segments (Figure 1) to optimize particle properties for reduced variance. The basic sequence is to read an escape file, filter or transform the particle parameters and to write a modify file to be used as the source for a following calculation segment. In the past, the operations were controlled interactively by the user in a program window. We have recently added a script capability for autonomous production runs. Here is a summary of the new script commands:

READ FPrefix.SRC
Load an escape file

WRITE Fprefix.SRC
Write a file of transformed particle parameters, applying any filter conditions that have been set.

AXIS [X,Y,Z]
Set a reference axis for evaluating transverse velocity in transformations and filters.

FILTER Condition Value
Set any number of filter conditions. The set of conditions is the same as those in the GamBet EscapeFilter command.

XFORM GENERAL XShift YShift ZShift XRotate YRotate ZRotate
Move or rotate the particles to match the coordinate system of the next computational segment.

XFORM UNIDIST Dist
XFORM NORMPLANE Pos
XFORM CLOSETOLINE HLine YLine
Move particles in ballistic orbits following their velocity vectors. The options are shifts 1) a uniform distance backward or forward, 2) to a plane normal to the current axis or 3) to their positions closest to a line parallel to the current axis. One main application is to find the effective radius of a bremsstrahlung source for X-ray imaging applications by back-projection to the target.

BEAMSECT2D NThet
BEAMSECT3D X 0 X X
The first command converts a distribution from 2D cylindrical calculation to one suitable for a 3D calculation. The second command limits 3D beam distributions to specific transverse quadrants or mirrors a beam distribution.

SCALE Fact
Change the size of a particle beam or convert spatial units to match calculation segments. For example, a source calculation may use units of mm, but units of m may be more appropriate for the detector calculation.

For more details, use these links to download the GamBet and GenDist manuals. Free updates to the new programs are available to GamBet and Xenos users.

Footnotes

[1] The term standard deviation reduction would be a better choice, but it is a less compelling phrase.

[2] Contact us : techinfo@fieldp.com.

[3] Field Precision home page: www.fieldp.com.

 

Moving thermal sources in HeatWave

The generation of X-rays via bremsstrahlung is energetically inefficient at electron beam energies of interest for medical applications. Therefore, high-power electron beams are necessary for procedures such as CT scans. Dissipation of energy to prevent target melting is a major issue. One solution is a rapidly-rotating target. Alternatively, we recently worked with a customer using a liquid jet of gallium alloy as the X-target.

HeatWave, a component of our Xenos suite for X-ray source design, performs 3D thermal conduction calculations. Energy deposition profiles determined by the GamBet Monte Carlo code can be imported as thermal source distributions. The distributions can be modulated with user-specified functions to represent pulsed or periodic beams. This capability can provide, at best, an rough approximation to a rotating target.

The HeatWave solution is performed on a stationary mesh — it would be difficult and unwieldy to introduce a moving mesh. On the other hand, it is relatively easy to move the heat source through the mesh, achieving the same result. Accordingly, we have introduced a new capability in HeatWave to generate exact moving-target simulations. This article summarizes its operation.

Power distributions from the programs Aether (microwave fields), RFE3 (time-dependent electric fields) and GamBet (Monte Carlo X-rays and electrons in matter) can be imported into HeatWave with the SourceFile command. The only restriction is that the field and particle solutions must be performed on the same mesh as that used for the thermal solution. The result is that the elements of the HeatWave solution have assigned source power densities. File sources may be used inb both static and dynamic solutions. In the dynamic mode, a time variation may be associated with the file source. The variation may be defined by a table of [t,f(t)] values using the SourceMod command. Here, f(t) is a multiplication factor. The factor may also be defined by an algebraic expression.

In the new version of HeatWave, the element power distribution is initially copied to a reference mesh variable. A time-dependent displacement vector is defined with the new commands XDisp, YDisp and/or ZDisp. As with the modulation functions, the time-dependent vector components may be defined by either a table of values (e.g., [t, x(t)]) or an algebraic function. The element power densities from the reference mesh are periodically mapped to the computational mesh using the current value of the displacement vector [x(t), y(t), z(t)]. The mapping algorithm conserves energy. Depending on the displacement, some source energy may be mapped outside the source region. As an example, heating of a rotating target could be modeled using a sawtooth function for x(t).

Figure 1 shows the results of a demonstration calculation. A 4 MeV electron beam of radius 5.0 mm with current 23.9 μA impinges on an aluminum target of thickness 10.0 mm. The top section shows the temperature distribution in the center of the target (z = 5.0 mm) with a total irradiation time of 10 s and a fixed source distribution. In contrast, the lower section shows a solution where the source moves the right 20.0 mm during the 10 s interval. Note that the peak temperature drops from 505.1 ºC t0 338.2 ºC.

Electron beam heating demonstration

Figure 1. Electron beam heating demonstration. Top: Stationary source. Bottom: Source moved 20.0 mm to the right over the 10 s calculation period.

HeatWave with the expanded functionality will be included in all future distributions of Xenos. A free upgrade is available to current users.

Footnotes

[1] Contact us : techinfo@fieldp.com.

[2] Field Precision home page: www.fieldp.com.

Magnetic saturation: educated guessing

In finite-element solutions at low magnetic fields, the properties of ferromagnetic materials are not critically important. The quantity used in the solutions is γ= 1/μr. In this case, both 1/800 and 1/8000 are close to zero compared to unity; therefore, the details of the μr(B) curve make little difference in the solution. On the other hand, the variation of relative permeability may be quite important at high field levels where the materials are driven into saturation. Unfortunately, the available data in tabulations of μr(B) are generally the iverse of what we need: there is extensive information at low B but most tabulations end well below the saturation region. The subject of this article is how to guess μr(B) at high field levels from the low-field behavior. Sounds impossible, but its actually easy thanks to two factors:

  • Magnetic saturation occurs smoothly (i.e., there are no abrupt changes).
  • We know the variation at extremely high fields.

The key is choosing a good method to plot the data.

First, some background. The typical measurement setup is a torus of material with a drive coil winding. The quantity H (in A/m) is the magnetic field produced by the coil in the absence of the material. A useful quantity is B0 = μ0*H (in Tesla), the magnetic flux density produced by the coil inside the torus with no material. The quantity B is the total flux density in the torus with the material present. In this case, alignment of atomic currents adds to the field value so that B > B0. In a soft magnetic material (i.e., no permanent magnetization), both B0 and B equal zero when there is no drive current. The alignment of magnetic domains increases as the drive current increases; therefore, B grows faster than B0. The relative magnetic permeability is defined as μr = B/B0. At high values of drive current, all the material domains have been aligned. In this state, the material makes a maximum contribution to the total flux density of Bs (the saturation flux density). This contribution does not change with higher drive current. For high values of B0, the total flux density is approximated by

B ≅ B0 + Bs.  (1)

To illustrate the estimation procedure, we’ll consider the specific example of Magnifer 50 RG, a nickle alloy with a high value of Bs. Figure 1 shows a graph from a data sheet supplied by VDM Metals. The sheet lists the saturation flux density as Bs = 1.55 T. The plot shows B (in mT) versus H (mA/cm) at several frequencies. Because we are interested in the static properties, we’ll consider only curve 1. The data extend to a peak value of H = 200 A/m. At this point, μr > 1000, so that the material is well below saturation. I have an application where the material is driven well into to saturation by applied fields up to H = 5000 A/m, about 400 times the highest known value! Is it possible to make calculations with confidence?

Figure 1. Data on Magnifer 50 RG

Figure 1. Data from VDM Metals on Magnifer 50 RG, B versus H.

The first step is convert the graphics data to a number set. The FP Universal Scale is the ideal tool for this task. After setting the correct log scales, I could record a set of points with simple mouse clicks, including the conversion factors to create a list of B versus B0 in units of Tesla. In this case, the relative magnetic permeability is the ratio μr = B/B0.

The key to estimating the missing values is to create plots of the material behavior at the two extremes: the tabulated values at low B0 and predictions from Eq. 1 at high B0. To ensure the validity of Eq. 1, I picked B0 values corresponding to highly saturated material: 0.1, 0.2 and 0.5 T. The art is picking the right type of plot. Figure 2 shows B versus B0 with log-log scales. With the requirement of a smooth variation, clearly the unknown values must lie close to the dashed red line connecting the data extremes. Accordingly, I used the Universal Scale to find several points along the line. I combined the interpolated values with the low field tabulated values and the high-field predictions to build a data set that spans the complete range of behavior for Magnifer 50. The new data are available on our magnetic materials page.

B versus B0 for Magnifer 50 showing values at high and low field

Figure 2. Plot of B versus B0 for Magnifer 50 showing values at high and low field with a proposed interpolation (dashed red line).

Finally, Figure 3 shows alternate plots to Fig. 2: B0 versus μr and B versus μr. In all cases, the variation over the unknown saturation region is well approximated by a simple straight-line fit.

Alternative plots and interpolations Magnifer 50: B0

Figure 3. Alternative plots and interpolations for Magnifer 50 RG: B0 and B versus relative magnetic permeability.

Footnotes

[1] Contact us : techinfo@fieldp.com.

[2] Field Precision home page: www.fieldp.com.