The APDL is a scripting language that can be used to automate common tasks or build a finite element model (FEM) using parameters (variables). The code has many features such as repeating a command, macros, if-then-else branching, do-loops, and scalar, vector and matrix operations [2]. Its main difference from the ANSYS Workbench environment is the complete control it offers users of the model – both at the pre-processing (data input) and post-processing (results extrapolation) phases – with the possibility of writing the entire program in batch mode, as opposed to the click mode of a graphical user interface (GUI). This feature means that the APDL applications are limited only by the user’s mind and no longer by the GUI functionality.
Particle beams lose energy as heat when interacting with a target material, generating a (thermal) power density inside the latter. Particles, in our case electrons, are not distributed uniformly along the transverse sections of the beam, but according to a Gaussian distribution; consequently, the power density generated is not spatially uniform on the target surface either, but Gaussian as well. Moreover, when the beam is pulsed, the power density generated is not constant over the time, but has a cyclic trend, going from zero to a maximum peak. In our case study, the electron beam was not only pulsed but also had a very high frequency of 100 Hz, which was the same as that of the power density generated.
The ANSYS Workbench GUI does not allow users to simulate this kind of thermal problem, whereas the development of an opportune APDL script does. Moreover, an APDL script allows the parametrization of the physical radiation variables and material characteristics. In our case, the parametrization was done by writing external custom data files containing the information to be used by the APDL script, depending on the type of simulation, while repeated actions (such as the load steps for simulating the power density trend), do-loops, mathematical formulas and matrixes were programmed directly inside the main scripting code. Furthermore, writing the entire program in batch mode allowed optimal management of memory usage by saving only the useful data for each simulation.
Finally, the choice to use the APDL environment was driven by the distribution of the particles of the beam to be simulated in the case study, since they reflected a peculiar spatial and time distribution of the power density in the target material. The APDL script developed may also easily be adapted to suit other applications involving different specific beams, material and setting geometries.

The APDL script was developed during the investigation of the interactions between the OTR screens and electron beams designed for the ELI-NP Gamma Beam System (GBS) facility – an advanced source of up to ≈ 20 MeV gamma rays based on Compton back scattering, i.e. the collision of an intense high-power laser beam and a high-brightness electron beam with a maximum kinetic energy of about 740 MeV. In this case study, the simulated electron beam consists of trains of 32 electron bunches of 250 pC each, separated by 16 ns and distributed along a 0.5 μs RF pulse with a repetition rate of 100 Hz [3].
The analysis was developed starting from previous theoretical analysis results [4] that were used as input for the FE simulations to evaluate the thermal and mechanical performance of two candidate materials for the target (i.e. aluminum and monocrystalline silicon), especially under a high number of thermal cycles, with the objective of choosing the best-performing one.
This methodology takes into account the most demanding operating conditions of the GBS machine: indeed, along the linear accelerator (LINAC), the electron beam is symmetric but in a particular position, at more or less 13 meters from the gun, the electron beam becomes asymmetric with an elliptical shape, where σx is the beam size along the x-axis equal to 47.5x10-6 m and σy is the beam size along the y-axis equal to 109.0x10-6 m.
In this condition, the beam has a high charge density and the impact on the target produces a different behavior in the temperature i.e.a continuous oscillating change of the temperature distribution. In order to evaluate this distribution during the temperature heating and decreasing (thermal cycle) of the OTR, used for the beam diagnostic, several thermal transient analyses were performed. The OTR target’s geometry (30 mm long for each edge and 1 mm thick) was modelled with 3D SOLID70 elements [2] using a dedicated FE code in ANSYS APDL [2]. Moreover, a refinement of the mesh (Fig.1, darker area) was applied close to the electron beam spot where the heat generation is concentrated (minimum size of the mesh elements 6x10-6 m) to obtain a correct distribution of the power density resulting from the beam-target interaction.
The ANSYS APDL code applies the thermal load to those mesh elements that correspond to the portion of the OTR target that significantly interacts with the electron beam. The worst-case scenario was considered to be an elliptical beam spot that corresponded to the more focused beam (Fig.1).

The OTR is the radiation emitted when a fast particle crosses the boundary between two media with a different dielectric constant. This light is emitted with a characteristic angular distribution that depends on the particle energy and the angular momentum [5]. The particle’s electro-magnetic field has a certain configuration in the vacuum in front of the foil, which differs from the one of the field inside the media, because the foil’s dielectric constant is different from vacuum’s. When approaching the foil, the particle’s electro-magnetic field creates a time-dependent polarization at the foil boundary. The change of this polarization causes the radiation with a characteristic intensity and angular distribution.
A typical setup of an OTR measurement is shown in Fig. 2. In most cases, the foil is inserted at a 45° angle to the beam path. The foil is typically made of aluminized silicon with a thickness in the order of tens of micrometers. The light is emitted both in the forward direction and at 90° because the metallic surface acts as a mirror. Depending on the particle energy, the angular distribution is peaked at the angle θ=1⁄γ where γ is the relativistic Lorentz factor.
The light that is emitted at 90° was acquired with appropriate optics, and an image of the foil was recorded with a charge-coupled device (CCD) camera. This is useful to measure some important properties of the beam: indeed, one can directly measure the transverse spot size. Moreover, using a quadrupole scan technique, it is possible to measure the beam’s emittance; while an RF Deflector makes it possible to retrieve the bunch length and, combined with a dipole, to generate the complete 6D characterization of the beam. Additionally, by acquiring the OTR angular distribution (focus on infinity), it is possible to measure the beam’s energy and divergence.

The APDL script is structured in three main parts. The first refers to the definition of the electron beam’s properties and the OTR target’s material and geometry. The second spatially models the power density resulting from the interaction between the electron beam and the target. The third models the power density over time, assuming a pulsed operation mode for the beam described in section 3 and by simulating a high number of thermal cycles on the target material using a specific configuration of the ANSYS solver code (i.e. load steps in Fig. 3). These three parts are described in detail below, although the well-known parts of post-processing a FEM simulation have been omitted.

Definition of beam and OTR target properties

The APDL script supports three kinds of parameters: “variable”, “array” and “table”. This latter is a special type of multi-dimensional array that allows ANSYS to calculate the values between the array elements through linear interpolation.
Every physical property, from the electron beam’s characteristics to the geometric dimensions and the material properties of the OTR target, becomes a parameter of the simulation and is modelled as an ANSYS variable. These variables are defined in the pre-processing section of the APDL script, as is usual in FEM analysis. The parametrization of the physical problem allows one to easily tailor the simulations to specific user cases and to execute a sensitivity analysis to tune the beam and/or the equipment properties to achieve the best system performance.

Spatial modelling of the power density due to the beam target interaction

As previously stated, the electrons along the transverse sections of the beam have a Gaussian distribution; consequently, resulting power density also has a Gaussian distribution along the target surface. We assumed the beam particles to be symmetric along the beam’s longitudinal axis (z): this is more than reasonable when the range of the beam’s particles in the target material is larger than the thickness of the target, as in our case. On assuming the coordinates of the target surface (x,y) to be geometrical, the power density Q can be expressed analytically in the following way (1):

The most efficient way to model Q(x,y) is to define an ANSYS 2-dimensional table parameter, hereafter called Q_TABLE. Each (i,j) cell of the table matrix corresponds to the power density Q calculated in the x-coordinate of the centroid of the i-mesh element, and in the y-coordinate of the centroid of the j-mesh element. The dimensions of the table matrix depend on the number of mesh that discretized the FEM problem. To reduce the computational time and resources, the APDL script circumscribes the number of “powered” mesh to a set, the so-called “hotspot”, around the beam-target interaction area. In this way, it avoids generating huge matrixes with many 0-value elements corresponding to the mesh elements geometrically far from the beam-target interaction area and where the power density is negligible.

Simulation of a high number of thermal cycles: do-loops and load steps

The simulation of a high number of thermal cycles required the writing of a specific do-loop with two different load steps: the first one regarding the heating of the portion of the material hit by the radiation, and the second one representing the consequent temperature decrease. Simulations were done in transient mode. An extract of the APDL code used in the solution part of the script follows (Fig. 3). The main challenge was the time discretization in relation to the physical characteristics of a generic radiation (the T, DT, PULSE_ DISTANCE, PULSE_LENGTH, BUNCH_PULSE, BUNCH_SPACING variables). For this reason, the time discretization of the heating phase was done using a custom time-step for the entire load step solution (AUTOTS,OFF and DELTIM,ON). On the other hand, for the second load step, the time discretization did not have particular requirements, and was managed directly by the solver (AUTOTS,ON). The Gaussian power density distribution is a relevant aspect of the script and it is represented by means of the Q_TABLE parameter described above. For each load step of the solution, the solver increases the global number of *DO cycles, involving several Gigabyte of storage memory.
By setting the OUTRES result to the LAST parameter (Fig.3), the solver saves only the FE results for the latest sub-step of the transient solution, saving memory and including only the useful data results.

In order to evaluate the variation of temperature in the target during the interaction with the beam, the thermal transient analysis was conducted to simulate the power density described above. Fig. 4 compares the “hotspot” area between the analytical power density based on formula (1) (green graph) and the power density as simulated in ANSYS environment through the Q TABLE parameter, as explained in the previous section (blue graph). The dots of the blue graph represent the simulated power density corresponding to the mesh centroids, whereas the lines represent the linear interpolation applied by the Q_TABLE between the power density values in the centroids. Fig. 4 highlights the ANSYS environment’s good approximation of the analytical power density through the use of the Q TABLE parameter and the mesh settings illustrated above. The main properties of aluminum (Al) and monocrystalline silicon (Si) assumed in the FE simulations are reported in Table 1. These properties are necessary to eventually evaluate the fatigue damage of the material. In fact, once the temperature evolution in the OTR screen over time had been calculated, the equivalent Von Mises stress state [6] could be obtained for the steady temperature reached during the heating (ANSYS first load step) and the temperature decrease (ANSYS second load step).

Table 1 - aluminum and monocrystalline silicon properties [7]

The first thermal boundary condition applied to the OTR screen was the initial temperature of 295.15 K for all nodes of the mesh, corresponding to Gamma Beam System/- (GBS) room temperature. An additional boundary condition was the fixed temperature of 295.15 K along the OTR target edges that were in contact with the frame support and the screws (Fig. 1). Fig. 5 represents the spatial thermal distribution (x, y) when the beam hotspot reaches the maximum temperature for both materials in the first thermal cycle. The Si has a higher and a more spatially concentrated ΔT respect to the Al. The Al has a higher specific heat capacity and thus a lower maximum temperature than the Si for the same amount of deposited beam power.
In order to evaluate the number of cycles needed to reach constant temperature, a dedicated transient analysis was set up through the ANSYS APDL code (600 thermal cycles). Table 2 shows the steady maximum temperature for the aluminum that was reached after 80 cycles (0.80 s) and is equal to 345.3 K, whereas for the silicon it was reached after 92 cycles (0.92 s) and is equal to 358.9 K. The different steady temperature of the two materials is due to their different specific heat capacity as well.
The maximum and the minimum stress induced by thermal deformation – and calculated with this code – were then used to evaluate fatigue and related lifetime. Al reaches a Goodman alternate stress of 30.20 MPa whereas the Si reaches 7.74 MPa. After one hour of exposition (36,000 cycles) – taking into account the related Wöhler curves for each material – the Al reached a cumulative fatigue damage of 0.59, while the Si did not reach the fatigue life limit at all.

Table 1: aluminum and monocrystalline silicon properties [7]