- Email: [email protected]
An iterative approach to curved ray reconstruction of birefringent fields
Optics and Lasers in Engineering 22 (1995) 347-312 Copyright Q 1995 Elsevier Science Limited Printed in Northern Ireland. All rights reserved 0143-8166/95/$9.50 ELSEVIER
014%8166(94)ooo34-4
An Iterative Approach to Curved Ray Reconstruction of Birefringent Fields Allan T. Dolovich, University
of Saskatchewan,
Zhaoqing
Wang,
Department of Mechanical Engineering, Canada. S7N OWO
Saskatoon,
G. M. L. Gladwell University
of Waterloo,
Solid Mechanics Division, Waterloo,
Canada, N2L 3Gl
ABSTRACT The transmission of polarized light through non-homogeneous birefringent fields is examined, and an iterative approach for incorporating ray curvature into the analysis of tomographic data is proposed. The method extends strain-gradient theory and applies iterative algorithms which have been used by other researchers to reconstruct fields which are non-homogeneous but optically isotropic. Optical anisotropy is incorporated by using a light propagation model which includes both bending and separation of initially coincident light ray pairs. To illustrate the method, data are generated computationally and are analysed to determine the bending moment in a prismatic beam and the forces acting on a diametrically loaded disc.
INTRODUCTION Transmission methods of optical stress analysis are often based on the phenomenon of artificial birefringence (also known as temporary double refraction). ‘J This phenomenon is exhibited by certain transparent, polymeric materials which may be used to model stressed bodies of different geometry. If an initially isotropic, homogeneous model is subjected to an applied load, the resulting stress/strain field causes the material to become optically anisotropic. Consequently, a ray of linearly polarized light (that is, a light wave with a single, constant 347
348
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
plane of vibration) entering the field is split into two rays which are linearly polarized at right angles to each other. Each ray propagates with a different velocity and therefore a different refractive index. The refractive indices can be mathematically related to the stress state within the medium. Thus, measured changes in phase and direction of probing light rays can be used to extract information about the field. If sufficient information is obtained, the unknown stress field can be reconstructed. Classical transmission photoelasticity is perhaps the most wellknown method which makes use of this phenomenon. It relies, however, upon a substantially simplified view of light propagation within a non-homogeneous birefringent medium. For example, probing light rays are assumed to travel in straight lines. This assumption implies that originally coincident slow and fast light rays remain collinear, and their vibrations are superimposed at the data plane. If the stress field is assumed to be two-dimensional, the straight ray model also implies that: the axis of vibration for each ray remains constant in orientation through the thickness of the body, the components of stress influencing a given ray remain constant along that ray, and a point in the photoelastic body directly projects onto its corresponding point in the data plane (i.e. in the plane of the observed fringe pattern). These simplifications lead to an elementary description of the relationship between photoelastic fringe data N and principal stress differences uz. The dimensionless interferometric data N are considered to be (-rl a direct measure of the phase difference between collinear fast and slow rays. Moreover, a linearized stress-birefringence relationship or sfressoptic law is used, and it is assumed that the principal directions of stress and birefringence coincide. In the resulting mathematical model, fringe data N are simply proportional to principal stress differences g1 - c2, and the inverse problem of determining g1 - a2 from data N is mathematically trivial. These simplifications also facilitate the interpretation of isoclinics and the separation of principal stresses; at each point in the field, the isoclinic parameter is regarded as a direct measure of the principal stress orientation at that point. This simplified theory has been extremely useful within the context of its applications. However, a number of methods now demand a more complete description of the interaction between light and birefringent media. Examples include isodyne methods4*5 and procedures for tomographically reconstructing optically anisotropic fields.6 Applications include determination of residual stress in glass and plastics7 and the assessment of local effects.8 In particular, the reconstruction of three-dimensional birefringent
Reconstruction of birefringent fields
349
fields from tomographic data poses a tremendous challenge. Here, we use the term ‘tomography’ in a very general way to include any approach where an unknown field is reconstructed from data generated by field-induced alterations of probing radiation? The tomography of scalar, optically isotropic fields has been studied in detail and is well understood.lO*” In contrast, the tomography of stress-optic fields is much more complex, and a number of factors must be incorporated into the physical models used for data analysis. Images produced in the data plane represent extremely complicated integrated effects where both the magnitude and direction of the secondary principal stresses change continuously along any given light path. The method of integrated photoelasticity’* has been used to study these effects under the assumption of rectilinear light propagation, with axisymmetric problems of particular interest.1S’5 There is experimental evidence, however, that the ramifications of non-rectilinear light propagation cannot be ignored. In particular, it has been shown that the essential character of integrated photoelastic data depends on the object plane of the imaging system.16 Pindera and Hecker” have developed a strain-gradient theory which relates the derivatives of principal stresses to the bending of light in photoelastic specimens. The field-induced deflection of polarized laser beams through a specimen are measured, and the theory is used primarily to determine: the gradients of the sums and differences of principal stresses, stress-optic coefficients over a wide spectral range, and the errors which may result in classical or integrated methods which methods ignore the strain-gradient effect.‘* In addition, strain-gradient have shown some promise as stress reconstruction techniques. For example, Aben” has proposed a method for determining the complete stress state in cylinders by supplementing tomographic photoelastic data with ray deflection measurements. The contention here, however, is that strain-gradient theory should be incorporated into the theoretical fabric of tomographical methods which currently ignore ray bending. We believe that a reasonable objective is to work towards a post-processing approach for analysing integrated birefringence data in such a way that refraction effects are taken into account. When refraction effects are included, even the so-called direct problem of predicting data produced for a known stress field demands a higher level of sophistication regarding the mathematics and physics used to model the complex phenomena involved. Qualitative effects due to ray bending include the physical separation of initially collinear rays, the rotation of vibration axes, and the distortion of propagating
350
Allan T. Dolouich, Z. Wang, G. M. L. Gladwell
wavefronts. The effects of ray separation and wavefront distortion necessitate that non-classical parameters, such as the characteristics of the light source and imaging system, be included in the analysis. Specifically, the interaction between non-parallel rays of light depends on the plane of observation, and therefore the tomographic data set obtained for any particular case depends on the object plane of the imaging system. Also, because associated light components (originating at the same point of incidence) diverge and propagate along different paths, the data, in general, cannot be an exact measure of the interference between collinear rays, and coherence of the light source is a major issue. These complications are compounded by the inherent mathematical nonlinearity associated with integrated effects along curvilinear paths. Thus, the direct problem presents a considerable challenge. Accordingly, the inverse problem of reconstructing stress fields from tomographic data appears to be completely unmanageable. Since ray paths are functions of the unknown field, they are themselves unknown quantities, and the complexity of the direct problem would seem to preclude direct data analysis to simultaneously determine ray paths and field solutions. The spatial relation between a point in the data plane and its associated stress points in the body is unknown without some information regarding the rays; that is, the displacement of interferometric fringes obscures the field-data correspondence. Clearly then, at least an approximation to the ray paths is a prerequisite for inversion incorporating the effects due to refraction. The solution proposed here is an iterative approach. If, for a given problem, the straight-ray reconstruction is regarded as but a first approximation to the stress field, it may be used in conjunction with strain-gradient theory to determine an estimate of light ray paths, and a curved-ray inversion may be performed to yield (what is hoped to be) an improved approximation to the field. Then, in turn, this approximation may be used to redefine the rays and these can form the basis for the next inversion and so on until a suitable termination criterion is satisfied. This recursive process is illustrated in Fig. 1. This approach is conceptually identical to algorithms which have been used in ultrasound tomography,2’ geophysical imaging,” and holographic interferometry. 23 Within the context of holographic interferometry, it has been described as Curved Ray Algebraic Inversion (CRAI). It has been successfully applied to the reconstruction of optically isotropic fields, and its development has been followed by the introduction of other refraction-correction algorithms (see, for example, Refs 11 and 24). The objective here is to take a first step towards
Reconstruction of birefringent fields
351
Straight-Ray Reid Appmxlmation
Ray Appmximatbn
Fold Approximation
Fig.
1. Iterative scheme for curved-ray data analysis.
extending the CRAI approach to problems involving artificially birefringent media. It should be emphasized, however, that we are addressing but one issue within a larger scheme. Refraction-related processes are accompanied by other important phenomena which have been identified in the recent literature. For example, one of the fundamental assumptions underlying the Ramachandran-Ramaseshan model of artificial birefringence,25 and therefore the linearized stress-optic law, is that the optical axes of vibration coincide with the secondary principal directions of stress. The optical axes corresponding to the splitting of light rays in an anisotropic refractive index field can be predicted from first principles and Maxwell’s equations for electromagnetism.’ However, the co-directionality of these axes with those of secondary principal stress is assumed on the basis of a speculative, phenomenological approach.25 Thus, the classical stress-birefringence relations may not be appropriate in cases where the angular separation between the two sets of axes is significant. Advanced physical models of artificial birefringence should also include wavelength of the light source, A, as a major parameter. The of optical axes, and stress-optic refractive indices, orientation coefficients are all functions of A, and therefore the corresponding effects of dispersion, dispersion of optical axes, and dispersion of birefringence all contribute to the real response of photoelastic bodies.26-28 Furthermore, the effects of creep and stress relaxation in polymers must be considered. Linear viscoelasticity29 is often a reasonable assumption, but the limits of linearity must be established in terms of stress, strain, and birefringence. In summary, a complete description
352
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
of light propagation through a non-homogeneous, birefringent medium involves a complex interaction between these and other phenomena. By performing computer simulations, however, we can gain some experience with individual aspects of the system. In this paper, we present a data generation model which emphasizes the effects due to refraction. We then present an iterative approach to curved-ray analysis, and we demonstrate the approach by using two elementary examples: a beam in pure bending and a diametrically loaded disc. In each example, we use an analytical reference solution for the stress field together with our data generation model to numerically produce an interferometric data set. We then compare straight-ray analyses of these data to curved-ray solutions using the iterative method. In using computer simulation to address the inverse problem of curved-ray inversion, we take a preliminary step towards establishing its feasibility for practical applications. Mathematical inversion of the tomographic process is not the only aspect of curved-ray tomography, but it is a necessary condition for the successful analysis of real data obtained in physical experiments.
DATA
GENERATION
MODEL
To illustrate the concept of curved ray inversion, we may use a relatively simple example. Consider the propagation of polarized light within a plane of symmetry in a two-dimensional birefringent field, as shown in Fig. 2. The light source (not shown in Fig. 2) is a laser, followed by a pinhole placed at the focal point of a collimating lens. Each ray of light is linearly polarized, passes through an immersion liquid, and then enters the body at a point along z = 0 within the y-z plane. The refractive index, nf, of the fluid is selected to match that of the specimen (for the wavelength of light used in the experiment) to minimize refraction due to surface distortion and deformation. For this example, the stress field is assumed to be two-dimensional, with shearing stresses rzYand r,, equal to zero, normal stress uz equal to zero, and normal stresses CJ-,and a,, independent of coordinate position z. Because the field is symmetrical about the y-z plane, rXv= 0 (within the plane) and a, and a,, are principal stresses. The direct problem of predicting data for a given stress field involves two steps: description of light ray propagation through the system, and determination of phase difference intensities at sample points p in the plane of observation. In the following two sections we consider each of these steps, respectively.
Reconstruction
of birefringent fields
353
plane
Fig.
2.
-
ray component with vibrations initially in the x-direction
-
ray component with vibrations initially in the y-direction
Light propagation
model for plane of symmetry within two-dimensional
field.
Light ray paths Each ray enters the field parallel to the z-axis, and is resolved component rays with vibrations assumed to be collinear secondary principal stress directions. We also assume that induced optical anisotropy can be related to the secondary stresses by the linearized Ramachandran-Ramaseshan model n,=n,+c,a,
+cza2
into two with the the fieldprincipal given by (1)
and n2 = n, + ClCT2+ C2CTl
(2)
where
and c is the speed of light in vacuum, u1 is the speed of light corresponding to ray 1 (with vibrations in the direction of secondary principal axis l), up is the speed of light corresponding to ray 2 (with vibrations in the direction of secondary principal axis 2), no is the refractive index of the material before loading, c,(A) and c,(A) are
354
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
stress-optic coefficients, and u1 and a2 are secondary principal stresses. ((TVand up are defined by two orthogonal axes which are normal to the direction of light propagation, and which corresponds to a shearing stress of zero.‘) The bending of each ray is governed by Fermat’s principle of minimum traveltime, and can be described mathematically by the well-known ray equation3’
where n can represent n, or n 2, r is the position vector, and s is the distance along a ray. By noting that dr/ds is the unit tangent vector, and by taking the dot product of eqn (4) with principal normal vector i,,, we obtain Vn .i, K=(5) n
where K is the curvature along a given ray path.30 This indicates that the curvature of a ray, at any particular point, is given by the normal component of the refractive index gradient divided by the refractive index at that point. Thus, ray bending is especially pronounced in regions where the stress components are rapidly changing in a direction normal to the direction of light propagation. In addition, the y-z plane is a plane of symmetry, so the xcomponent of eqn (4) becomes
(6) or d-x
n - = constant ds Atz=O,
dx/ds=O,
(7)
so
CLX -= 0 (8) d&Y throughout the medium, and trajectories remain within the y-z plane. Now, it would appear that eqns (l), (2), and (4) could be used directly to solve for light ray paths in terms of secondary principal stresses cl and u2. The actual situation, however, is somewhat more complex. Refer again to Fig. 2, and consider the ray entering the field at point A. Initially, at z = 0 within the body, one component vibration is in the x-direction, while the other is in the y-direction. Also, depending on
Reconstruction of birefringent fields
355
the field, the secondary principal stresses at z = 0 are given by U, = a, and u2 = u,,, or by u1 = aY and u2 = 0,. For the sake of discussion, let crl = a, and u2 = a,,. As each ray bends, the plane normal to the ray rotates by an angle 8. In the case of ray 1, vibrations remain in the x-direction (since the x-direction remains perpendicular to the ray) but for ray 2, the vibrations rotate by 8 and are no longer parallel to the y-axis. Also, by performing a transformation of coordinates, it can be shown that the secondary principal stresses are given by u1 = a,
(9)
and u2
=
uy
cos2
8
00)
In general, changes in the secondary principal stresses due to ray rotation cannot be neglected. That is, the coupling between ray bending and refractive index must be incorporated into ray tracing schemes. To determine a ray path by eqn (4) the refractive index field must be known; but the refractive index field is defined as in eqns (1) and (2) by secondary principal stresses which, in turn, depend on ray rotation. For the thin specimen of Fig. 2, however, we assume that 8 is small so that cos2 8 = 1, u2= a,,, and n2 =n,,. Accordingly, we may use the stress-birefringence relations 111= n, = Ro + c,u, + czu,
(II)
n2=n, ~n,+c,u,+c,u,
(12)
and These relations, together with the assumption of plane stress, imply that the refractive index field within the plane of symmetry is onedimensional, with n, = n,.(y) and n,, = n,(y). Lira and Vest’O have shown that for n = n(y) and small 8, solving eqn (4) yields a parabolic description of ray paths y = y(z) given by y
‘yo+;$
(13)
0
where y. is the y coordinate
of a ray at z = 0, and
(14) with 11’= dn/dy. If n’(y,) = 0, then y = yo, and the ray remains horizontal. As n’(yo) increases, so does the magnitude of correction term z2/2to. Ray paths within the specimen may be related to the stress field by
356
A&an T. Dolovich, Z. Wang, G. M. L. Gladwell
using eqns (11) and (12). For a ray with vibrations substituting eqn (11) into eqn (14) gives
in the x-direction,
where al = do,/dy and 0; = drr,/dy. Also, the quantity c,cr, + c2ay is very small in comparison with no. Although stress-induced anisotropy is suffkient to produce interferometric data, coefficients cl and c2 are generally of the order of 1O-‘o m’/N, and the actual changes in n (from n,) are relatively sma11.17Thus, ray trajectories with vibrations in the x-direction can be written as y, =
yo+;p
OSZSh
ox
where
(
n0
fox= c,u:+ c*o; >
(17) y’yo
and similarly, a ray trajectory with vibrations initially in the y-direction can be written as Yy
=yo+;;
OSZ";h
oy
(18)
where
The ray components
exit the specimen at y = h, where (Y&i
=yo+$
OX
(20)
and (Yy)h=Yo+;~
(21)
w Now, to completely describe each ray from z = 0 to z = h f L, we need to dist~guish between angle 8 and angles 1y and p, as shown in Fig. 3. Each angle is measured between the unit tangent vector and the z-axis, but each refers to a different part of the ray: 8 refers to 0 I z zs h, a refers to h I z s w, and p refers to w I z s h + L, where w is the distance between the entry surface of the specimen and exit surface of the region fluid. By using subscripts x and y to disting~sh between rays with vibrations initially in the x and y directions, respectively, we can identify six angles, & e,,, cys, (Ye, & and py for each ray pair.
357
Reconstruction of birefringent fields
plane
Fig. 3.
---
straight-ray theory
-
curvfJd-ray theory
Single ray path in plane of symmetry.
Angles 0, and 0, are not constants, but vary with position. They can be defined by differentiating eqns (16) and (18), respectively, to give 0, -tan
0, =%=F
(22) OX
and e,=tan
0, =%=F QY
(23)
The maximum values occur at z = h, where (24) and (25) At the boundary, z = h, between the specimen and the immersion fluid, each light ray refracts according to Snell’s law, and then remains at constant inclination (Y within the fluid. If we continue with our assumption that sin 8 - 8 and n = no, then no(eX), = nf sin (Y,
(26)
358
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
or
and similarly
(If the immersion fluid is selected so that nf = no, then cy, = (f!&)hand % = (UJ Once angles CY,and Q,, have been determined, ray paths within the fluid may be given by yX= (yA + (z - h) tan a,
h=zsw
(29)
yY= (y,L + (z - h) tan aY
hszsw
(30)
and where (Y~)~and (y,),, are defined by eqns (20) and (21), respectively. By a similar process, ray paths between the immersion fluid and the data plane can be given by Y, = (YA + (z - w) tan Px
wszsL+h
(31)
y, = (Y,L + (z - w) tan PY
wrzrL+h
(32)
and where (YJ, = (YA + (w - h) tan a,
(33)
(YA = (r,)h + (w - h) tan aY
(34)
5 sin/3,=-sina, n,
(35)
sin fly = -nf sin Q, n,
and n, is the ambient refractive index. Together, eqns (16)-(36) define the propagation of light rays through the system shown in Fig. 2. In the next section, we describe how the information obtained along these rays results in a tomographic data set. Note that Pindera and Hecker17 have used eqns (24), (25), (27), (28), (35) and (36) as the principle of ray deflection measurements. In developing their strain-gradient theory, however, they begin with very general expressions which are reduced in a number of stages to versions of increasing specificity. This process reveals the separate effect of each simplifying assumption, and establishes the conditions under which linearization of the Ramachandran-Ramaseshan model is reasonable.
Reconstruction of birefringentfields
359
This more detailed approach may be used in future analyses as curved ray inversion techniques are applied to three-dimensional fields. Interferometric data For the configuration shown in Fig. 2, the data set consists of an interference pattern observed at z = h + L.Because initially coincident rays bend and separate, however, the intensity at each data point, p, is a measure of the phase-difference between two rays with different points of origin along z = 0. This is in contrast to straight-ray theory where the fringe data are assumed to be caused by the superposition of collinear rays. The specific correspondence between rays at z = h + L depends on a number of system parameters including h, L,w,nf, and, of course, the stress field itself. Let A, and B, denote two rays which meet at a point p in the data plane. Ray A,,originates at point A and initially has vibrations in the y-direction; ray B, originates at point B and initially vibrates in the x-direction (Fig. 2). The two rays can produce interference at p because the analyser resolves the two perpendicular vibrations onto a common axis (not within the y-z plane). At point p, the dimensionless phase difference, iV, between the two rays is given by
(37) where uY is the velocity corresponding to ray A,, LJ, is the velocity corresponding to ray B,,and T is the period of the light source. Equation (37) can be rewritten as
(38) where (39) (40)
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
360
and A is the wavelength of the light source. Integrals 4,+ and 4B, are the optical path lengths of rays A,, and B,, respectively. Now, let YAYbe that part of ray A, within the specimen. Then
s,=/n,h+ lyAY
n,(w-h)+n,(L+h-w)
cmP,
cos cry
(41)
where LYEand p,, are given by eqns (28) and (36), respectively, with y, = yA. For the case d epicted in Fig. 2, the first term on the right-hand side of eqn (41) involves integration along paths of the form given by eqn (18). Consequently, it may be shown that”
I
n,,=..,[l+~(~~]
(42)
WA,
where nA = $(yA), and toy is given by eqn (19), with y, = yA. Equation (41) therefore becomes
h=nAh[l+;($)‘] +nf(w - h) cos(yy
n,(L + h - w) +
cospy
(43)
and similarly
(44) where n, = nx(yB), and cy,, fix, and tax are given by eqns (27), (35), and (17), respectively, with y0 = ye. Together, eqns (38), (43), and (44) define the interferometric data set, and complete our description of the direct problem. For each data point, p, the corresponding rays A,, and B, can be determined by using a numerical search in conjunction with the ray path theory developed in the previous section. Given a value of y,, the entry point of ray A, can be found by varying the value of y. in eqns (18), (30) and (32) until y,(h f L) = y,. Similarly, the entry point of ray B, can be obtained by varying the value of y. in eqns (16), (29) and (31) using a systematic approach. Alternatively, a pre-set number of ‘evenly spaced’ rays together with an interpolation scheme may be used. Finally, note that our physical model has been greatly simplified by assuming a collimated light source, and restricting the analysis to a plane of symmetry. Models for diffused light would be much more involved. Similarly, the direct problem for light propagation through a general stress field will pose a tremendous challenge. Potentially, it will embody the theories of integrated photoelasticity, and will require the
Reconstruction of birefringent fields
use of numerical dimensions.31*32
techniques
DATA
for ray tracing
ANALYSIS
361
in two-
and
three-
PROCEDURE
To incorporate the effects of ray curvature into the analysis of tomographic-birefringence data, an iterative approach is proposed. Since the light ray trajectories are unknown a priori, the first step in the procedure is to obtain the straight-ray solution to the field. Then, this may be used together with ray tracing theory to obtain a first approximation to the curved ray paths. Straight-ray inversion For a two-dimensional stress field, the classical theory assumes that yA = yB = y,, as shown in Fig. 3. Consequently, the fringe value N at data point p is taken to be a direct measure of the phase-difference between collinear rays in the specimen, and
(45) or, referring to eqns (11) and (12),
The values of c,(A) and c2(h) may be obtained from the literature,17 and inverting eqn (46) gives
(47) The individual quantities uX and ay may be determined by the method of oblique incidence, or by any other method for separating principal stresses.3 In the general tomographic case, eqn (46) is replaced by an expression involving integrals of stress along straight-line paths through the field. The unknown stress field is discretized and represented by a finite number of parameters, $; these parameters might be the pixel values within a mesh or the coefficients of an assumed form. By
362
Allan T. Dolovich, 2. Wang, G. M. L. Gladwell
discretizing the stress field and sampling the interferometric data set, the governing integral transform becomes a system of linear equations Af=b
(48)
where matrix A is known, and data vector b is measured. Also, to ensure that the system given by eqn (48) is overdetermined, data are generated for several viewing directions. The system is then solved for the unknown parameters. (In the classical theory of scalar field tomography, an unknown field is reconstructed from data corresponding to all straight-line paths through the field. Field reconstruction from partial or incomplete data, therefore, is an extremely important issue.) For each case to be considered in this paper, the stress field is one-dimensional, and the unknown parameter is the applied load. This allows us to demonstrate the data generation and data inversion processes using a single viewing direction. Moreover, the inversion formula given by eqn (47) immediately gives the straight-line solution to the problem, and therefore the zeroth iterate of the recursive procedure. Curved-ray inversion Let as’)(y) and @(y) be the stresses obtained using the straight-line inversion formula of eqn (47). By replacing CT,and a,, with u$‘) and a$“), respectively, in eqns (16)-(36), we can obtain a first approximation to the curved rays illustrated in Figs 2 and 3. That is, for each data point we can define kinematic quantities such as y$‘), yg), f&t),f$,), a:‘), a$!‘), p(O) where the superscript (0) indicates approximation based X 9 and p$O)“‘, on the zeroth iterate. Now, consider an .approximation to eqns (38), (43) and (44), given by
(49)
and
(51) In eqns (50) and (51), each kinematic quantity is defined on the basis of a$.‘) and OF). Substituting eqns (50) and (51) into (49) yields a single linear equation in terms of unknown quantities n$) and ng). Here, the
Reconstruction
of birefringent fields
363
superscript (1) indicates iterate number one. By writing a similar equation for each sample point p, we obtain a system of linear equations where the unknown quantities are refractive index values along the y-axis. If we represent n, = nX(y) and n, =n,(y) with assumed forms (such as 12,= a0 + a,y + a2y2 and IZ, = b, + b,y + b2y2, for example), we can write all values in terms of a common set of parameters (such as ao, a,, a2, bo, bl, and b2). Then, we may solve the system of linear equations to determine the refractive index components nk’) and n$?), and ultimately the stress components ai’) and ay). Similarly, the field components up) and a$“) corresponding to any iteration, k, can be obtained by defining all ray path parameters on the basis of (rik-‘) and a$‘-‘), and then using
and
to assemble a system of linear equations. In general, discretization of the curved-ray non-linear system of equations given by
problem
Af=b
leads to a
(55)
where each element of A is related to kinematic parameters, and is therefore a function of unknowns i. Within this context, each step of the proposed iterative approach involves solution of a linear system A(k-1$(k)
=
b
(56)
where the elements of matrix Ack-‘) are defined on the basis of previous iterate ftk-‘).
RESULTS
AND DISCUSSION
To investigate the effects of ray curvature and demonstrate the potential of iterative reconstruction methods, we consider two examples: a beam in pure bending and a diametrically loaded disc.
364
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell P
Y,
P
Y
7
I:
d
P
P
--id--
Fig. 4.
Example l-Beam
14, v = f
ii
A beam subjected to equal and opposite end-moments.
in pure bending
Consider a beam subjected to equal and opposite couples, as shown in Fig. 4. The beam is 250 mm long, has a depth of 2d = 60 mm, and a thickness of h = 10 mm. It is assumed to be made of Polyester Resin Palatal P6 (BASF), with no = 1.56, cl = -2.6 X lo-” m*/N, and c2 = -5.5 X 10-l’ m2/N, at A = 514.5 nm.” The beam is placed in the configuration shown in Fig. 2, with Dow Corning Pump Fluid #704 used as the immersion fluid, and 5 = 1.565 at A = 514.5 nm. The total width of the immersion tank is 100 mm, and w = 55 mm. At midspan, the beam is in pure bending, and the y-z plane is a plane of symmetry, see Fig. 4. Principal stresses gX and aY are given by
a=_Klf x
I
(57)
and uv = 0
(58)
respectively, where M is the bending moment and I is the area moment of inertia about the neutral axis. For the beam shown in Fig. 4, M = Pb, I = 2hd3/3, and r = y - d. The theory developed in the previous sections can be applied directly to real data obtained for planes of symmetry within two-dimensional fields. That is, the iterative approach for curved ray inversion is independent of analytical stress field solutions. Equations (57) and (58), however, may be used in conjunction with our data generation model to simulate physical experiments and artificially produce data. This is advantageous in that the iterative approach may be tested by comparing the reference solution (given by eqns (57) and (58)) to iterates produced by the curved ray method. Thus, we are able to focus
Reconstruction of birefringent fields
365
on difficulties associated with mathematical inversion. By starting with a reference solution, using our model to produce data, and then analysing the data on the basis of the same model, we are checking whether the curved ray tomographic operation is invertible. Successful solution of the mathematical inverse problem is necessary (but not sufficient) for accurate data analysis. It must, however, be coupled with the development of physical models from first principles, and experimental verification. A further advantage of computer simulation is the ease with which quantitative results may be produced for a wide range of parameters. These results may be used to design physical experiments and perform sensitivity analyses. For the beam shown in Fig. 4, data sets were generated numerically for a range of bending moments and two values of L: O-1m and l.Om. (Recall that L is the distance between the specimen and the data plane.) Values of M less than 240 N-m were used so that stresses would remain within the proportional limit of the material: upI = 40 MPa. The data were then analyzed under the a priori assumption of pure bending so that M could be taken as the only unknown parameter in the inversion process. The straight line solution, M(O), was obtained by rearranging eqn (47) to give M(O) =
NAI yjm
(59)
- 4
As can be seen from eqn (59), only one data point is theoretically needed to obtain straight-line approximation M(O)(although, in practice, a number of points would be used together with some form of regression analysis). In Fig. 5, the discrepancy between M(O)and applied bending moment, M, is plotted as a function of M, for a data point located at Yp= 0.2d (relative to the neutral axis). The plots show that the relative error increases with both M and L, and these trends can be explained on a physical basis.
0
60
120
Bending Moment M
Fig. 5.
Straight-ray
reconstruction
180
240
(N.m)
error for beam in pure bending with &Id = 0.2.
366
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
The magnitude
of M is a direct measure of the stress gradient, given
bY da M xdy --T
(60)
Also, recall from eqns (l), (2) and (5) that ray curvature is directly proportional to this gradient. By increasing M, ray deflection is increased, and the difference between a ray’s origin, yo, and its associated data point, y,, becomes more pronounced. The difference between y. and y, is also affected by L, which is a crude measure of imaging system characteristics. As L is increased, the rays have more opportunity to spread out and deflect before reaching the data plane. Thus, increasing either M or L will reduce the accuracy of straight-ray inversion. For this particular simulation, the error is relatively small. For L = O-1 m, it is less than l%, and for L = 1-Om, it is less than 10%. Also, these errors can be reduced by using a proper imaging system focused on the midplane of the specimen. Using a diffused light source might achieve some error reduction as well. (A curved-ray model incorporating these aspects should be developed.) We note, however, that the definition of acceptable accuracy depends on context. While an error of 10% might seem small within a classical context, it is unacceptable for contemporary techniques demanding precision measurement and data interpretation. To reduce the error numerically, the proposed iterative approach was applied to data generated for M = 240 Nem. For each iteration, k, substituting eqns (53) and (54) into eqn (52) produced a single expression with M w the only unknown. Plots of reconstruction error versus iteration are shown in Fig. 6, for L = O-1m and L = 1-Om. Each plot shows that the iterates quickly converge to a value which
2T
3
6
9
12
15
Number of Iterations Fig.
6.
Curved-ray
reconstruction
error for beam in pure bending and j$Jd = 0.2.
with M = 240 N.m
Reconstruction of birefringent fields
367
constitutes a significant improvement over M(O).The final error is very small and may be attributed to numerical round-off, only. Thus, the iterative method can potentially eliminate the modelZing error introduced by straight-ray inversion. In practice, this can represent a significant portion of the total error, which is influenced by a number of factors including the effects of noisy data. Example 2-Diametrically
loaded disc
To further test the curved ray approach, we consider a diametrically loaded disc as shown in Fig. 7. The disc has a diameter, D, of 80mm and a thickness, h, of 10mm. It is made of the same material as the beam in Example 1, and is surrounded by the same immersion fluid. For stresses within the y-z plane, the two-dimensional elasticity solution is given by33 3D
(61) (62) and rXy= 0
(63)
where P is the applied load per unit thickness. Equations (61) and (62) were used to numerically
Fig. 7.
A diametrically
loaded disc.
generate
data for
368
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
0.6
0.7
0.6
0.9
1
ii.lR Fig. 8.
Straight-ray
reconstruction
error for diametrically
loaded disc.
P = 90000 N/m. Then, in solving the inverse problem, we used the form given by eqns (61) and (62), but with applied load P unknown. This enabled us to obtain straight-ray solution P(O) and curved ray iterates P in much the same way we obtained M(O) and Mk) in the previous example. For each iteration, we used a single data point at y,, and solved for the unknown parameter. In this example, however, the nature of the field is quite different. The stress gradient in a beam in pure bending is uniform, but in a diametrically loaded disc the stress gradient increases with proximity to applied load P. We would therefore expect the straight-ray reconstruction error to vary significantly with the location of y,. As shown in Fig. 8, the error in P(O) increases dramatically as we use data closer to the top edge of the disc. At &/R = 0.95, the error is over 15% for L = O-1m and over 40% for L = 1-Om. Even for L = O-1m, the largest error is greater than 60%. The largest ray angle in this simulation is 0.16”. Thus, our assumption of small 8 is justified, and we can conclude that even small rotations can result in large reconstruction error if ray bending effects are ignored. We must be cautious in evaluating these results, however, since the elasticity solution deviates significantly from experimentally measured values near the applied load. 34 Theoretically, the stress field is twodimensional, and the magnitude of a,, is infinitely large at y = R, where the load is concentrated. This is contrary to physical experiments where the load P is applied over a finite area, and the stresses vary through the thickness of the disc. Our results, therefore, should be interpreted within the context of eqns (61) and (62). Errors of similar magnitude may be expected for straight-ray field solutions with gradients close to those implied by eqn (62). Also, we observe that the field points considered in this simulation are not as close to the applied load as the data points y,. As a result of ray bending, yA and y, are significantly
Reconstruction
0
3
369
of birefringent fields
6
9
12
15
Number of Iterations
Fig. 9.
Curved-ray
reconstruction
error for diametrically 0.9875.
loaded
disc with yp/R =
lower than y, (Fig. 3). Consequently, stress components a, and aY are less than the proportional limit, uPI, for all results shown. In Fig. 9, error plots for the curved-ray inversion process are presented. Again, the iterative reconstruction technique is effective, and iterates Pck) converge very rapidly. These simulations, however, should be regarded as only preliminary. In general, a tomographic analysis involves full-field discretization with many degrees of freedom, numerical integration techniques and ray tracing algorithms, specialized solution procedures for large systems of equations, and data sampling strategies for many viewing directions. At this scale, the problem becomes ill-posed, and numerical instability can lead to divergent reconstruction attempts. Even the first iterate can be extremely sensitive to small errors in the measured data. Fortunately, many of these issues have been addressed within the context of scalar field tomography and integrated photoelasticity. For example, regularization techniques and the singular value decomposi-. tion have been applied with success.‘3*35 In addition, strategies for convergence enhancement have been proposed.36 Thus, future application to strongly refracting birefringent media may take advantage of the experience already gained in these related areas.
CONCLUSION In this paper, we propose an iterative approach for reconstructing birefringent fields with correction for refraction. The equations presented are valid for planes of symmetry within two-dimensional fields, but the iterative concept is general and may be extended to tomography-based techniques such as integrated photoelasticity.
370
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
ACKNOWLEDGEMENT The authors gratefully acknowledge the funding provided by the Natural Sciences and Engineering Research Council of Canada.
REFERENCES 1. Guenther, R., Modern Optics. John Wiley & Sons, New York, USA, 1990, pp. 596-601. 2. Hecht, E., Optics (2nd edn). Addison-Wesley, Don Mills, 1987, pp. 3146. 3. Dally, J. W. & Riley, W. F., Experimental Stress Analysis (2nd edn). McGraw-Hill, New York, USA, 1978, pp. 406-531. 4. Pindera, M. J., Pindera, J. T. & Ji, X., Three-dimensional effects in beams: isodyne assessment of a plane solution. Exper. Mechanics, 29(l) (1989), 23-31. 5. Pindera, J. T. & Krasnowski, B. R., Theory of Elastic and Photoelastic Isodynes (SM Paper No. 184). Solid Mechanics Division, University of Waterloo, Waterloo, Canada, 1983. 6. Aben, H. K. & Kell, K. J., Integrated photoelasticity as tensor field tomography. Zeitschriji fur Angewandte Mathematik und Mechanik (ZAMM), 66(4) (1986), T118-9. 7. Aben, H. K., Idnurm, S. J., Josepson, J. I., Kell, K. E. & Puro, A. E., Integrated photoelasticity for residual stresses in glass specimen of complicated shape. SPZE Proc., 1554A (1991), 298-309. 8. Pindera, J. T., Advanced experimental mechanics in modem engineering science and technology. Trans. CSME, ll(3) (1987), 125-38. 9. Nievergelt, Y., Elementary inversion of Radon’s transform. SIAM Rev., 28( 1) (1986), 79-84. 10. Lira, I. H. & Vest, C. M., Refraction correction in holographic interferometry and tomography of transparent objects. Appl. Optics, 26(18) (1987), 3919-28. 11. Cha, S. & Vest, C. M., Tomographic reconstruction of strongly refracting fields and its application to interferometric measurement of boundary layers. Appl. Optics, 20(16) (1981), 2787-94. 12. Aben, H. K., Integrated Photoelasticity. McGraw-Hill, London, UK, 1979, pp. l-203. 13. Godbole, P. B. & Murty Gutta, R. K., Algorithms for integrated photoelasticity-axisymmetric cases. Exper. Mechanics, 29(l) (1989), 408. 14. Godbole, P. B., Chaudhari, U. M. & Bhave, S. K., A new integrated photoelasticity method for axisymmetric stress distribution. Exper. Mechanics, 27(l) (1987), 31-6. 15. Doyle, J. F. & Danyluk, H. T., Integrated photoelasticity for axisymmetric problems. Exper. Mechanics, H(6) (1978), 215-20. 16. Aben, H. K., Krasnowski, B. R. & Pindera, J. T., Nonrectlinear light
Reconstruction of birefringent fields
propagation CSME,
in integrated
8(4) (1984),
photoelasticity
of axisymmetric
371 bodies.
Trans.
195-200.
17. Pindera, J. T. & Hecker, 18.
19. 20. 21.
22.
23.
24. 25.
26.
27. 28.
F. W., Basic theory and experimental techniques of the strain-gradient method. Exper. Mechanics, 27(3) (1987), 314-27. Hecker, F. W., Kepich, T. Y. & Pindera, J. T., Non-rectilinear optical effects in photoelasticity caused by stress gradients. In Optical Methods in Mechanics of Solids, ed. A. Lagarde. Sijthoff and Noordhoff, Alpen aan den Rijn, The Netherlands, 1981, pp. 123-34. Aben, H. K., Integrated photoelasticity of axisymmetric bodies. Optical Engng, 21(4) (1982), 689-95. Post, D., Optical analysis of photoelastic polariscopes. Exper. Mechanics, 10(l) (1970), 15-23. Schomberg, H., An improved approach to reconstructive ultrasound tomography (Letters to the Editor). J. Phys. D: Appl. Phys., ll(15) (1978), L181-5. Lytle, R. J. & Dines, K. A., lterative ray tracing between boreholes for underground image reconstruction. IEEE Trans. Geosci. Remote Sensing, GE-18(3) (1980), 234-40. Lira, I. H., Correcting for refraction effects in holographic interferometry and tomography of transparent objects. PhD thesis, University of Michigan, Ann Arbor, MI, USA, 1987. Berryman, J. G., Stable iterative reconstruction algorithm for nonlinear traveltime tomography. Inverse Problems, 6 (1990), 21-42. Ramachandran, G. N. & Ramaseshan, S., Crystal optics. Handbuch der Physik (Vol. 25), ed. S. Flugge. Springer-Verlag, Berlin, Germany, 1961, pp. 202-17. Pindera, J. T., Foundations of experimental mechanics: principles of modelling observation and experimentation. In New Physical Trends in Experimental Mechanics, ed. J. T. Pindera. CISM Courses and Lectures No. 264, International Centre for Mechanical Sciences. Springer-Verlag, New York, USA, 1981, pp. 199-327. Pindera, J. T. & Cloud, G., On dispersion of birefringence of photoelastic materials. Exper. Mechanics, 6(9) (1966), 470-80. Pindera, J. T. & Straka, P., On physical measures of rheological responses of some materials in wide ranges of temperature and spectral frequency.
Rheol. Acta, 13 (1974), 338-51. 29. Dally, J. W. & Riley, W. F., Experimental
Stress Analysis (2nd edn). McGraw-Hill, New York, USA, 1978, pp. 475-8. 30. Born, M. & Wolf, E., Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. Pergamon Press, Oxford, UK, 1975, pp. 121-4. 31. Andersen, A. H. & Kak, A. C., Digital ray tracing in two-dimensional refractive fields. J. Acoust. Sot. Am., 72(5) (1982), 1593-606. 32. Trolinger Jr, J. D., Chipman, R. A. & Wilson, D. K., Polarization ray tracing in birefringent media. Optical Engng, 30(4) (1991), 461-6. 33. Volterra, E. & Gaines, J. H., Advanced Strength of Materials. PrenticeHall, Englewood Cliffs, NJ, USA, 1971, p. 186. 34. Pindera, J. T. & Pindera, M. J., Zsodyne Stress Analysis. Kluwer Academic Publishers, Boston, MA, USA, 1989, pp. 195-208.
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
372 35. Dubovikova,
E. A. & Dubovikov, M. S., Regularization, experimental and accuracy estimation in tomography and interferometry. J. Optical Sm. Am. A., 4(11) (1987), 2033-8. 36. Dolovich, A. T. & Gladwell, G. M. L., A generalized iterative approach to curved-ray tomography. Optics Lasers Engng, 17 (1992), 147-65. errors,
014%8166(94)ooo34-4
An Iterative Approach to Curved Ray Reconstruction of Birefringent Fields Allan T. Dolovich, University
of Saskatchewan,
Zhaoqing
Wang,
Department of Mechanical Engineering, Canada. S7N OWO
Saskatoon,
G. M. L. Gladwell University
of Waterloo,
Solid Mechanics Division, Waterloo,
Canada, N2L 3Gl
ABSTRACT The transmission of polarized light through non-homogeneous birefringent fields is examined, and an iterative approach for incorporating ray curvature into the analysis of tomographic data is proposed. The method extends strain-gradient theory and applies iterative algorithms which have been used by other researchers to reconstruct fields which are non-homogeneous but optically isotropic. Optical anisotropy is incorporated by using a light propagation model which includes both bending and separation of initially coincident light ray pairs. To illustrate the method, data are generated computationally and are analysed to determine the bending moment in a prismatic beam and the forces acting on a diametrically loaded disc.
INTRODUCTION Transmission methods of optical stress analysis are often based on the phenomenon of artificial birefringence (also known as temporary double refraction). ‘J This phenomenon is exhibited by certain transparent, polymeric materials which may be used to model stressed bodies of different geometry. If an initially isotropic, homogeneous model is subjected to an applied load, the resulting stress/strain field causes the material to become optically anisotropic. Consequently, a ray of linearly polarized light (that is, a light wave with a single, constant 347
348
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
plane of vibration) entering the field is split into two rays which are linearly polarized at right angles to each other. Each ray propagates with a different velocity and therefore a different refractive index. The refractive indices can be mathematically related to the stress state within the medium. Thus, measured changes in phase and direction of probing light rays can be used to extract information about the field. If sufficient information is obtained, the unknown stress field can be reconstructed. Classical transmission photoelasticity is perhaps the most wellknown method which makes use of this phenomenon. It relies, however, upon a substantially simplified view of light propagation within a non-homogeneous birefringent medium. For example, probing light rays are assumed to travel in straight lines. This assumption implies that originally coincident slow and fast light rays remain collinear, and their vibrations are superimposed at the data plane. If the stress field is assumed to be two-dimensional, the straight ray model also implies that: the axis of vibration for each ray remains constant in orientation through the thickness of the body, the components of stress influencing a given ray remain constant along that ray, and a point in the photoelastic body directly projects onto its corresponding point in the data plane (i.e. in the plane of the observed fringe pattern). These simplifications lead to an elementary description of the relationship between photoelastic fringe data N and principal stress differences uz. The dimensionless interferometric data N are considered to be (-rl a direct measure of the phase difference between collinear fast and slow rays. Moreover, a linearized stress-birefringence relationship or sfressoptic law is used, and it is assumed that the principal directions of stress and birefringence coincide. In the resulting mathematical model, fringe data N are simply proportional to principal stress differences g1 - c2, and the inverse problem of determining g1 - a2 from data N is mathematically trivial. These simplifications also facilitate the interpretation of isoclinics and the separation of principal stresses; at each point in the field, the isoclinic parameter is regarded as a direct measure of the principal stress orientation at that point. This simplified theory has been extremely useful within the context of its applications. However, a number of methods now demand a more complete description of the interaction between light and birefringent media. Examples include isodyne methods4*5 and procedures for tomographically reconstructing optically anisotropic fields.6 Applications include determination of residual stress in glass and plastics7 and the assessment of local effects.8 In particular, the reconstruction of three-dimensional birefringent
Reconstruction of birefringent fields
349
fields from tomographic data poses a tremendous challenge. Here, we use the term ‘tomography’ in a very general way to include any approach where an unknown field is reconstructed from data generated by field-induced alterations of probing radiation? The tomography of scalar, optically isotropic fields has been studied in detail and is well understood.lO*” In contrast, the tomography of stress-optic fields is much more complex, and a number of factors must be incorporated into the physical models used for data analysis. Images produced in the data plane represent extremely complicated integrated effects where both the magnitude and direction of the secondary principal stresses change continuously along any given light path. The method of integrated photoelasticity’* has been used to study these effects under the assumption of rectilinear light propagation, with axisymmetric problems of particular interest.1S’5 There is experimental evidence, however, that the ramifications of non-rectilinear light propagation cannot be ignored. In particular, it has been shown that the essential character of integrated photoelastic data depends on the object plane of the imaging system.16 Pindera and Hecker” have developed a strain-gradient theory which relates the derivatives of principal stresses to the bending of light in photoelastic specimens. The field-induced deflection of polarized laser beams through a specimen are measured, and the theory is used primarily to determine: the gradients of the sums and differences of principal stresses, stress-optic coefficients over a wide spectral range, and the errors which may result in classical or integrated methods which methods ignore the strain-gradient effect.‘* In addition, strain-gradient have shown some promise as stress reconstruction techniques. For example, Aben” has proposed a method for determining the complete stress state in cylinders by supplementing tomographic photoelastic data with ray deflection measurements. The contention here, however, is that strain-gradient theory should be incorporated into the theoretical fabric of tomographical methods which currently ignore ray bending. We believe that a reasonable objective is to work towards a post-processing approach for analysing integrated birefringence data in such a way that refraction effects are taken into account. When refraction effects are included, even the so-called direct problem of predicting data produced for a known stress field demands a higher level of sophistication regarding the mathematics and physics used to model the complex phenomena involved. Qualitative effects due to ray bending include the physical separation of initially collinear rays, the rotation of vibration axes, and the distortion of propagating
350
Allan T. Dolouich, Z. Wang, G. M. L. Gladwell
wavefronts. The effects of ray separation and wavefront distortion necessitate that non-classical parameters, such as the characteristics of the light source and imaging system, be included in the analysis. Specifically, the interaction between non-parallel rays of light depends on the plane of observation, and therefore the tomographic data set obtained for any particular case depends on the object plane of the imaging system. Also, because associated light components (originating at the same point of incidence) diverge and propagate along different paths, the data, in general, cannot be an exact measure of the interference between collinear rays, and coherence of the light source is a major issue. These complications are compounded by the inherent mathematical nonlinearity associated with integrated effects along curvilinear paths. Thus, the direct problem presents a considerable challenge. Accordingly, the inverse problem of reconstructing stress fields from tomographic data appears to be completely unmanageable. Since ray paths are functions of the unknown field, they are themselves unknown quantities, and the complexity of the direct problem would seem to preclude direct data analysis to simultaneously determine ray paths and field solutions. The spatial relation between a point in the data plane and its associated stress points in the body is unknown without some information regarding the rays; that is, the displacement of interferometric fringes obscures the field-data correspondence. Clearly then, at least an approximation to the ray paths is a prerequisite for inversion incorporating the effects due to refraction. The solution proposed here is an iterative approach. If, for a given problem, the straight-ray reconstruction is regarded as but a first approximation to the stress field, it may be used in conjunction with strain-gradient theory to determine an estimate of light ray paths, and a curved-ray inversion may be performed to yield (what is hoped to be) an improved approximation to the field. Then, in turn, this approximation may be used to redefine the rays and these can form the basis for the next inversion and so on until a suitable termination criterion is satisfied. This recursive process is illustrated in Fig. 1. This approach is conceptually identical to algorithms which have been used in ultrasound tomography,2’ geophysical imaging,” and holographic interferometry. 23 Within the context of holographic interferometry, it has been described as Curved Ray Algebraic Inversion (CRAI). It has been successfully applied to the reconstruction of optically isotropic fields, and its development has been followed by the introduction of other refraction-correction algorithms (see, for example, Refs 11 and 24). The objective here is to take a first step towards
Reconstruction of birefringent fields
351
Straight-Ray Reid Appmxlmation
Ray Appmximatbn
Fold Approximation
Fig.
1. Iterative scheme for curved-ray data analysis.
extending the CRAI approach to problems involving artificially birefringent media. It should be emphasized, however, that we are addressing but one issue within a larger scheme. Refraction-related processes are accompanied by other important phenomena which have been identified in the recent literature. For example, one of the fundamental assumptions underlying the Ramachandran-Ramaseshan model of artificial birefringence,25 and therefore the linearized stress-optic law, is that the optical axes of vibration coincide with the secondary principal directions of stress. The optical axes corresponding to the splitting of light rays in an anisotropic refractive index field can be predicted from first principles and Maxwell’s equations for electromagnetism.’ However, the co-directionality of these axes with those of secondary principal stress is assumed on the basis of a speculative, phenomenological approach.25 Thus, the classical stress-birefringence relations may not be appropriate in cases where the angular separation between the two sets of axes is significant. Advanced physical models of artificial birefringence should also include wavelength of the light source, A, as a major parameter. The of optical axes, and stress-optic refractive indices, orientation coefficients are all functions of A, and therefore the corresponding effects of dispersion, dispersion of optical axes, and dispersion of birefringence all contribute to the real response of photoelastic bodies.26-28 Furthermore, the effects of creep and stress relaxation in polymers must be considered. Linear viscoelasticity29 is often a reasonable assumption, but the limits of linearity must be established in terms of stress, strain, and birefringence. In summary, a complete description
352
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
of light propagation through a non-homogeneous, birefringent medium involves a complex interaction between these and other phenomena. By performing computer simulations, however, we can gain some experience with individual aspects of the system. In this paper, we present a data generation model which emphasizes the effects due to refraction. We then present an iterative approach to curved-ray analysis, and we demonstrate the approach by using two elementary examples: a beam in pure bending and a diametrically loaded disc. In each example, we use an analytical reference solution for the stress field together with our data generation model to numerically produce an interferometric data set. We then compare straight-ray analyses of these data to curved-ray solutions using the iterative method. In using computer simulation to address the inverse problem of curved-ray inversion, we take a preliminary step towards establishing its feasibility for practical applications. Mathematical inversion of the tomographic process is not the only aspect of curved-ray tomography, but it is a necessary condition for the successful analysis of real data obtained in physical experiments.
DATA
GENERATION
MODEL
To illustrate the concept of curved ray inversion, we may use a relatively simple example. Consider the propagation of polarized light within a plane of symmetry in a two-dimensional birefringent field, as shown in Fig. 2. The light source (not shown in Fig. 2) is a laser, followed by a pinhole placed at the focal point of a collimating lens. Each ray of light is linearly polarized, passes through an immersion liquid, and then enters the body at a point along z = 0 within the y-z plane. The refractive index, nf, of the fluid is selected to match that of the specimen (for the wavelength of light used in the experiment) to minimize refraction due to surface distortion and deformation. For this example, the stress field is assumed to be two-dimensional, with shearing stresses rzYand r,, equal to zero, normal stress uz equal to zero, and normal stresses CJ-,and a,, independent of coordinate position z. Because the field is symmetrical about the y-z plane, rXv= 0 (within the plane) and a, and a,, are principal stresses. The direct problem of predicting data for a given stress field involves two steps: description of light ray propagation through the system, and determination of phase difference intensities at sample points p in the plane of observation. In the following two sections we consider each of these steps, respectively.
Reconstruction
of birefringent fields
353
plane
Fig.
2.
-
ray component with vibrations initially in the x-direction
-
ray component with vibrations initially in the y-direction
Light propagation
model for plane of symmetry within two-dimensional
field.
Light ray paths Each ray enters the field parallel to the z-axis, and is resolved component rays with vibrations assumed to be collinear secondary principal stress directions. We also assume that induced optical anisotropy can be related to the secondary stresses by the linearized Ramachandran-Ramaseshan model n,=n,+c,a,
+cza2
into two with the the fieldprincipal given by (1)
and n2 = n, + ClCT2+ C2CTl
(2)
where
and c is the speed of light in vacuum, u1 is the speed of light corresponding to ray 1 (with vibrations in the direction of secondary principal axis l), up is the speed of light corresponding to ray 2 (with vibrations in the direction of secondary principal axis 2), no is the refractive index of the material before loading, c,(A) and c,(A) are
354
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
stress-optic coefficients, and u1 and a2 are secondary principal stresses. ((TVand up are defined by two orthogonal axes which are normal to the direction of light propagation, and which corresponds to a shearing stress of zero.‘) The bending of each ray is governed by Fermat’s principle of minimum traveltime, and can be described mathematically by the well-known ray equation3’
where n can represent n, or n 2, r is the position vector, and s is the distance along a ray. By noting that dr/ds is the unit tangent vector, and by taking the dot product of eqn (4) with principal normal vector i,,, we obtain Vn .i, K=(5) n
where K is the curvature along a given ray path.30 This indicates that the curvature of a ray, at any particular point, is given by the normal component of the refractive index gradient divided by the refractive index at that point. Thus, ray bending is especially pronounced in regions where the stress components are rapidly changing in a direction normal to the direction of light propagation. In addition, the y-z plane is a plane of symmetry, so the xcomponent of eqn (4) becomes
(6) or d-x
n - = constant ds Atz=O,
dx/ds=O,
(7)
so
CLX -= 0 (8) d&Y throughout the medium, and trajectories remain within the y-z plane. Now, it would appear that eqns (l), (2), and (4) could be used directly to solve for light ray paths in terms of secondary principal stresses cl and u2. The actual situation, however, is somewhat more complex. Refer again to Fig. 2, and consider the ray entering the field at point A. Initially, at z = 0 within the body, one component vibration is in the x-direction, while the other is in the y-direction. Also, depending on
Reconstruction of birefringent fields
355
the field, the secondary principal stresses at z = 0 are given by U, = a, and u2 = u,,, or by u1 = aY and u2 = 0,. For the sake of discussion, let crl = a, and u2 = a,,. As each ray bends, the plane normal to the ray rotates by an angle 8. In the case of ray 1, vibrations remain in the x-direction (since the x-direction remains perpendicular to the ray) but for ray 2, the vibrations rotate by 8 and are no longer parallel to the y-axis. Also, by performing a transformation of coordinates, it can be shown that the secondary principal stresses are given by u1 = a,
(9)
and u2
=
uy
cos2
8
00)
In general, changes in the secondary principal stresses due to ray rotation cannot be neglected. That is, the coupling between ray bending and refractive index must be incorporated into ray tracing schemes. To determine a ray path by eqn (4) the refractive index field must be known; but the refractive index field is defined as in eqns (1) and (2) by secondary principal stresses which, in turn, depend on ray rotation. For the thin specimen of Fig. 2, however, we assume that 8 is small so that cos2 8 = 1, u2= a,,, and n2 =n,,. Accordingly, we may use the stress-birefringence relations 111= n, = Ro + c,u, + czu,
(II)
n2=n, ~n,+c,u,+c,u,
(12)
and These relations, together with the assumption of plane stress, imply that the refractive index field within the plane of symmetry is onedimensional, with n, = n,.(y) and n,, = n,(y). Lira and Vest’O have shown that for n = n(y) and small 8, solving eqn (4) yields a parabolic description of ray paths y = y(z) given by y
‘yo+;$
(13)
0
where y. is the y coordinate
of a ray at z = 0, and
(14) with 11’= dn/dy. If n’(y,) = 0, then y = yo, and the ray remains horizontal. As n’(yo) increases, so does the magnitude of correction term z2/2to. Ray paths within the specimen may be related to the stress field by
356
A&an T. Dolovich, Z. Wang, G. M. L. Gladwell
using eqns (11) and (12). For a ray with vibrations substituting eqn (11) into eqn (14) gives
in the x-direction,
where al = do,/dy and 0; = drr,/dy. Also, the quantity c,cr, + c2ay is very small in comparison with no. Although stress-induced anisotropy is suffkient to produce interferometric data, coefficients cl and c2 are generally of the order of 1O-‘o m’/N, and the actual changes in n (from n,) are relatively sma11.17Thus, ray trajectories with vibrations in the x-direction can be written as y, =
yo+;p
OSZSh
ox
where
(
n0
fox= c,u:+ c*o; >
(17) y’yo
and similarly, a ray trajectory with vibrations initially in the y-direction can be written as Yy
=yo+;;
OSZ";h
oy
(18)
where
The ray components
exit the specimen at y = h, where (Y&i
=yo+$
OX
(20)
and (Yy)h=Yo+;~
(21)
w Now, to completely describe each ray from z = 0 to z = h f L, we need to dist~guish between angle 8 and angles 1y and p, as shown in Fig. 3. Each angle is measured between the unit tangent vector and the z-axis, but each refers to a different part of the ray: 8 refers to 0 I z zs h, a refers to h I z s w, and p refers to w I z s h + L, where w is the distance between the entry surface of the specimen and exit surface of the region fluid. By using subscripts x and y to disting~sh between rays with vibrations initially in the x and y directions, respectively, we can identify six angles, & e,,, cys, (Ye, & and py for each ray pair.
357
Reconstruction of birefringent fields
plane
Fig. 3.
---
straight-ray theory
-
curvfJd-ray theory
Single ray path in plane of symmetry.
Angles 0, and 0, are not constants, but vary with position. They can be defined by differentiating eqns (16) and (18), respectively, to give 0, -tan
0, =%=F
(22) OX
and e,=tan
0, =%=F QY
(23)
The maximum values occur at z = h, where (24) and (25) At the boundary, z = h, between the specimen and the immersion fluid, each light ray refracts according to Snell’s law, and then remains at constant inclination (Y within the fluid. If we continue with our assumption that sin 8 - 8 and n = no, then no(eX), = nf sin (Y,
(26)
358
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
or
and similarly
(If the immersion fluid is selected so that nf = no, then cy, = (f!&)hand % = (UJ Once angles CY,and Q,, have been determined, ray paths within the fluid may be given by yX= (yA + (z - h) tan a,
h=zsw
(29)
yY= (y,L + (z - h) tan aY
hszsw
(30)
and where (Y~)~and (y,),, are defined by eqns (20) and (21), respectively. By a similar process, ray paths between the immersion fluid and the data plane can be given by Y, = (YA + (z - w) tan Px
wszsL+h
(31)
y, = (Y,L + (z - w) tan PY
wrzrL+h
(32)
and where (YJ, = (YA + (w - h) tan a,
(33)
(YA = (r,)h + (w - h) tan aY
(34)
5 sin/3,=-sina, n,
(35)
sin fly = -nf sin Q, n,
and n, is the ambient refractive index. Together, eqns (16)-(36) define the propagation of light rays through the system shown in Fig. 2. In the next section, we describe how the information obtained along these rays results in a tomographic data set. Note that Pindera and Hecker17 have used eqns (24), (25), (27), (28), (35) and (36) as the principle of ray deflection measurements. In developing their strain-gradient theory, however, they begin with very general expressions which are reduced in a number of stages to versions of increasing specificity. This process reveals the separate effect of each simplifying assumption, and establishes the conditions under which linearization of the Ramachandran-Ramaseshan model is reasonable.
Reconstruction of birefringentfields
359
This more detailed approach may be used in future analyses as curved ray inversion techniques are applied to three-dimensional fields. Interferometric data For the configuration shown in Fig. 2, the data set consists of an interference pattern observed at z = h + L.Because initially coincident rays bend and separate, however, the intensity at each data point, p, is a measure of the phase-difference between two rays with different points of origin along z = 0. This is in contrast to straight-ray theory where the fringe data are assumed to be caused by the superposition of collinear rays. The specific correspondence between rays at z = h + L depends on a number of system parameters including h, L,w,nf, and, of course, the stress field itself. Let A, and B, denote two rays which meet at a point p in the data plane. Ray A,,originates at point A and initially has vibrations in the y-direction; ray B, originates at point B and initially vibrates in the x-direction (Fig. 2). The two rays can produce interference at p because the analyser resolves the two perpendicular vibrations onto a common axis (not within the y-z plane). At point p, the dimensionless phase difference, iV, between the two rays is given by
(37) where uY is the velocity corresponding to ray A,, LJ, is the velocity corresponding to ray B,,and T is the period of the light source. Equation (37) can be rewritten as
(38) where (39) (40)
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
360
and A is the wavelength of the light source. Integrals 4,+ and 4B, are the optical path lengths of rays A,, and B,, respectively. Now, let YAYbe that part of ray A, within the specimen. Then
s,=/n,h+ lyAY
n,(w-h)+n,(L+h-w)
cmP,
cos cry
(41)
where LYEand p,, are given by eqns (28) and (36), respectively, with y, = yA. For the case d epicted in Fig. 2, the first term on the right-hand side of eqn (41) involves integration along paths of the form given by eqn (18). Consequently, it may be shown that”
I
n,,=..,[l+~(~~]
(42)
WA,
where nA = $(yA), and toy is given by eqn (19), with y, = yA. Equation (41) therefore becomes
h=nAh[l+;($)‘] +nf(w - h) cos(yy
n,(L + h - w) +
cospy
(43)
and similarly
(44) where n, = nx(yB), and cy,, fix, and tax are given by eqns (27), (35), and (17), respectively, with y0 = ye. Together, eqns (38), (43), and (44) define the interferometric data set, and complete our description of the direct problem. For each data point, p, the corresponding rays A,, and B, can be determined by using a numerical search in conjunction with the ray path theory developed in the previous section. Given a value of y,, the entry point of ray A, can be found by varying the value of y. in eqns (18), (30) and (32) until y,(h f L) = y,. Similarly, the entry point of ray B, can be obtained by varying the value of y. in eqns (16), (29) and (31) using a systematic approach. Alternatively, a pre-set number of ‘evenly spaced’ rays together with an interpolation scheme may be used. Finally, note that our physical model has been greatly simplified by assuming a collimated light source, and restricting the analysis to a plane of symmetry. Models for diffused light would be much more involved. Similarly, the direct problem for light propagation through a general stress field will pose a tremendous challenge. Potentially, it will embody the theories of integrated photoelasticity, and will require the
Reconstruction of birefringent fields
use of numerical dimensions.31*32
techniques
DATA
for ray tracing
ANALYSIS
361
in two-
and
three-
PROCEDURE
To incorporate the effects of ray curvature into the analysis of tomographic-birefringence data, an iterative approach is proposed. Since the light ray trajectories are unknown a priori, the first step in the procedure is to obtain the straight-ray solution to the field. Then, this may be used together with ray tracing theory to obtain a first approximation to the curved ray paths. Straight-ray inversion For a two-dimensional stress field, the classical theory assumes that yA = yB = y,, as shown in Fig. 3. Consequently, the fringe value N at data point p is taken to be a direct measure of the phase-difference between collinear rays in the specimen, and
(45) or, referring to eqns (11) and (12),
The values of c,(A) and c2(h) may be obtained from the literature,17 and inverting eqn (46) gives
(47) The individual quantities uX and ay may be determined by the method of oblique incidence, or by any other method for separating principal stresses.3 In the general tomographic case, eqn (46) is replaced by an expression involving integrals of stress along straight-line paths through the field. The unknown stress field is discretized and represented by a finite number of parameters, $; these parameters might be the pixel values within a mesh or the coefficients of an assumed form. By
362
Allan T. Dolovich, 2. Wang, G. M. L. Gladwell
discretizing the stress field and sampling the interferometric data set, the governing integral transform becomes a system of linear equations Af=b
(48)
where matrix A is known, and data vector b is measured. Also, to ensure that the system given by eqn (48) is overdetermined, data are generated for several viewing directions. The system is then solved for the unknown parameters. (In the classical theory of scalar field tomography, an unknown field is reconstructed from data corresponding to all straight-line paths through the field. Field reconstruction from partial or incomplete data, therefore, is an extremely important issue.) For each case to be considered in this paper, the stress field is one-dimensional, and the unknown parameter is the applied load. This allows us to demonstrate the data generation and data inversion processes using a single viewing direction. Moreover, the inversion formula given by eqn (47) immediately gives the straight-line solution to the problem, and therefore the zeroth iterate of the recursive procedure. Curved-ray inversion Let as’)(y) and @(y) be the stresses obtained using the straight-line inversion formula of eqn (47). By replacing CT,and a,, with u$‘) and a$“), respectively, in eqns (16)-(36), we can obtain a first approximation to the curved rays illustrated in Figs 2 and 3. That is, for each data point we can define kinematic quantities such as y$‘), yg), f&t),f$,), a:‘), a$!‘), p(O) where the superscript (0) indicates approximation based X 9 and p$O)“‘, on the zeroth iterate. Now, consider an .approximation to eqns (38), (43) and (44), given by
(49)
and
(51) In eqns (50) and (51), each kinematic quantity is defined on the basis of a$.‘) and OF). Substituting eqns (50) and (51) into (49) yields a single linear equation in terms of unknown quantities n$) and ng). Here, the
Reconstruction
of birefringent fields
363
superscript (1) indicates iterate number one. By writing a similar equation for each sample point p, we obtain a system of linear equations where the unknown quantities are refractive index values along the y-axis. If we represent n, = nX(y) and n, =n,(y) with assumed forms (such as 12,= a0 + a,y + a2y2 and IZ, = b, + b,y + b2y2, for example), we can write all values in terms of a common set of parameters (such as ao, a,, a2, bo, bl, and b2). Then, we may solve the system of linear equations to determine the refractive index components nk’) and n$?), and ultimately the stress components ai’) and ay). Similarly, the field components up) and a$“) corresponding to any iteration, k, can be obtained by defining all ray path parameters on the basis of (rik-‘) and a$‘-‘), and then using
and
to assemble a system of linear equations. In general, discretization of the curved-ray non-linear system of equations given by
problem
Af=b
leads to a
(55)
where each element of A is related to kinematic parameters, and is therefore a function of unknowns i. Within this context, each step of the proposed iterative approach involves solution of a linear system A(k-1$(k)
=
b
(56)
where the elements of matrix Ack-‘) are defined on the basis of previous iterate ftk-‘).
RESULTS
AND DISCUSSION
To investigate the effects of ray curvature and demonstrate the potential of iterative reconstruction methods, we consider two examples: a beam in pure bending and a diametrically loaded disc.
364
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell P
Y,
P
Y
7
I:
d
P
P
--id--
Fig. 4.
Example l-Beam
14, v = f
ii
A beam subjected to equal and opposite end-moments.
in pure bending
Consider a beam subjected to equal and opposite couples, as shown in Fig. 4. The beam is 250 mm long, has a depth of 2d = 60 mm, and a thickness of h = 10 mm. It is assumed to be made of Polyester Resin Palatal P6 (BASF), with no = 1.56, cl = -2.6 X lo-” m*/N, and c2 = -5.5 X 10-l’ m2/N, at A = 514.5 nm.” The beam is placed in the configuration shown in Fig. 2, with Dow Corning Pump Fluid #704 used as the immersion fluid, and 5 = 1.565 at A = 514.5 nm. The total width of the immersion tank is 100 mm, and w = 55 mm. At midspan, the beam is in pure bending, and the y-z plane is a plane of symmetry, see Fig. 4. Principal stresses gX and aY are given by
a=_Klf x
I
(57)
and uv = 0
(58)
respectively, where M is the bending moment and I is the area moment of inertia about the neutral axis. For the beam shown in Fig. 4, M = Pb, I = 2hd3/3, and r = y - d. The theory developed in the previous sections can be applied directly to real data obtained for planes of symmetry within two-dimensional fields. That is, the iterative approach for curved ray inversion is independent of analytical stress field solutions. Equations (57) and (58), however, may be used in conjunction with our data generation model to simulate physical experiments and artificially produce data. This is advantageous in that the iterative approach may be tested by comparing the reference solution (given by eqns (57) and (58)) to iterates produced by the curved ray method. Thus, we are able to focus
Reconstruction of birefringent fields
365
on difficulties associated with mathematical inversion. By starting with a reference solution, using our model to produce data, and then analysing the data on the basis of the same model, we are checking whether the curved ray tomographic operation is invertible. Successful solution of the mathematical inverse problem is necessary (but not sufficient) for accurate data analysis. It must, however, be coupled with the development of physical models from first principles, and experimental verification. A further advantage of computer simulation is the ease with which quantitative results may be produced for a wide range of parameters. These results may be used to design physical experiments and perform sensitivity analyses. For the beam shown in Fig. 4, data sets were generated numerically for a range of bending moments and two values of L: O-1m and l.Om. (Recall that L is the distance between the specimen and the data plane.) Values of M less than 240 N-m were used so that stresses would remain within the proportional limit of the material: upI = 40 MPa. The data were then analyzed under the a priori assumption of pure bending so that M could be taken as the only unknown parameter in the inversion process. The straight line solution, M(O), was obtained by rearranging eqn (47) to give M(O) =
NAI yjm
(59)
- 4
As can be seen from eqn (59), only one data point is theoretically needed to obtain straight-line approximation M(O)(although, in practice, a number of points would be used together with some form of regression analysis). In Fig. 5, the discrepancy between M(O)and applied bending moment, M, is plotted as a function of M, for a data point located at Yp= 0.2d (relative to the neutral axis). The plots show that the relative error increases with both M and L, and these trends can be explained on a physical basis.
0
60
120
Bending Moment M
Fig. 5.
Straight-ray
reconstruction
180
240
(N.m)
error for beam in pure bending with &Id = 0.2.
366
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
The magnitude
of M is a direct measure of the stress gradient, given
bY da M xdy --T
(60)
Also, recall from eqns (l), (2) and (5) that ray curvature is directly proportional to this gradient. By increasing M, ray deflection is increased, and the difference between a ray’s origin, yo, and its associated data point, y,, becomes more pronounced. The difference between y. and y, is also affected by L, which is a crude measure of imaging system characteristics. As L is increased, the rays have more opportunity to spread out and deflect before reaching the data plane. Thus, increasing either M or L will reduce the accuracy of straight-ray inversion. For this particular simulation, the error is relatively small. For L = O-1 m, it is less than l%, and for L = 1-Om, it is less than 10%. Also, these errors can be reduced by using a proper imaging system focused on the midplane of the specimen. Using a diffused light source might achieve some error reduction as well. (A curved-ray model incorporating these aspects should be developed.) We note, however, that the definition of acceptable accuracy depends on context. While an error of 10% might seem small within a classical context, it is unacceptable for contemporary techniques demanding precision measurement and data interpretation. To reduce the error numerically, the proposed iterative approach was applied to data generated for M = 240 Nem. For each iteration, k, substituting eqns (53) and (54) into eqn (52) produced a single expression with M w the only unknown. Plots of reconstruction error versus iteration are shown in Fig. 6, for L = O-1m and L = 1-Om. Each plot shows that the iterates quickly converge to a value which
2T
3
6
9
12
15
Number of Iterations Fig.
6.
Curved-ray
reconstruction
error for beam in pure bending and j$Jd = 0.2.
with M = 240 N.m
Reconstruction of birefringent fields
367
constitutes a significant improvement over M(O).The final error is very small and may be attributed to numerical round-off, only. Thus, the iterative method can potentially eliminate the modelZing error introduced by straight-ray inversion. In practice, this can represent a significant portion of the total error, which is influenced by a number of factors including the effects of noisy data. Example 2-Diametrically
loaded disc
To further test the curved ray approach, we consider a diametrically loaded disc as shown in Fig. 7. The disc has a diameter, D, of 80mm and a thickness, h, of 10mm. It is made of the same material as the beam in Example 1, and is surrounded by the same immersion fluid. For stresses within the y-z plane, the two-dimensional elasticity solution is given by33 3D
(61) (62) and rXy= 0
(63)
where P is the applied load per unit thickness. Equations (61) and (62) were used to numerically
Fig. 7.
A diametrically
loaded disc.
generate
data for
368
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
0.6
0.7
0.6
0.9
1
ii.lR Fig. 8.
Straight-ray
reconstruction
error for diametrically
loaded disc.
P = 90000 N/m. Then, in solving the inverse problem, we used the form given by eqns (61) and (62), but with applied load P unknown. This enabled us to obtain straight-ray solution P(O) and curved ray iterates P in much the same way we obtained M(O) and Mk) in the previous example. For each iteration, we used a single data point at y,, and solved for the unknown parameter. In this example, however, the nature of the field is quite different. The stress gradient in a beam in pure bending is uniform, but in a diametrically loaded disc the stress gradient increases with proximity to applied load P. We would therefore expect the straight-ray reconstruction error to vary significantly with the location of y,. As shown in Fig. 8, the error in P(O) increases dramatically as we use data closer to the top edge of the disc. At &/R = 0.95, the error is over 15% for L = O-1m and over 40% for L = 1-Om. Even for L = O-1m, the largest error is greater than 60%. The largest ray angle in this simulation is 0.16”. Thus, our assumption of small 8 is justified, and we can conclude that even small rotations can result in large reconstruction error if ray bending effects are ignored. We must be cautious in evaluating these results, however, since the elasticity solution deviates significantly from experimentally measured values near the applied load. 34 Theoretically, the stress field is twodimensional, and the magnitude of a,, is infinitely large at y = R, where the load is concentrated. This is contrary to physical experiments where the load P is applied over a finite area, and the stresses vary through the thickness of the disc. Our results, therefore, should be interpreted within the context of eqns (61) and (62). Errors of similar magnitude may be expected for straight-ray field solutions with gradients close to those implied by eqn (62). Also, we observe that the field points considered in this simulation are not as close to the applied load as the data points y,. As a result of ray bending, yA and y, are significantly
Reconstruction
0
3
369
of birefringent fields
6
9
12
15
Number of Iterations
Fig. 9.
Curved-ray
reconstruction
error for diametrically 0.9875.
loaded
disc with yp/R =
lower than y, (Fig. 3). Consequently, stress components a, and aY are less than the proportional limit, uPI, for all results shown. In Fig. 9, error plots for the curved-ray inversion process are presented. Again, the iterative reconstruction technique is effective, and iterates Pck) converge very rapidly. These simulations, however, should be regarded as only preliminary. In general, a tomographic analysis involves full-field discretization with many degrees of freedom, numerical integration techniques and ray tracing algorithms, specialized solution procedures for large systems of equations, and data sampling strategies for many viewing directions. At this scale, the problem becomes ill-posed, and numerical instability can lead to divergent reconstruction attempts. Even the first iterate can be extremely sensitive to small errors in the measured data. Fortunately, many of these issues have been addressed within the context of scalar field tomography and integrated photoelasticity. For example, regularization techniques and the singular value decomposi-. tion have been applied with success.‘3*35 In addition, strategies for convergence enhancement have been proposed.36 Thus, future application to strongly refracting birefringent media may take advantage of the experience already gained in these related areas.
CONCLUSION In this paper, we propose an iterative approach for reconstructing birefringent fields with correction for refraction. The equations presented are valid for planes of symmetry within two-dimensional fields, but the iterative concept is general and may be extended to tomography-based techniques such as integrated photoelasticity.
370
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
ACKNOWLEDGEMENT The authors gratefully acknowledge the funding provided by the Natural Sciences and Engineering Research Council of Canada.
REFERENCES 1. Guenther, R., Modern Optics. John Wiley & Sons, New York, USA, 1990, pp. 596-601. 2. Hecht, E., Optics (2nd edn). Addison-Wesley, Don Mills, 1987, pp. 3146. 3. Dally, J. W. & Riley, W. F., Experimental Stress Analysis (2nd edn). McGraw-Hill, New York, USA, 1978, pp. 406-531. 4. Pindera, M. J., Pindera, J. T. & Ji, X., Three-dimensional effects in beams: isodyne assessment of a plane solution. Exper. Mechanics, 29(l) (1989), 23-31. 5. Pindera, J. T. & Krasnowski, B. R., Theory of Elastic and Photoelastic Isodynes (SM Paper No. 184). Solid Mechanics Division, University of Waterloo, Waterloo, Canada, 1983. 6. Aben, H. K. & Kell, K. J., Integrated photoelasticity as tensor field tomography. Zeitschriji fur Angewandte Mathematik und Mechanik (ZAMM), 66(4) (1986), T118-9. 7. Aben, H. K., Idnurm, S. J., Josepson, J. I., Kell, K. E. & Puro, A. E., Integrated photoelasticity for residual stresses in glass specimen of complicated shape. SPZE Proc., 1554A (1991), 298-309. 8. Pindera, J. T., Advanced experimental mechanics in modem engineering science and technology. Trans. CSME, ll(3) (1987), 125-38. 9. Nievergelt, Y., Elementary inversion of Radon’s transform. SIAM Rev., 28( 1) (1986), 79-84. 10. Lira, I. H. & Vest, C. M., Refraction correction in holographic interferometry and tomography of transparent objects. Appl. Optics, 26(18) (1987), 3919-28. 11. Cha, S. & Vest, C. M., Tomographic reconstruction of strongly refracting fields and its application to interferometric measurement of boundary layers. Appl. Optics, 20(16) (1981), 2787-94. 12. Aben, H. K., Integrated Photoelasticity. McGraw-Hill, London, UK, 1979, pp. l-203. 13. Godbole, P. B. & Murty Gutta, R. K., Algorithms for integrated photoelasticity-axisymmetric cases. Exper. Mechanics, 29(l) (1989), 408. 14. Godbole, P. B., Chaudhari, U. M. & Bhave, S. K., A new integrated photoelasticity method for axisymmetric stress distribution. Exper. Mechanics, 27(l) (1987), 31-6. 15. Doyle, J. F. & Danyluk, H. T., Integrated photoelasticity for axisymmetric problems. Exper. Mechanics, H(6) (1978), 215-20. 16. Aben, H. K., Krasnowski, B. R. & Pindera, J. T., Nonrectlinear light
Reconstruction of birefringent fields
propagation CSME,
in integrated
8(4) (1984),
photoelasticity
of axisymmetric
371 bodies.
Trans.
195-200.
17. Pindera, J. T. & Hecker, 18.
19. 20. 21.
22.
23.
24. 25.
26.
27. 28.
F. W., Basic theory and experimental techniques of the strain-gradient method. Exper. Mechanics, 27(3) (1987), 314-27. Hecker, F. W., Kepich, T. Y. & Pindera, J. T., Non-rectilinear optical effects in photoelasticity caused by stress gradients. In Optical Methods in Mechanics of Solids, ed. A. Lagarde. Sijthoff and Noordhoff, Alpen aan den Rijn, The Netherlands, 1981, pp. 123-34. Aben, H. K., Integrated photoelasticity of axisymmetric bodies. Optical Engng, 21(4) (1982), 689-95. Post, D., Optical analysis of photoelastic polariscopes. Exper. Mechanics, 10(l) (1970), 15-23. Schomberg, H., An improved approach to reconstructive ultrasound tomography (Letters to the Editor). J. Phys. D: Appl. Phys., ll(15) (1978), L181-5. Lytle, R. J. & Dines, K. A., lterative ray tracing between boreholes for underground image reconstruction. IEEE Trans. Geosci. Remote Sensing, GE-18(3) (1980), 234-40. Lira, I. H., Correcting for refraction effects in holographic interferometry and tomography of transparent objects. PhD thesis, University of Michigan, Ann Arbor, MI, USA, 1987. Berryman, J. G., Stable iterative reconstruction algorithm for nonlinear traveltime tomography. Inverse Problems, 6 (1990), 21-42. Ramachandran, G. N. & Ramaseshan, S., Crystal optics. Handbuch der Physik (Vol. 25), ed. S. Flugge. Springer-Verlag, Berlin, Germany, 1961, pp. 202-17. Pindera, J. T., Foundations of experimental mechanics: principles of modelling observation and experimentation. In New Physical Trends in Experimental Mechanics, ed. J. T. Pindera. CISM Courses and Lectures No. 264, International Centre for Mechanical Sciences. Springer-Verlag, New York, USA, 1981, pp. 199-327. Pindera, J. T. & Cloud, G., On dispersion of birefringence of photoelastic materials. Exper. Mechanics, 6(9) (1966), 470-80. Pindera, J. T. & Straka, P., On physical measures of rheological responses of some materials in wide ranges of temperature and spectral frequency.
Rheol. Acta, 13 (1974), 338-51. 29. Dally, J. W. & Riley, W. F., Experimental
Stress Analysis (2nd edn). McGraw-Hill, New York, USA, 1978, pp. 475-8. 30. Born, M. & Wolf, E., Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. Pergamon Press, Oxford, UK, 1975, pp. 121-4. 31. Andersen, A. H. & Kak, A. C., Digital ray tracing in two-dimensional refractive fields. J. Acoust. Sot. Am., 72(5) (1982), 1593-606. 32. Trolinger Jr, J. D., Chipman, R. A. & Wilson, D. K., Polarization ray tracing in birefringent media. Optical Engng, 30(4) (1991), 461-6. 33. Volterra, E. & Gaines, J. H., Advanced Strength of Materials. PrenticeHall, Englewood Cliffs, NJ, USA, 1971, p. 186. 34. Pindera, J. T. & Pindera, M. J., Zsodyne Stress Analysis. Kluwer Academic Publishers, Boston, MA, USA, 1989, pp. 195-208.
Allan T. Dolovich, Z. Wang, G. M. L. Gladwell
372 35. Dubovikova,
E. A. & Dubovikov, M. S., Regularization, experimental and accuracy estimation in tomography and interferometry. J. Optical Sm. Am. A., 4(11) (1987), 2033-8. 36. Dolovich, A. T. & Gladwell, G. M. L., A generalized iterative approach to curved-ray tomography. Optics Lasers Engng, 17 (1992), 147-65. errors,