|
|
UNDER CONSTRUCTION
|
|
|
> UNDER CONSTRUCTION
|
|
|
|
|
|
Photoelatic image: inverse method analysis
|
|
|
==========================================
|
|
|
|
|
|
In this section we mostly refer to Jonathan Kollmer [Git project](https://github.com/jekollmer/PEGS).
|
|
|
## 1. Overview
|
|
|
The inverse method analysis provide a quantitative estimation of the force network inside a granular media made in a photoelastic material. With such analysis the contact forces magnitude and orientation can be determined under some assumptions. The main idea of inverse method is to generate a numerical photoelastic picture that matches the experimental one, like illustrated below (from [1].
|
|
|
|
|
|
![Example](uploads/c4fd38c19a5c2c6263c786a28a3ff9bd/Example.png)
|
|
|
|
|
|
|
|
|
To achieve a good matching an optimization procedure is used, this procedure required a numerical computation of the photoelastic signal in the granular media. For cylindrical particles such an analytical expression of this signal can be obtained using the theory of elasticity and the theory of photoelasticity.
|
|
|
|
|
|
The different steps of this methods are:
|
|
|
1. Detecting the particles;
|
|
|
2. Getting an initial guess of contact forces;
|
|
|
3. Fitting the numerical and experimental photoelastic pictures.
|
|
|
|
|
|
|
|
|
## 2. Inverse method procedure
|
|
|
|
|
|
About the first step every details are given in the dedicated page and nothing has to be added. For the initial guess of the contact force, the $G^2$ value defined in the [gradient analysis](https://git-xen.lmgc.univ-montp2.fr/PhotoElasticity/Main/wikis/gradient-analysis) is used. Depending on the quality of the pictures, some manipulation based on this $`G^2`$ can be done in order to distribute the force over the contacts. A good initial force guessing increases the optimization procedure in speed and [results quality](https://git-xen.lmgc.univ-montp2.fr/PhotoElasticity/Main/wikis/gradient-analysis#1-overview). The forces determination occurs with an optimization procedure. This procedure is described hereafter as it is specific to the inverse method analysis.
|
|
|
|
|
|
### 2.1 Analytical expression of the photoelastic signal
|
|
|
|
|
|
All the steps needed to retrieve the stress field inside a cylindrical particle (called disks for the remainder of this document) submitted to multiple loads are available in different sources [1,2]. Here we briefly summarized the assumptions, the principle and give the final results.
|
|
|
|
|
|
The assumptions required to obtain the stress field inside the disk are:
|
|
|
1. Condition of equilibrium (force and torque balance);
|
|
|
2. Linear elasticity;
|
|
|
3. Plane strain;
|
|
|
4. Conditions of compatibility;
|
|
|
5. No body-forces (gravity...);
|
|
|
|
|
|
The stress field inside the disk is obtained by the superposition of:
|
|
|
1. A simple radial stress distribution for each load $`i`$: $`\sigma_{r_i} = -\frac{2f_{i}}{h\pi}\frac{\cos\theta_i}{r_i}`$;
|
|
|
2. An uniform tension to ensure a stress-free boundary: $`\sigma_{xx} = \sigma_{yy} = \sum_{i=1}^{N}\frac{f_{i}}{\pi h d}\cos\left(\theta_{1,i}+\theta_{2,i}\right)`$.
|
|
|
|
|
|
Of course for the summation, the stress tensors have to be expressed in the same coordinates system.
|
|
|
|
|
|
In a [circular polariscope](https://git-xen.lmgc.univ-montp2.fr/PhotoElasticity/Main/wikis/reflection-photoelasticity#21-photoelasticity-in-the-reflective-polariscope), the light intensity is a function of the difference of the two principal stress. This difference is equal to $`\sigma_{1}-\sigma_{2} = \sqrt{\left(\left(\sigma_{xx}-\sigma_{yy}\right)^2+\left(2\sigma_{xy}\right)^2\right)}`$ which is easely derivate from the [Mohr's circle](https://en.wikipedia.org/wiki/Mohr%27s_circle) for 2D stress.
|
|
|
|
|
|
### 2.2 Guessing initial forces value
|
|
|
|
|
|
From the position and dimensions of the disk, it is easy to construct a possible contact network based on geometrical consideration. A contact occurs if the distance between the two disk centers is lower than the sum of the two radius increased by a small tolerance to overcome imprecision from detection step. This provide the location of all possible contact forces acting on each disk, but it includes some false-contacts (due to the tolerance) or _non-force-bearing_ (which does not transmit any load).
|
|
|
|
|
|
A correct result of optimization procedure required a good initial force distribution. The [gradient analysis](https://git-xen.lmgc.univ-montp2.fr/PhotoElasticity/Main/wikis/gradient-analysis) provide the $`G^2`$ value which is proportional to the sum of the contact force magnitude $`\sum_i|F_i|`$ acting on the disk. To provide the initial guess of force distribution, two methods are available depending on the quality of the pictures.
|
|
|
|
|
|
__In the case of high quality pictures__ (resolution and contrast), a $`G^2_i`$ at each contact $`i`$ can be computed using a small region of interest located near the contact point. The false and _non-force-bearing_ contacts are eliminated if their $`G^2_i`$ value is lower than a threshold value on both contacting disks. The sum of the contact force magnitude is then distributed on the remaining contacts in proportion to the value of $`G^2_i`$. [1]
|
|
|
|
|
|
__In the case of low quality pictures__, each possible combination of active contact is provided to the optimization procedure. The sum of the contact force magnitude is equally distributed over the active contacts. This approach is slower than the previous one and only gives good results when the contact force are almost equally distributed over the active contacts. To overcome the second drawback the fitting procedure will consider all the granular media to propagate the fitted forces on the adjacent disk when the optimization procedure succeed (PhD Thesis of O. Lantsoght, currently defending).
|
|
|
|
|
|
### 2.3 Optimization procedure
|
|
|
|
|
|
The optimization procedure uses the forces magnitude and orientation as parameters and try to minimize the difference between the computed photoelastic signal and the picture of disk. Different methods exist to compare two images, among them we have the Mean Squared Error (MSE) or the Structural Similarity Index (SSIM) [3,4]. The main advantage of the MSE is the speed of computation but it considers the image globally. On the opposite he SSIM is slower but considers region of the image to take into account its structure. Presently, for this application, both comparison criteria are similar, they are globally good but for some case failed to detect that the images do not match at all.
|
|
|
|
|
|
> Here I would like to comute MSE and SSIM on the three particles given as example on the first figure.
|
|
|
|
|
|
The optimization algorithm must be able to handle non-linear problem and possibly constraints to impose force and torque balance. If the chosen algorithm does not handle constraints, the equilibrium condition can either be impose by reducing the free parameters by 3 [1], or simply be skipped and assuming that the fitted pictures impose the torque and force balance.
|
|
|
|
|
|
### 2.4 Force propagation
|
|
|
|
|
|
In the case of low quality pictures, the initial guess of the forces distribution may be not precise enough to converge toward a correct results. A correct result is identified by a threshold value to reach by the MSE or SSIM value.
|
|
|
|
|
|
To improve the initial guess of the forces provided to the optimizer, the forces identified on some disks are applied on the adjacent disk using the Newton's third law. Then the optimization procedure is called with all possible combination of active contact (the propagated contact force being always active). If an optimization procedure succeed, the newly identified forces are propagated to they neighbors.
|
|
|
|
|
|
> I have to work on the following images to explain propagation steps with particle 139 and/or 480.
|
|
|
<img src="uploads/cadedc89a2659bbcb5d0d927d5767e85/Picture.png" width="250">
|
|
|
<img src="uploads/093ba577f1853595e148c32d324fb0da/synthetic0.png" width="250">
|
|
|
<img src="uploads/7ec3e96162cda9ad086b560426877142/synthetic1.png" width="250">
|
|
|
<img src="uploads/5cf649d69e2a1180a8a59511b339ceb9/synthetic2.png" width="250">
|
|
|
<img src="uploads/7e5f9fa31da2d739f34f372157221460/synthetic3.png" width="250">
|
|
|
<img src="uploads/a0ea651b60494eb790cef8ef30a2dac9/synthetic4.png" width="250">
|
|
|
|
|
|
|
|
|
## 3. Examples
|
|
|
|
|
|
### Force propagation
|
|
|
The following pictures show the picture of an experiment (left) and the corresponding photoelastic signal (right). In this experiment the frame is vertically vibrated while the blender (stadium black shape at the center) oscillate by rotating around its upper extremity (that does not move vertically). In this experiment a circular bright-field polariscope is used, this is why the background is white.
|
|
|
|
|
|
<img src="uploads/6cb8360723ac75a52b43600e86d74bde/Wiki_Global.png" width="250">
|
|
|
<img src="uploads/31712778bfe2de3263c13d30a0b615cc/PhotoElast.png" width="250">
|
|
|
|
|
|
The quality of the picture does only allows to correctly identify good enough initial forces to converge the photoelastic signal on disks with simple force distribution (top left figure). The force propagation procedure allows to reach the convergence on more disks with more complex loading.
|
|
|
|
|
|
<img src="uploads/d5fa3d72a01cb149cfafe889ed594540/synthetic0.png" width="250">
|
|
|
<img src="uploads/ffb7930bcf8f1ab5fc19bebfb92726ad/synthetic1.png" width="250">
|
|
|
|
|
|
<img src="uploads/381567fa79c5355828309dfb539bbcae/synthetic3.png" width="250">
|
|
|
<img src="uploads/1f5ec507b8fef180cf598736c1f02c26/synthetic6.png" width="250">
|
|
|
|
|
|
Eventually the procedure did not succeed to reach convergences on all disks (bottom right figure) for reasons related to picture quality as well as assumption violation:
|
|
|
* Some disks show so close fringes that the picture of the signal is just gray, leading to completely wrong forces estimation;
|
|
|
* Some disks show larges deformations (up to 10% of their radius), in this situation the analytical expression if the stress is no more valid;
|
|
|
* The experiment is dynamic, some disks show residual constraints from a previous contact.
|
|
|
|
|
|
To improve the current result, it may be possible to manually provide some initial force for the optimization procedure for some disks. If the optimization succeed, the force propagation procedure can be restarted from the newly converged disks.
|
|
|
|
|
|
## 4. Reference
|
|
|
[1] K. E. Daniels, J. E. Kollmer, and J. G. Puckett, "Photoelastic force measurements in granular materials." Review of Scientific Instruments, vol. 88, no. 5, May 2017.
|
|
|
[2] S. Timoshenko and J. N. Goodier, Theory of Elasticity, 3rd ed. New-York: McGraw-Hill, 1970.
|
|
|
[3] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment : from error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, pp. 600–612, Apr. 2004.
|
|
|
[4] Z. Wang and A. C. Bovik, "Mean squared error: love it or leave it? - A new look at signal fidelity measures," IEEE Signal Processing Magazine, vol. 26, no. 1, pp. 98-117, Jan. 2009.
|
|
|
|
|
|
### 4.1 Implementation
|
|
|
|
|
|
Photo-elastic Grain Solver (PEGS):
|
|
|
1. Matlab implementation: [Git project](https://github.com/jekollmer/PEGS) [1].
|
|
|
2. Python implementation: [Git project](https://git.immc.ucl.ac.be/olantsoght/pegs_py) (based on [1]).
|
|
|
|
|
|
|
|
|
|
... | ... | |