**V. Kravcenkiene**^{1}
,
**S. Aleksiene**^{2}
,
**M. Vaidelys**^{3}

^{1, 2, 3}Kaunas University of Technology, Studentu str. 50-222, LT-51368 Kaunas, Lithuania

^{1}Corresponding author

Journal of Measurements in Engineering, Vol. 1, Issue 2, 2013, p. 106-112.

Received 11 April 2013; accepted 7 June 2013; published 30 June 2013

Copyright © 2013 Vibroengineering

A strategy for the selection of optimal smoothing parameter for hybrid experimental-numerical photoelastic experiment is proposed in this paper. This is not a straightforward procedure for photoelastic analysis since conventional finite element techniques are based on the approximation of nodal displacements, not stresses. Although the displacements are continuous at inter-element boundaries, the calculated stresses are discontinuous due to the operation of differentiation.

**Keywords: ** photoelasticity, smoothing, fringes.

Hybrid numerical–experimental photoelastic analysis is an attractive methodology evaluating the quality of the structural design in terms of the stress distribution. But conventional finite element analysis techniques are based on the approximation of nodal displacements (not stresses) via the shape functions. Ramesh et al. [1] have noted that photoelastic fringes can be effectively used for the detection of FEM meshing problems. Several different strategies for the selection of the photoelastic smoothing parameter have been proposed (a short overview of these strategies is given in the next section). Nevertheless a subjective and interactive human assessment of the reconstructed field of photoelastic fringes is required for all these methodologies before the final decision regarding the acceptable smoothing is taken. Therefore there exists a need for an effective objective data driven strategy for the selection of optimal smoothing parameter and this paper gives a qualitative description of such a technique.

A finite element norm [2] can be exploited for the adaptive selection of the smoothing parameter, which characterizes the error of reconstructed stress field. First, using the least squares method the nodal stress values can be identified in the global domain by minimizing the differences between the interpolated stress fields calculated on the basis of the values of nodal stresses and the theoretically calculated stress values based on the displacement field. Since the minimization is performed in the global domain and the interpolation is performed in the local domain of each element, the direct stiffness procedure can be applied for the described problem.

The common stress field calculation method is based on the global minimization of error calculated as the integral of the squared difference between the interpolated and factual stress fields:

(1)

$\sum _{D.S.}\left(\iint {\left(\left[N\right]\left\{s\right\}-\sigma \right)}^{2}dxdy\right),$
here $D.S.$ means FEM direct stiffness procedure; integrals are calculated in the domain of each finite element; $\left[N\right]$ is the shape function series of the finite element; $\left\{s\right\}$ is the column of the nodal values of the finite element stress components, $x$ and $y$ are the coordinates of the two-dimensional coordinate system.

The calculated nodal stress values, the finite element norms are calculated for each element as the average error of the field reconstruction, obtained by interpolating these nodal values. The first step of the calculation of nodal stress value creates a continuous stress field in the global domain. Nevertheless this field is purely suitable for visualization procedures, because the partial derivatives of the field are discontinuous on the inter-element boundaries, therefore during visualization the fringes become broken. The smoothing method proposed in [3] is based on the extension of the minimized error by a penalty for a rapid change of stress value in any analyzed direction:

(2)

${\lambda}_{i}\left({\left(\frac{\partial \sigma}{\partial x}\right)}^{2}+{\left(\frac{\partial \sigma}{\partial y}\right)}^{2}\right),i=1,...,p,$
here ${\lambda}_{i}>0$ are smoothing parameters selected individually for each finite element, $p$ – the number of finite elements in the global structure. The choice of the smoothing parameter is performed interactively from the qualitative view of the digital photoelastic images. When the parameter is too small the images are of unacceptable quality because of the unphysical behavior of the stress fields as a result of their calculation from the displacement formulation. When the parameter is too big an over-smoothed image is obtained which may look acceptable but be far from the real photoelastic image. So the best value of the parameter might be considered when most of the image is of acceptable quality without the unphysical behavior produced by the approximation. The smoothing parameter can be interpreted as a small penalty parameter to prevent the growth of derivatives of each component of the stresses:

(3)

$\frac{1}{2}\iint \left({\left(\left[N\right]\left\{{\delta}_{x}\right\}-{\sigma}_{x}\right)}^{2}+\lambda \left({\left(\frac{\partial {\sigma}_{x}}{\partial x}\right)}^{2}+{\left(\frac{\partial {\sigma}_{x}}{\partial y}\right)}^{2}\right)\right)dxdy$
$=\frac{1}{2}\iint \left({\left(\left[N\right]\left\{{\delta}_{x}\right\}-{\sigma}_{x}\right)}^{2}+\lambda {\left\{{\delta}_{x}\right\}}^{T}{\left[{B}^{*}\right]}^{T}\left[{B}^{*}\right]\left\{{\delta}_{x}\right\}\right)dxdy,$

$\frac{1}{2}\iint \left({\left(\left[N\right]\left\{{\delta}_{y}\right\}-{\sigma}_{y}\right)}^{2}+\lambda {\left\{{\delta}_{y}\right\}}^{T}{\left[{B}^{*}\right]}^{T}\left[{B}^{*}\right]\left\{{\delta}_{y}\right\}\right)dxdy,$

$\frac{1}{2}\iint \left({\left(\left[N\right]\left\{{\delta}_{xy}\right\}-{\tau}_{xy}\right)}^{2}+\lambda {\left\{{\delta}_{xy}\right\}}^{T}{\left[{B}^{*}\right]}^{T}\left[{B}^{*}\right]\left\{{\delta}_{xy}\right\}\right)dxdy,$

where $\lambda $ is the smoothing parameter; $\left\{{\delta}_{x}\right\}$ – the global vector of nodal values of $\left\{{\sigma}_{x}\right\}$; $\left\{{\delta}_{y}\right\}$ – the global vector of nodal values of $\left\{{\sigma}_{y}\right\}$; $\left\{{\delta}_{xy}\right\}$ – the global vector of nodal values of $\left\{{\tau}_{xy}\right\}$; $\left[N\right]$ – the row of the shape functions of the finite element; and $\left[{B}^{*}\right]$ – the matrix of the derivatives of the shape functions:

(4)

$\left[{B}^{*}\right]=\left[\begin{array}{ccc}\frac{\partial {N}_{1}}{\partial x}& \frac{\partial {N}_{2}}{\partial x}& \dots \\ \frac{\partial {N}_{1}}{\partial y}& \frac{\partial {N}_{2}}{\partial y}& \dots \end{array}\right].$
It can be noted that when the parameter of smoothing is too small the reconstructed images are of unacceptable quality because of the non-physical behavior of the stress field as a result of its calculation from the displacement formulation. When the parameter is too big an oversmoothed image is obtained, which may look acceptable but is far from the real photoelastic image. The components of the stresses can be interpolated from their nodal values by using the shape functions of the finite element. Then the components of strains ${\epsilon}_{x}$, ${\epsilon}_{y}$, and ${\gamma}_{xy}$ are obtained using those values of stresses and the matrix of elastic constants:

(5)

$\left\{\epsilon \right\}=\left\{\begin{array}{c}{\epsilon}_{x}\\ {\epsilon}_{y}\\ {\gamma}_{xy}\end{array}\right\}={\left[D\right]}^{-1}\left\{\begin{array}{c}{\sigma}_{x}\\ {\sigma}_{y}\\ {\tau}_{xy}\end{array}\right\}.$
Then the relative error norm for the $i$-th finite element can be calculated as [4]:

(6)

${\psi}_{i}=\frac{\underset{{e}_{i}}{\iint}{\left(\left\{\epsilon \right\}-\left[B\right]\left\{{\delta}_{0}\right\}\right)}^{T}\left[D\right]\left(\left\{\epsilon \right\}-\left[B\right]\left\{{\delta}_{0}\right\}\right)dxdy}{\underset{{e}_{i}}{\iint}{\left\{\epsilon \right\}}^{T}\left[D\right]\left\{\epsilon \right\}dxdy}.$
For those parts of the image where the relative error norms of the finite elements are too large, the image is expected to be of insufficient quality and a finer meshing or smoothing of the image with a larger value of the smoothing parameter may be required. The selection of the smoothing parameter $\lambda $ is then straightforward:

(7)

$\lambda =f\left(\sum _{i}{\psi}_{i}\right),$
where function $f$ can be selected as a linear growing function, $f\left(0\right)=0$. The slope of the function depends on the meshing and particularly on the type of the finite element used.

We will illustrate the proposed smoothing technique based on a simple one-dimensional example. Let us consider the following structure comprising 3 finite elements shown in Fig. 1.

**Fig. 1.**
The illustrative example comprising 3 quadratic finite elements

Nodal coordinates are set as follows: ${x}_{k}=k-\text{1}$, $k=\text{1,}\text{}\text{2}\text{,}\text{}\text{...,}\text{}\text{7}$. Then the shape functions for the $i$-th element read:

${N}_{1}^{\left(i\right)}\left(x\right)=\frac{1}{2}{\left(x-2i+\frac{1}{2}\right)}^{2}-\frac{1}{8},\mathrm{}{N}_{2}^{\left(i\right)}\left(x\right)=-{\left(x-2i+1\right)}^{2}+1,$

${N}_{3}^{\left(i\right)}\left(x\right)=\frac{1}{2}{\left(x-2i+\frac{3}{2}\right)}^{2}-\frac{1}{8},$

${N}_{3}^{\left(i\right)}\left(x\right)=\frac{1}{2}{\left(x-2i+\frac{3}{2}\right)}^{2}-\frac{1}{8},$

where $i=\text{1,}\text{}\text{...,3}$. The distribution of strain in the $i$-th element is expressed as:

(8)

${\epsilon}^{\left(i\right)}\left(x\right)={N}_{1}^{\left(i\right)}\left(x\right){\epsilon}_{2i-1}+{N}_{2}^{\left(i\right)}\left(x\right){\epsilon}_{2i}+{N}_{3}^{\left(i\right)}\left(x\right){\epsilon}_{2i+1}.$
Thus, stresses in the domain of the $i$-th element read:

(9)

${\sigma}^{\left(i\right)}\left(x\right)={B}_{1}^{\left(i\right)}\left(x\right){\epsilon}_{2i-1}+{B}_{2}^{\left(i\right)}\left(x\right){\epsilon}_{2i}+{B}_{3}^{\left(i\right)}\left(x\right){\epsilon}_{2i+1},$
where:

${B}_{1}^{\left(i\right)}\left(x\right)=\frac{\partial {N}_{1}^{\left(i\right)}\left(x\right)}{\partial x}=\frac{1}{2}-2i+x,\mathrm{}{B}_{2}^{\left(i\right)}\left(x\right)=\frac{\partial {N}_{2}^{\left(i\right)}\left(x\right)}{\partial x}=4i-2-2x,\mathrm{}$

${B}_{3}^{\left(i\right)}\left(x\right)=\frac{\partial {N}_{3}^{\left(i\right)}\left(x\right)}{\partial x}=\frac{3}{2}-2i+x.$

${B}_{3}^{\left(i\right)}\left(x\right)=\frac{\partial {N}_{3}^{\left(i\right)}\left(x\right)}{\partial x}=\frac{3}{2}-2i+x.$

Let us denote the unknown nodal values of stresses as ${s}_{k}$. Then:

$\left[N\right]=\left[{N}_{1}^{\left(i\right)}\left(x\right);{N}_{2}^{\left(i\right)}\left(x\right);{N}_{3}^{\left(i\right)}\left(x\right)\right],\left[C\right]=\left[{B}_{1}^{\left(i\right)}\left(x\right);{B}_{2}^{\left(i\right)}\left(x\right);{B}_{3}^{\left(i\right)}\left(x\right)\right],$

and:

(10)

$\underset{2\left(i-1\right)}{\overset{2i}{\int}}{\left[N\right]}^{T}\left[N\right]dx=\left[\begin{array}{ccc}\frac{4}{15}& \frac{2}{15}& -\frac{1}{15}\\ \frac{2}{15}& \frac{16}{15}& \frac{2}{15}\\ -\frac{1}{15}& \frac{2}{15}& \frac{4}{15}\end{array}\right],\mathrm{}\underset{2\left(i-1\right)}{\overset{2i}{\int}}{\left[C\right]}^{T}\left[C\right]dx=\left[\begin{array}{ccc}\frac{7}{6}& -\frac{4}{3}& \frac{1}{6}\\ -\frac{4}{3}& \frac{8}{3}& -\frac{4}{3}\\ \frac{1}{6}& -\frac{4}{3}& \frac{7}{6}\end{array}\right],$
$\underset{2\left(i-1\right)}{\overset{2i}{\int}}{\left[N\right]}^{T}{\sigma}^{\left(\text{i}\right)}dx=\left[\begin{array}{c}-\frac{1}{2}{\epsilon}_{2i-1}+\frac{2}{3}{\epsilon}_{2i}-\frac{1}{6}{\epsilon}_{2i+1}\\ -\frac{2}{3}{\epsilon}_{2i-1}+\frac{2}{3}{\epsilon}_{2i+1}\\ \frac{1}{6}{\epsilon}_{2i-1}-\frac{2}{3}{\epsilon}_{2i}+\frac{1}{2}{\epsilon}_{2i+1}\end{array}\right].$

Nodal values of stresses are computed from the linear algebraic system (Eq. (3)). Then the smoothed field of stresses in the domain of the $i$-th element reads:

(11)

${S}^{\left(i\right)}\left(x,\lambda \right)={s}_{2i-1}\left(\lambda \right){N}_{1}^{\left(i\right)}\left(x\right)+{s}_{2i}\left(\lambda \right){N}_{2}^{\left(i\right)}\left(x\right)+{s}_{2i+1}\left(\lambda \right){N}_{3}^{\left(i\right)}\left(x\right).$
**Fig. 2.**
Discontinuous theoretical stress field (the solid line) and the regularized stress field (the dashed line)

Let us assume that $\epsilon =\left(\begin{array}{ccccccc}-0.01& 0& 0.01& -0.02& 0.04& 0.02& -0.06\end{array}\right)\text{.}$ Then the graphical representation of the theoretical and regularized fields of stress is shown in Figure 2. Note that the regularized field of stress is constructed as the interpolation of nodal values of stresses, computed as the arithmetic average of the values of the theoretical field of stresses at corresponding nodes.

Now we minimize the residual $R$ computed as the difference between the nodal values of the regularized field (denoted as ${\sigma}_{reg}$) and the smoothed nodal values ${S}_{\lambda}$:

(12)

$\underset{\lambda >0}{\mathrm{m}\mathrm{i}\mathrm{n}}R=\mathrm{m}\mathrm{i}\mathrm{n}\sqrt{\frac{1}{n}\sum _{i}^{n}{({\sigma}_{reg}-{S}_{\lambda})}^{2}}.$
The smoothed fields of stresses at $\lambda =$0 and $\lambda =$0.032 are illustrated in Figure 3.

Finally it can be noted that $\lambda =$0.032 is the optimal value of the photoelastic smoothing – the variation of $R$ with respect to $\lambda $ is shown in Figure 4.

**Fig. 3.**
The smoothed field of stresses at $\lambda =$0 (the thin solid line) and at $\lambda =$0.032 (the thick solid line) shown together with the regularized field (the dashed line)

**Fig. 4.**
Minimation of the residual $R$ with respect to the smoothing parameter $\lambda $

The proposed technique for the optimal selection of the smoothing parameter can be effectively implemented for a more complex FEM analysis. We will illustrate this technique using a 2D FEM structure comprising 9 Lagrange finite elements with quadratic interpolation. The theoretical field of stresses is continuous in the domain of every element, but discontinuous at inter-element boundaries (Figure 5).

**Fig. 5.**
The theoretical field of stresses is continuous in the domain of each finite element, but is discontinuous at inter-element boundaries

The variation of the residual $R$ (Eq. (12)) with respect to $\lambda $ is shown in Figure 6; the optimal value of the smoothing coefficient is ${\lambda}_{opt}=$0.124.

**Fig. 6.**
The minimal value of the residual $R$ is reached at ${\lambda}_{opt}=$0.124

Finally the smoothed fields of stresses at $\lambda =$0 and ${\lambda}_{opt}=$0.124 are illustrated in the Figures 7a and 7b.

**Fig. 7.**
The smoothed fields of stresses at $\lambda =$0 (part a) and ${\lambda}_{opt}=$0.124 (part b)

a)

b)

The digital smoothing of photoelastic stresses produced by the finite element method is an important part of hybrid experimental-numerical methods for measurement of engineering structures. A number of different approaches for qualitative smoothing of the photoelastic fringes had been proposed during the last decade. Our approach belongs to the same group of problems, but the value of the smoothing parameter is selected not on the basis of the subjective visual interpretation of the smoothed image, but on the minimization of the residual assessing differences among nodal values of theoretical and regularized stresses. Such an approach fits well into the modular structure of finite element algorithms and can be effectively exploited in different photoelastic measurement scenarios.

**Ramesh K., Pathak P. M.**Role of photoelasticity in evolving discretization schemes for FE analysis. Journal of Experimental Techniques, Vol. 23, Issue 4, 1999, p. 36-38.**Bathe K. J.**Finite Element Procedures in Engineering Analysis. Prentice-Hall, New Jersey, 1982.**Ragulskis M., Ragulskis L.**Conjugate approximation with smoothing for hybrid photoelastic and FEM analysis. Communications in Numerical Methods in Engineering, Vol. 20, Issue 8, 2004, p. 647-653.**Ragulskis M., Ragulskis L.**Plotting isoclinics for hybrid photoelasticity and finite element analysis. Experimental Mechanics, Vol. 44, Issue 3, 2004, p. 235-240.