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.