younes-elhachi commented 3 years ago

I compared strains from DeformationGradientTensor to strains from xfab/tools/ubi_to_u_and_eps. It fails even in the reference co-ordinate system. Below is an example of my 1st test.

import numpy as np
from ImageD11.grain import grain
from ImageD11.finite_strain import e6_to_symm, DeformationGradientTensor as Ftensor
from xfab.tools import ubi_to_cell, ubi_to_u_and_eps

def strain_xfab(g1, g0):
    uc = ubi_to_cell(g0)
    U, E6 = ubi_to_u_and_eps(g1, uc)
    E = e6_to_symm(E6)
    e = np.dot(U, np.dot(E, U.T))
    return E, e

def test_strain(csys, n):
    ubis0 = [grain(np.random.random((3,3))).ubi for i in range(n)]
    ubis1 = [grain(np.random.random((3,3))).ubi for i in range(n)]
    E_Ftensor = [Ftensor(ubis1[i],ubis0[i]).finite_strain_ref() for i in range(n)]
    e_Ftensor = [Ftensor(ubis1[i],ubis0[i]).finite_strain_lab() for i in range(n)]
    E_xfab = [strain_xfab(ubis1[i],ubis0[i])[0] for i in range(n)]
    e_xfab = [strain_xfab(ubis1[i],ubis0[i])[1] for i in range(n)]
    if csys == 'reference':
        for i in range(n):
            assert np.allclose(E_Ftensor[i], E_xfab[i])
        return 'seems the strains in %s co-ordinate system are OK'%csys
    elif csys == 'lab':
        for i in range(n):
            assert np.allclose(e_Ftensor[i], e_xfab[i])
        return 'seems the strains in %s co-ordinate system are OK'%csys
    else:
        return 'csys is reference or lab'

test_strain('reference', 1000)
test_strain('lab', 1000)

jonwright commented 3 years ago

Sure. See : ImageD11/docs/DeformationGradientTensor.ipynb ImageD11/docs/NotesOnNewDeformationGradientTensor.ipynb They should match as infinitesimal strain, but otherwise they will not.

jonwright commented 3 years ago

This test was patched : ImageD11/test/eps_sig/test_eps.py

younes-elhachi commented 3 years ago

Sorry the infinitesimal displacement condition went right out of my mind. The two above documents are very informative, thanks! Same for the test you already did.

So maybe the only optional enhancement to make here is adding stress prediction from strain in the case of elastic law as it is the most common case. I am now thinking to add a separate module to go back and forth between strain and stress, giving the sample engineering parameters (Young modulus in each direction, Poisson ratios) or elastic stiffness parameters. I'll do that for the cases of:

[ ] isotropic,
[ ] orthotropic,
[ ] fully anisotropic material.

By the way, how to label something in GitHub as an enhancement rather than issue? Thank you.

jonwright commented 3 years ago

There's a menu on the right for tags. I think your existing code does the fully anisotropic case via the FitAllB formstiffnessMV thing? This can also be interesting for #77 - fitting an average stress state for a sample seems more physical than an average strain state.

jhektor commented 3 years ago

I started looking at the DeformationGradient notebooks tonight and I think there is something wrong in the Eulerian strains. There is a comment that "The finite strain page has the opposite sign for the Eulerian strain." It is not just the sign that is different, there is also an inverse on the term with the deformation gradients. I think the correct expression for the Eulerian strain is e=1/2*(I-F^(-T)F^(-1)). I haven't checked if that expression is used in the code somewhere, but it is at least not correct in the strains_from_deformation function in DeformationGradientTensor.ipynb

jonwright commented 3 years ago

Thanks Johan! I am afraid this could impact the code. The calculation in git master is currently doing:

grain.py -> computes "F" from finite_strain.py:DeformationGradientTensor
finite_strain_lab and finite_strain_ref are using F.F^T or F^T.F and powers of "m" from the Seth-Hill things
there is no version mentioning F^{-1}, but I think you can supply m = -1.

So if we are lucky there is a way to get this out somehow, but I have no idea how to check it. Do you have any reference numbers to use for testing and validation? I guess these things are well known in engineering mechanics but with small strains and small rotations I find it very hard to see if the numbers are right or wrong.

I am kind of hoping to validate/fix this stuff before doing a new release (which is getting a bit overdue now - the version on pypi is pretty dated)

jhektor commented 3 years ago

To me it looks like the finite_strain_lab with m=-1 will be correct, assuming that the np.linalg.matrix_power works with negative exponents.

I can try to dig out some old continuum mechanics book and look for examples with Eulerian strains. I can also do some testing, if I find something suitable.

jhektor commented 3 years ago

I haven't checked with a text book example yet, but I can verify that finite_strain_lab with m=-1 gives the same result as a my expression above ( e=1/2*(I-F^(-T)F^(-1))) on my data, and so it seem to be the correct Eulerian strain.

Another comment: Perhaps it would be better to default to m=1 in the finite_strain_ref (Green-Lagrange strain tensor) and m=-1 in finite_strain_lab (Almansi strain tensor).

jonwright commented 3 years ago

That sounds promising - it would be good to document this better. There are some testcases in test/test_finite_strain.py which contains the wonderfully cryptic comment "this is a bit misleading - negative m puts you in the other frame". Sadly I no longer remember what I had in mind when I wrote that, but it sounds like the thing you have run into?

The 0.5 was picked as a default in grain.py with an intention to try to pick the one which matches xfab/FitAllB more closely (where strain appears to get defined here https://github.com/FABLE-3DXRD/xfab/blob/fd3015dbd4145daf6a471765887c3aa139eb8bab/xfab/tools.py#L539). It seemed like it should be closer to the 0.5 power, but maybe I was confused. Manually testing now, I see that (m=1) is different to (m=0.5) from the xfab point of view, and xfab is apparently matching (m=0.5) for the example in the test_xfab_strain.

jadball commented 1 year ago

Hi @jonwright - I believe this has ramifications for stress tensors too. Just to confirm what I have in mind - if you use the ImageD11 finite_strain.py library to calculate grain strains, with m=0.5 as the default it'll give you the Biot strain? If so, what ramifications does this have for conversion to stress? Would you end up with the Biot stress tensor?

jonwright commented 1 year ago

Hi @jadball - I think @AxelHenningsson has been putting some of this into xfab recently here: https://github.com/FABLE-3DXRD/xfab/issues/42 Perhaps this is more fresh in his mind than mine right now.

Google could find a couple of relevant papers:

Acta Mech, S. N. Korobeynikov, "Basis-free expressions for families of objective strain tensors, their rates, and conjugate stress tensors" doi: 10.1007/s00707-017-1972-7

IJSS, Guo & Man, Conjugate stress and tensor equation doi:10.1016/S0020-7683(99)00209-7

My thinking on this is that any 'stress' measure depends on an accurate estimate for 'd-zero'. For typical metallurgy experiments: if the choice of 'm' for the strain tensor choice gives an important difference, then perhaps the d-zero is incorrect. These are usually small strains.

For cases where there is a large finite strain, then you might also need to account for non-linearity in the elastic constants? Perhaps this comes up for high pressure community (@smerkel ?). Going from 'strain' to 'stress' for a sample in a diamond anvil cell could be one example where the finite strain effects come up. I guess the lattice parameters and equation of state are more independent from the various finite strain measure definitions.

AxelHenningsson commented 1 year ago

It looks to me that finite_strain.py will give you the Green strain for m=0.5. (it uses m2 = int(round(m*2)) internally)? This is probably something that should be fixed in the docstring which implies that "the m you give is the m you get" (unless I am mistaken).

For different stress measures you are correct @jadball. For each strain measure, $\boldsymbol{E}$, there is a conjugate stress, $\boldsymbol{T}$, measure such that an increment in the stored energy of the body, $w=w( E )$, can be written as

$dw = \boldsymbol{T}^{(m)} d\boldsymbol{E}^{(m)}$

where $m$ denotes the selected strain. Thus the conjugate stress for the Biot strain is not the same as that of Almansi’s strain or the nominal strain or the Swainger’s strain or the Green’s strain etc, etc. They all must have their own conjugate stresses. When the deformation is very small, however, they should all reduce to the Cauchy small strain and stress tensors (else the measure of strain makes no sense).

jadball commented 1 year ago

@jonwright thanks for your advice here - it's reassuring to know any associated error is small for small applied strains. I played around in a Jupyter Notebook last night and seemed to confirm what you're saying.

@AxelHenningsson I'm not sure about your interpretation here. If m = 0.5, the following path should be used in finite_strain.py:

V,R,S = self.VRS # == F
Um = np.linalg.matrix_power(S, m2)
Em = (Um - np.eye(3))/m2
return Em

Given that m2 = 1:

V,R,S = self.VRS # == F
Um = S
Em = Um - np.eye(3)

This formula matches the Biot strain tensor definition here.

The challenge then is to calculate S (the right stretch tensor, also called U on the finite strain wiki) from F. Jon uses polar decomposition to do this (via F = R.S). You could also use the fractional_matrix_power method from scipy:

eps_biot = fractional_matrix_power(np.matmul(F.T, F), 0.5) - np.eye(e)

It's good to know that the stress tensor conjugates are also all slightly different! I think for small strains it should be negligible.

AxelHenningsson commented 1 year ago

Ah. Yeah. Sorry. I mistook the stretch, U, for the right Cauchy green. M=0.5 should be Biot like you say. 🙂

jonwright commented 1 year ago

Hi @AxelHenningsson : there should be some check in the code that (2m) is an integer, otherwise I wasn't sure how to compute these things.

Going back to the points from @jhektor above - I have a feeling there is still something in the code here that might not be right - but I don't know what it is. Something about the sign of the strain in the lab convention perhaps?

It would be good to clean up whatever else is wrong and then close this issue...

smerkel commented 1 year ago

Hi,

Regarding large strains as we would get in the high pressure community, John is correct that estimating d-zero is a key element.

We try to separate the stress (strain) tensors into an hydrostatic and non hydrostatic component as follows

The hydrostatic component is accounted for by adjusting the unit cell parameters, from which calculate a d-zero.
We evaluate the hydrostatic pressure using a known equation of state and unit cell volume.
We calculate the elastic constants at the relevant pressure, which we then use to convert deviatoric strains into deviatoric stresses.

This procedure is safer than trying to use second order elastic constants and trying to account for "large strains".

Deducing "hydrostatic" unit cell parameters, is not always straightforward. It works if we can make simple assumptions on the shape of the stress / strain tensor.

Sébastien

FABLE-3DXRD / ImageD11

Strain conventions #115

$dw = \boldsymbol{T}^{(m)} d\boldsymbol{E}^{(m)}$