Obtaining local reciprocal lattice vectors from finite-element analysis

Finite-element analysis is frequently used by engineers at synchrotron beamlines to calculate the elastic deformation of a single crystal undergoing mechanical bending or thermal load. ANSYS^® Workbench™ software is widely used for such simulations. However, although ANSYS^® Workbench™ software provides useful information on the displacements, strains and stresses within the crystal, it does not yield the local reciprocal lattice vectors that would be required for X-ray diffraction calculations. To bridge this gap, a method based on the shape functions and interpolation procedures of the software itself has been developed. An application to the double-crystal bent Laue monochromator being designed for the I12 (JEEP) wiggler beamline at the Diamond Light Source is presented.

1. Introduction

Nearly all synchrotron beamlines have monochromators used to select a wavelength that is desirable for the experiment being performed. Those beamlines built for hard X-rays normally use Bragg-diffracting perfect crystals, usually of silicon but sometimes of diamond or germanium. However, the high radiation flux produced by synchrotrons, particularly at beamlines that use insertion devices, adds considerable heat to the first monochromator crystal struck by the beam. This crystal is thus cooled with either water or liquid nitrogen to keep it from overheating, but the temperature of the crystal remains higher within the volume that absorbs the beam than in the part of the crystal close to the coolant because the crystal's thermal conductivity is not infinite. Because the warmer and cooler parts of the crystal undergo different degrees of thermal expansion, the crystal is distorted, with serious effects on the efficiency of the diffraction. This problem has attracted much interest. The finite-element method has proven to be one of the most effective tools for predicting the distortion of a crystal caused by a given heat load. This allows the designer to choose the best among various proposed cooling schemes. Finite-element analysis also allows mechanical benders to be simulated, thus ensuring that a crystal is bent to a desired shape. In both cases, finite-element analysis yields the displacements, strains and stresses at the nodes of the small volumes into which the crystal is divided for the calculation.

However, it is not trivial to determine the resulting effects on a crystal's ability to diffract X-rays, even if finite-element analysis has provided the displacements and strains. This is because the local deformations of the Bragg-diffracting atomic planes, or in other words the spatial variation of the reciprocal lattice vector, cannot be derived from these quantities alone. Some, for example Hart (1990), avoided this problem by choosing instead to derive analytical expressions for the local Bragg angle, but these are not easily extended to other designs. Others who did simulate the rocking curves of crystals under heat load using displacements determined by finite-element analysis (Chrzas et al., 1990; Freund et al., 1997; Zhang et al., 2001; Mocella et al., 2001; Hoszowska et al., 2001; Tajiri et al., 2001; Zhang et al., 2003) derived slope errors from the given displacements, but only within the diffraction plane. Moreover, they did not mention the lattice spacing variation, which also affects the local Bragg angles. For those designs that are run near room temperature where the thermal expansion of silicon is significant, this is a serious issue whose neglect may partly explain the discrepancies between theory and experiment that some of these papers describe. This paper will complement these previous works by providing a systematic method for deriving all components of the spatial variation of the reciprocal lattice vector from the displacements determined by finite-element analysis.

We will apply our calculations to the I12 (JEEP) wiggler beamline of Diamond, which will use high-energy (50–150 keV) X-rays for imaging and diffraction experiments. The low Bragg angles for X-rays at this energy make Bragg-case monochromators unfeasible; therefore, a double-crystal Laue-case monochromator (Fig. 1) is to be used instead. Both crystals will be bent in order to increase their bandpass as shown by Suortti & Thomlinson (1988). To keep the exit beam fixed in position, the distance between the two crystals will need to be adjusted as the selected energy is varied. In addition, however, the radii of curvature of the crystals must be made adjustable to keep the wiggler source on their Rowland circles, and so a mechanical bender must be designed along with the cooling system.

Figure 1 Schematic of the double-crystal bent Laue monochromator being designed for the I12 (JEEP) beamline at the Diamond Light Source. θB is the Bragg angle.

The Laue geometry spreads the absorption of the incident beam throughout the crystal volume, rather than concentrating it in a thin layer near the surface as the Bragg geometry does. Therefore, the thermal gradient may be gentler, but on the other hand the crystal interacts with the beam throughout its whole thickness. The main part of the heat load is imposed by the low-energy background of the wiggler; therefore, an effective upstream filter needs to be designed. Moreover, to increase the field of view for imaging experiments, the size of the beam, 45 mm × 13.5 mm, is to be larger than that used in any other high-energy beamline using Laue-case monochromators, such as at the ESRF. To reduce the thermal distortion as much as possible, the crystals are to be cooled to about 120 K, where the thermal expansion coefficient of silicon is close to zero (Shah & Straumanis, 1972). Finally, as mentioned above, adjustable benders will be required to control the radii of curvature of the crystals. Because the cooling system will need to cope with the greater heat load of the larger beam, and because the benders must deform the crystal as nearly cylindrically as possible, a simulation of the beam diffraction of this monochromator is essential for designing the filter, the cooling system and the monochromator itself.

A diagram of the first crystal of the I12 monochromator, a single perfect 3.0 mm-thick silicon crystal, is shown in Fig. 2. It will use the (111) Bragg reflection at an asymmetry angle of −44°, thus diffracting at grazing exit. For the following simple trial the crystal was assumed to be influenced solely by the heat load and the cooling system. Mechanical bending was not taken into account; the effect of this on the diffraction will be reserved for later. The heat load is determined by the construction of the I12 wiggler (45 dipole magnets of alternating polarity with 4.2 T field strength and 48 mm periodicity), by the storage-ring current and by the beamline windows and filters. If the ring current is at its design maximum of 500 mA, then the total power emitted by the wiggler is 56 kW, of which the front-end window takes a section of 1.0 mrad horizontal × 0.3 mrad vertical divergence with a total power of 9 kW. Two filters, the first of diamond and the second of 8 mm SiC, remove the low-energy background, thus reducing the transmitted power to 1.6 kW, of which the first crystal of the monochromator absorbs 443 W. In the finite-element analysis of the crystal, this incident power is not assumed to be uniform across the crystal, but instead the volume of the crystal penetrated by the incident beam is divided into the seven parallelepipeds shown in Fig. 2. Each of these is 2.73 mm wide. From left to right the heat assumed to be generated within each parallelepiped was 7.5 × 10⁻², 0.115, 0.155, 0.35, 0.155, 0.115 and 7.5 × 10⁻² W mm⁻³. The cooling system is treated in a similarly simple way. The shaded surfaces of the crystal in Fig. 2, and the corresponding surfaces at the bottom that are not visible in the figure, are assumed to be in contact with the liquid nitrogen flow at 77 K (−196°C), with a heat transfer rate of 0.005 W mm⁻² K⁻¹). ANSYS^® Workbench™ (version 11.0) was then used to perform the finite-element analysis. It was assumed for simplicity that the crystal was elastically isotropic with Young's modulus 1.558 × 10⁵ MPa and Poisson's ratio 0.2152. This yielded the temperature distribution (Fig. 3) and the displacements (Figs. 4, 5 and 6) derived within the volume of the crystal that the beam is to traverse. It should be noted that a proper accounting of the elastic anisotropy will alter the values of the displacements, but it will not affect the method discussed in the next section for determining the local reciprocal lattice vectors at the element nodes. The finite-element software is capable of including the elastic anisotropy, although the computation is then of course more complex; this is planned for future refinements.

Figure 2 Diagram of the first crystal of the I12 monochromator, assuming that the crystal is indirectly cooled by a liquid-nitrogen current close to its sides. The incident beam actually strikes the center of the crystal, but only the back half of the crystal is shown here. The region in which the incident beam strikes the crystal is shown by the volume markings on the forward end. The slanted planes drawn within the crystal volume indicate the reflecting (111) atomic planes.

Figure 3 Temperature distribution within the region of the first crystal penetrated by the incident beam under preliminary assumptions about the low-energy filtering and the cooling system. The temperatures are given in degrees C.

Figure 4 Displacements along the X direction of points within the region of the crystal penetrated by the incident beam.

Figure 5 Displacements along the Y direction of points within the region of the crystal penetrated by the incident beam.

Figure 6 Displacements along the Z direction of points within the region of the crystal penetrated by the incident beam.

2. Finite-element method

2.1. Treatment of deformations

The treatment of this and the next section is based on the information set out in the Theory Reference for ANSYS and ANSYS Workbench (ANSYS Inc., 2007). Defining a Cartesian coordinate system in which the position R of a point in the undistorted crystal is given by , we then obtain the resulting displacements u = of the nodes of the finite elements. ANSYS^® Workbench™ software additionally provides the normal strains ∊_ii and the engineering shear strains γ_ij, which are defined as follows if the deformation is small,¹

ANSYS^® Workbench™ also calculates the corresponding stresses within the crystal.

We wish to calculate the effect of the crystal deformation on a Bragg reflection from some family of atomic planes defined in the undistorted crystal by

where m is an integer and h is a reciprocal lattice vector. The deformation u transforms each point R in the undistorted crystal to a new point r = R + u. Thus, after the deformation, the atomic planes in (7) are transformed into new surfaces described by

If the deformation is small, u may be treated as a function of r rather than R, so that a local reciprocal lattice vector h′ = h − ∇_r(u·h) may be defined. The Cartesian components of h′ may be rewritten in matrix form as follows,

where I is the identity matrix and U_N is the matrix of the derivatives of the displacement. The subscript N indicates that the quantity is evaluated in the neighborhood of one of the finite-element nodes. Thus it is clear that a complete knowledge of all derivatives of the displacement is required to calculate the local reciprocal lattice vector. Unfortunately, ANSYS^® Workbench™ software is not able to provide this in its present form. The strains are not sufficient, as they provide only the sums of the off-diagonal terms in U_N. The following treatment corrects this shortcoming, allowing the local reciprocal lattice vectors to be determined in the neighborhood of each node in the finite-element analysis.

2.2. Derivation of the displacement derivatives

ANSYS^® Workbench™ software uses shape functions to interpolate to the displacement at any point within each element from the displacements of the nodes of that element. Each node i = 1,…, n, where n is the total number of nodes of the element, has an associated shape function N_i given in terms of the element's `natural coordinates' r, s and t. These coordinates are a linear (but not necessarily orthogonal) set defined such that the surfaces r = ±1, s = ±1 and t = ±1 form the element's boundaries. From the above we obtain the components of the displacement at some arbitrary (r, s, t),

where (u_i, v_i, w_i) is the displacement at node i.

The matrix of displacement derivatives U at some point (r, s, t) within a given element can thus be found by differentiating (10), (11) and (12),

In order to obtain the derivatives of the shape functions with respect to the global Cartesian coordinates (x, y, z), as opposed to the element's natural coordinates (r, s, t), we need further information. This is provided by the assumption that all elements are isoparametric; that is, the global Cartesian coordinates of point (r, s, t) of any element are interpolated by the shape functions in just the same way as the displacement at the same point,

where (x_i, y_i, z_i) are the global Cartesian coordinates of the element's ith node.

Now, using the chain rule, we find

where J is the Jacobian matrix given by

The right-hand side can be found by differentiating (14), (15) and (16). With this knowledge we can obtain the derivatives of the displacement with respect to the global Cartesian coordinates,

2.3. Details of the finite-element calculation

ANSYS^® Workbench™ software provides numerous types of elements of different shapes and with different numbers of nodes. For our calculations we used the 20-node brick SOLID95. Eight of the nodes lie on the corners [(r, s, t) = (±1, ±1, ±1)], and the remaining twelve lie at the midpoints of the edges [(r, s, t) = (0, ±1, ±1), (±1, 0, ±1), (±1, ±1, 0)].

The shape functions of this element are defined in ANSYS^® Workbench™ software as follows.

For corner node i at (r, s, t) = (r_i,s_i,t_i), where r_i,s_i,t_i = ±1,

For edge node i at (r, s, t) = (0,s_i,t_i), where s_i,t_i = ±1,

For edge node i at (r, s, t) = (r_i,0,t_i), where r_i,t_i = ±1,

For edge node i at (r, s, t) = (r_i,s_i,0), where r_i,s_i = ±1,

ANSYS^® Workbench™ software calculates the strains at the Gaussian integration points of each element and then extrapolates these values to the element's nodes, as will be explained in more detail in the next section. We follow the same procedure in calculating the derivatives of the displacement. Each type of element has a defined number of integration points situated at defined natural coordinates. SOLID95 contains 14 integration points, of which eight lie close to the corners and six near the centers of the faces. The natural coordinates of these points are as follows.

Corner points: (r, s, t) = (±A, ±A, ±A), where A = 0.75868 69106 39328.

Face center points: (r, s, t) = (±B, 0, 0), (0, ±B, 0) and (0, 0, ±B), where B = 0.79582 24257 54222.

3. Calculation of displacement derivatives from displacements

The first task is to calculate the displacement derivatives at the integration points of each element. ANSYS^® Workbench™ software provides as input for our calculations the global Cartesian coordinates and the displacements of the nodes of each element. Using the appropriate shape functions given in (20)–(23), we use (19) to calculate the matrix of the displacement derivatives U(r_int, s_int, t_int) with respect to the global Cartesian coordinates at each integration point (r_int, s_int, t_int).

The next step is to extrapolate all of the displacement derivatives to the nodes. For this, ANSYS^® Workbench™ software uses only the eight corner integration points (r, s, t) = (±A, ±A, ±A), so that the values at each node (r_i, s_i, t_i) are given by

So far we have treated only a single element. However, in general each node will be shared by more than one element. Frequently, two different elements sharing the same node will not give the same U for that node after the extrapolation is carried out. Therefore, for each node I = 1,…, N, where N is now the total number of nodes in the entire crystal, we average all the values of U calculated at that node over all the elements that include it. We call the result U_I^avg. It is this that will be used in all subsequent calculations.

3.1. Calculation of lattice distortion

From the above we may now obtain the local reciprocal lattice vector h′ in the neighborhood of node I using (9), substituting U_I^avg for U_N. In comparing the local reciprocal lattice vectors of the deformed crystal with those of the undeformed crystal, it is useful to define the `diffraction plane' as the plane which contains the incident beam, the diffracted beam and the undistorted reciprocal lattice vector h. This must be distinguished from the `reflecting atomic plane' that is perpendicular to h and is responsible for the Bragg reflection. The following additional unit vectors may then be defined in relation to the diffraction plane,

= the normal to the diffraction plane,

= h ×

The effect of the local lattice distortion on the diffraction may be most readily described by the following:

|h′|/|h| − 1 ≃ −δd/d, where d is the interplanar spacing of the reflecting atomic planes of the crystal;

δθ_in, the rotation of the reflecting atomic planes within the diffraction plane. This is the angle between h and , the component of h′ that lies in the diffraction plane. It is positive if > 0;

δθ_out, the rotation of the reflecting atomic planes out of the diffraction plane. This is given by . It is positive if > 0.

Another useful parameter for estimating the degree of deviation from the Bragg condition at a given node is given by Zachariasen (1945),

χ₀ and χ_h are Fourier components of the crystal's electric susceptibility χ(r),² which of course has the same spatial periodicity as the lattice, so that

K is a polarization factor, equal to 1 if the incident beam is polarized perpendicular to the diffraction plane and to cos2θ_B (θ_B is the Bragg angle) if the beam's polarization lies within the diffraction plane. If the incident beam's wavevector is k₀, and the inward normal to the crystal's surface is , then³

In general, at the peak of the reflection in an unbent undistorted crystal, |y| ≤ 1 (Fig. 7). This defines the Darwin width of a Bragg reflection from a perfect crystal. Thus, if the local values of y at two different points in the crystal differ by 1, then the incidence angle (or wavelength) at which the Bragg condition is fulfilled at these points will differ by half the reflection's Darwin width. From the discussion above, it is clear that a bent crystal will diffract rays over a considerably larger range of y. For example, a flat crystal bent to a radius of 60 m, approximately the correct value for the I12 beamline, will diffract 50 keV rays over a total range Δy = 38.6 and 150 keV rays over a range Δy = 115.6, whereas an unbent flat crystal will diffract these rays only over a range Δy ≃ 2.

Figure 7 General dependence of the rocking curve of an unbent perfect crystal on the parameter y defined by equation (25). Note that the maximum reflection lies within the region |y| ≤ 1. y is directly dependent on the local reciprocal lattice vector as shown by equation (28), from which the rocking curve's dependence on δd, δθin or δθout may be determined.

4. Results

The method of the previous section was applied to a silicon crystal struck by a beam of X-ray photons from the I12 wiggler while being cooled by a flow of liquid nitrogen. The 50 keV photons in this beam are to be reflected by the crystal's (111) atomic planes, which lie at an asymmetry angle of −44° from the crystal's surface as shown in Fig. 2. These photons are thus diffracted at a Bragg angle of 2.266° in grazing exit. As in this figure, only the back half (X ≥ 27.5 mm) is displayed, since the front half (X ≤ 27.5 mm) is a mirror image of the back. A mesh of 0.5 mm was found to provide sufficient detail without requiring excessive computation time. Samples of the local reciprocal lattice vectors were taken across three planes parallel to the crystal surface: the central plane Z = 0 mm halfway between the entrance and exit surfaces of the crystal, the plane Z = +1 mm that lies 0.5 mm below the crystal's entrance surface, and the plane Z = −1 mm that lies 0.5 mm above the crystal's exit surface.

The local variations in δd/d displayed in Fig. 8 show no significant dependence on the depth inside the crystal, as would be expected from the temperature distribution shown in Fig. 3. On the other hand, the dependence of δd/d on the horizontal coordinates is significant, the total range of variation extending to eight parts per million even with the assumption made here that the crystal is not deliberately bent.

Figure 8 The relative deviation δd/d of the (111) interplanar spacing d in parts per million from the undistorted value at room temperature (see color scales at the right-hand corner of each figure) in the plane Z = +1 (top), Z = 0 (middle) and Z = −1 (bottom). The overall negative sign of δd/d results from the simulation of crystal cooling from room temperature by liquid nitrogen, which of course contracts the crystal. See the text of §4 for definitions of the Z planes.

Fig. 9 shows that δθ_in varies significantly not only horizontally within each Z plane but also with depth. The Y coordinate of the center of the region of maximum deformation within each plane lies at +4 mm for Z = +1 mm, at 0 mm for Z = 0 mm and at −4 mm for Z = −1 mm.

Figure 9 The deformation δθin in µrad (see color scales at the right-hand corner of each figure) in the plane Z = +1 (top), Z = 0 (middle) and Z = −1 (bottom). See the text of §4 for definitions of the Z planes.

Fig. 10 shows that |δθ_out| does not exceed 1 µrad. Since the diffraction from the (111) atomic planes is insensitive to warping of the planes in this direction, the effect of δθ_out is negligible.

Figure 10 The deformation δθout in µrad (see color scales at the right-hand corner of each figure) in the plane Z = +1 (top), Z = 0 (middle) and Z = −1 (bottom). See the text of §4 for definitions of the Z planes.

The effect of δd/d, δθ_in and δθ_out on the diffraction is summarized in Fig. 11. The variation in the Bragg parameter y extends to 1.418, or 0.709 times the Darwin width of the (111) rocking curve in the perfect crystal. Since δd/d and δθ_out vary little with depth, most of the variation of y with depth is due to the warping of the atomic planes given by δθ_in. We note that the variation in y caused by the heat load and the cooling system is much less than the range Δy over which a crystal bent to a 60 m radius will diffract X-rays. This encourages us to continue to refine the monochromator design.

Figure 11 The deviation δy of the Bragg parameter y from the center of the rocking curve y = 0 (see color scales at the right-hand corner of each figure) in the plane Z = +1 (top), Z = 0 (middle) and Z = −1 (bottom). Each unit change of δy corresponds to half the Darwin width of the Bragg reflection from an unbent perfect crystal. See the text of §4 for definitions of the Z planes.

5. Conclusions

The trial finite-element analysis reported here shows that the distortion of a silicon crystal caused by the thermal load of the incident beam from the Diamond I12 wiggler covers a significant fraction of the rocking curve of the asymmetric (111) Laue-case reflection in a perfect crystal, but only a small fraction of the total range of X-rays diffracted by the crystal when bent. The assumption that the crystal is cooled by liquid nitrogen close to that temperature of 120 K where the thermal expansion is nearly zero is thus already very useful. Both the variation in interplanar spacing and the warping of the diffracting atomic planes influence the deviation from the ideal Bragg condition. Further modeling calculations are now being carried out to find ways to reduce the thermal distortion and also to design a mechanical bender for the crystal.

The distribution of local reciprocal lattice vectors obtained by the method of this paper can be used as input for an X-ray diffraction program such as TRANSQ (Epelboin, 1996), or into some similar program based on the Takagi–Taupin (Takagi, 1969; Taupin, 1964) or Penning–Polder theories (Penning & Polder, 1961) developed to calculate X-ray diffraction in distorted crystals. It may also be used in a ray-tracing program based on kinematic diffraction theory, which is much simpler and appears particularly suitable for the diffraction of high-energy photons from a highly asymmetric bent crystal.

Footnotes

¹The engineering shear strains γ_ij are to be distinguished from ∊_ij, the off-diagonal components of the strain tensor. For , ∊_ij = γ_ij/2 (Dieter, 1988).

²χ₀ and χ_h are assumed not to be significantly altered if the displacements are small.

³Strictly speaking, we should use h′ in the formula for b, but the effect of using h instead is small. The effect on y is contained chiefly in α, where h′ must therefore be used.