## research papers

## Obtaining local

vectors from finite-element analysis**John P. Sutter,**

^{a}^{*}Thomas Connolley,^{a}Tim P. Hill,^{a}Houcheng Huang,^{a}Doug W. Sharp^{a}and Michael Drakopoulos^{a}^{a}Diamond Light Source Ltd, Diamond House, Harwell Science and Innovation Campus, Didcot, Oxfordshire OX11 0DE, United Kingdom^{*}Correspondence e-mail: john.sutter@diamond.ac.uk

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 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.

Keywords: finite element; X-ray; diffraction; lattice.

### 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

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 is not infinite. Because the warmer and cooler parts of the crystal undergo different degrees of 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 ), avoided this problem by choosing instead to derive analytical expressions for the local 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 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 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.

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 ). 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.

coefficient of silicon is close to zero (Shah & Straumanis, 1972A 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^{−2}, 0.115, 0.155, 0.35, 0.155, 0.115 and 7.5 × 10^{−2} W mm^{−3}. 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^{−2} K^{−1}). 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^{5} 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 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.

### 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,^{1}

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 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 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 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 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 *i*th 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 **h**′ in the neighborhood of node *I* using (9), substituting *U*_{I}^{avg} for *U*_{N}. In comparing the local 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 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

on the diffraction may be most readily described by the following:|**h**′|/|**h**| − 1 ≃ −δ*d*/*d*, where *d* is the 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),

χ_{0} and χ_{h} are Fourier components of the crystal's electric susceptibility χ(**r**),^{2} 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**_{0}, and the inward normal to the crystal's surface is , then^{3}

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.

### 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 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 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.

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.

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.

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.

### 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

is nearly zero is thus already very useful. Both the variation in 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 *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

^{1}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).

^{2}χ_{0} and χ_{h} are assumed not to be significantly altered if the displacements are small.

^{3}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.

### Acknowledgements

ANSYS, ANSYS Workbench, AUTODYN, CFX, FLUENT and any and all ANSYS Inc. brand, product, service and feature names, logos and slogans are registered trademarks or trademarks of ANSYS Inc. or its subsidiaries in the United States or other countries. All other brand, product, service and feature names or trademarks are the property of their respective owners.

### References

ANSYS Inc. (2007). ANSYS^{®} Workbench™, Release 11.0, Help System, *Theory Reference for ANSYS and ANSYS Workbench.*

Chrzas, J., Khounsary, A. M., Mills, D. M. & Viccaro, P. J. (1990). *Nucl. Instrum. Methods Phys. Res. A*, **291**, 300–304. CrossRef Web of Science

Dieter, G. E. (1988). *Mechanical Metallurgy*, p. 41, SI Metric Edition. London: McGraw-Hill.

Epelboin, Y. (1996). *J. Appl. Cryst.* **29**, 331–340. CrossRef CAS Web of Science IUCr Journals

Freund, A. K., Arthur, J. R. & Zhang, L. (1997). *Proc. SPIE*, **3151**, 216–226. CrossRef CAS

Hart, M. (1990). *Nucl. Instrum. Methods Phys. Res. A*, **297**, 306–311. CrossRef Web of Science

Hoszowska, J., Mocella, V., Zhang, L., Migliore, J.-S., Freund, A. K. & Ferrero, C. (2001). *Nucl. Instrum. Methods Phys. Res. A*, **467**–**468**, 631–634. Web of Science CrossRef CAS

Mocella, V., Ferrero, C., Freund, A. K., Hoszowska, J., Zhang, L. & Epelboin, Y. (2001). *Nucl. Instrum. Methods Phys. Res. A*, **467**–**468**, 414–417. Web of Science CrossRef CAS

Penning, P. & Polder, D. (1961). *Philips Res. Rep.* **16**, 419–440. CAS

Shah, J. S. & Straumanis, M. E. (1972). *Solid State Commun.* **10**, 159–162. CrossRef Web of Science

Suortti, P. & Thomlinson, W. (1988). *Nucl. Instrum. Methods Phys. Res. A*, **269**, 639–648. CrossRef Web of Science

Tajiri, G., Lee, W.-K., Fernandez, P., Mills, D. M., Assoufid, L. & Amirouche, F. (2001). *J. Synchrotron Rad.* **8**, 1140–1148. Web of Science CrossRef CAS IUCr Journals

Takagi, S. (1969). *J. Phys. Soc. Jpn*, **26**, 1239–1253. CrossRef CAS Web of Science

Taupin, D. (1964). *Bull. Soc. Fr. Minér. Crist.* **87**, 469–511. CAS

Zachariasen, W. H. (1945). *Theory of X-ray Diffraction in Crystals*, ch. III. New York: John Wiley and Sons.

Zhang, L., Hoszowska, J., Migliore, J.-S., Mocella, V., Ferrero, C. & Freund, A. K. (2001). *Nucl. Instrum. Methods. Phys. Res. A*, **467**–**468**, 409–413. Web of Science CrossRef CAS

Zhang, L., Lee, W.-K., Wulff, M. & Eybert, L. (2003). *J. Synchrotron Rad.* **10**, 313–319. Web of Science CrossRef CAS IUCr Journals

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.