Mesh-induced elastic anisotropy and distorting the quad. points as a way to palliate it

[David F]

Hello everyone,

I'm trying to solve a 2D solid mechanics homogenization problem, in which I have element-wise constant elastic properties, which are inhomogeneous and isotropic from element to element (i.e., I am assembling the system using the same 4-rank stiffness tensor for all the quadrature points of a certain element, but that tensor is different for each element). For this system, I would like to compute its effective elastic properties, which I do by loading the system under several different loading conditions as is done in standard homogenization approaches. The system should behave as an isotropic solid. However, I observe significant anisotropy (and clearly not due to random fluctuations that might arise because the element-to-element inhomogeneous properties are randomly distributed). I attribute this to a mesh dependency of the solution, since I have solved the same problem with a unstructured triangular mesh with another FEM package and I don't observe this issue. I believe the structured quadrilateral mesh induces some artificial elastic anisotropy, which is not there in the case of the unstructured triangular mesh due to its topological disorder.

I've thought of a way that might palliate this issue, which is to set different elastic properties at the quadrature points themselves (i.e., the properties are no longer element-wise constant). This seems to work to some extent since the system becomes less anisotropic, however it is not good enough.

Q1: is there a preferred way in dealII in which I could randomly distort a bit the location of the quadrature points? I think this extra distortion might help get rid of the mesh artifacts. Is is possible to do it with the in-built Lagrange linear FE or another type of FE is more suitable within dealII for this task? Basically I have no idea where to start from to do something like this, so any suggestion is welcome.

Q2: why the system behaves as anisotropic if its local inhomogeneous elastic properties are isotropic? If you have any comment or suggestion about the problem of mesh-induced elastic anistropy in FEM, I would like to know it.

Thanks in advance,

David.

[Bruno Turcksin]

David,

Unfortunately this is really not my domain of expertise but maybe you can randomize the cells https://dealii.org/current/doxygen/deal.II/namespaceGridTools.html#aec039d544c93a6e810427ad45ba0c84f instead of randomizing the quadrature points (which you cannot easily do).

Best,

Bruno

[Wolfgang Bangerth]

On 7/2/20 10:06 PM, David F wrote:
>
> *_Q2_:* why the system behaves as anisotropic if its local inhomogeneous

> elastic properties are isotropic? If you have any comment or suggestion about
> the problem of mesh-induced elastic anistropy in FEM, I would like to know it.

I don't know how exactly you choose your coefficient, but if you alternate
layers of isotropic materials, then you get an anisotropic material. Think
about layering styrofoam plates with steel plates -- the resulting stack of
layers is essentially incompressible under loads from the side (because the
steel plates provide the stiffness), but is quite compressible if you load it
from top and bottom (because the styrofoam layers will simply collapse).

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[David F]

I have a 2D system for which I create the stiffness tensor of an isotropic material, but for each finite element I create it with a different shear modulus. The shear modulus is random for each element (I use an exponential distribution, but any distribution leads to the same behavior as long as the std is high), with no structure such as layers or anything else. In this case, the system should clearly be macroscopically isotropic (up to statistical fluctuations due to the random properties) for symmetry reasons.

[Wolfgang Bangerth]

On 7/10/20 9:15 AM, David F wrote:
> I have a 2D system for which I create the stiffness tensor of an isotropic
> material, but for each finite element I create it with a different shear
> modulus. The shear modulus is random for each element (I use an exponential
> distribution, but any distribution leads to the same behavior as long as the
> std is high), with no structure such as layers or anything else. In this case,
> the system should clearly be macroscopically isotropic (up to statistical
> fluctuations due to the random properties) for symmetry reasons.

At least in the limit h->0 I agree. For finite mesh sizes, I would expect that
the material has a degree of anisotropy that goes to zero as you make the mesh
smaller. It is true that the axes of anisotropy should be oriented in random
ways for different realizations of the same experiment on the same mesh. When
you do your computations, have you checked (for different realizations of the
random process):
(i) whether the orientation of anisotropy is always the same, and always
related to the principal directions of the mesh?
(ii) how the magnitude of anisotropy behaves as you refine the mesh?

[David F]

I have made a somewhat extensive study on his issue and prepared some plots that will hopefully answer your questions, and also includes Bruno's suggestion about distorting the mesh. The basic setup is: I sheared the mesh along different orientations (see x-axis on the plots) and measured the shear modulus (y-axis). I have repeated the random process of setting the elastic properties many times to have good statistics (see errorbars on the plots). Each element has an isotropic stiffness tensor with a Poisson ratio of 1/3 and a shear modulus which is exponentially distributed with an average of 10. I use linear shape functions unless otherwise stated. If the picture are not big enough, you can find them in the links beneath them.

1) I change the resolution. By this I don't mean just a mesh with a bigger number of elements, but importantly each inhomogeneities is represented by a bigger number of elements. Therefore, we solve problems with exactly the same physical domain but with different resolution. In the legend, n means the resolution of the inhomogeneities. For n=1 each inhomogeneity is described by 1 FE. For n=2, by 2^2, and for n=4 by 4^2. We can see that for shearing with angle 0 (see pictures on the bottom for clarity) the shear modulus is minimum, while it is maximum for 45 degrees, when the principal axes are aligned with the mesh. The magnitude of the anisotropy is the difference between the maximum and the minimum. The difference decreases by increasing the resolution, but actually the relative difference is very similar, and it seems that by just increasing the resolution this problem won't go away. Finally, I have distorted the mesh, which doesn't change the behavior at all.

Link to the picture

2) I change the order of the shape functions. I use the original set up, in which each inhomogeneity is represented by 1 element. We see that increasing the shape function order has a somewhat similar effect as increasing the resolution of the inhomogeneity (expected, since in both cases we are increasing the number of dofs of each inhomogeneity). Therefore, increasing the order of the shape functions doesn't seem to be enough to fix the issue.

Link to the picture

3) I try different types of mesh. In this case, I compare the solution provided by dealII (i.e. quadrilateral mesh) with the solution obtained with a triangular mesh using the python FENICS package. Lastly, I solve using the triangular mesh but allocating the inhomogeneities in such a way that even if the mesh is triangular, the structure of the inhomogeneities looks quadrilateral. It seems that a triangular mesh, even it is structured, is able to provide an isotropic solution. Interestingly, the same mesh fails and behaves like a quadrilateral mesh if the structure depicted by the elastic properties look like those of the quadrilateral mesh.

Link to the picture

In summary, form all the tings that I tried, the only solution is to use a triangular mesh. Very interestingly, by arranging the elastic properties in a quad-like manner we can force the triangular mesh to give the wrong result of a quadrilateral mesh. Therefore, do you think it is possible to do the opposite, i.e., to make dealII's quadrilateral mesh "behave more triangular-like" by playing with elastic properties, quadrature points etc.?

Do you know the reason of the anisotropy in the quadrilateral mesh, or if there is some literature about this? I couldn't find anything.

Best,

David.

[David F]

There's a typo in the 3rd picture, the arrows of the triangular and the quad-like triangular should be interchanged. This is the right picture:

[Sebastian Stark]

Hi David,

I'm neither an expert, nor do I know the literature well, but looking on your pictures, I think, the situations you are studying are geometrically anisotropic. Just plot the distribution of angles the faces of your inhomogenities make with the x-axis. For the quad-case, you'll get two discrete peaks at 0 and 90 degree. for the triangular case, you get 0, 45 and 90 degree. So, from this, the results do not seem surprising to me (just consider the extreme case of cracks - if you have them only at 0 and 90 degree oriented, this is unlikely to be isotropic). The fact that you have randomly assigned elastic properties won't help to fix that.

A few examples (in 2d):

(1) equidistant circular inclusions in a matrix (matrix and inclusions two different isotropic linearly elastic materials) -> this should be isotropic if mesh size h->0

(2) equidistant square inclusions in a matrix, all aligned with x-axis (matrix and inclusions two different isotropic linearly elastic materials) -> probably anisotropic if mesh size h->0

(3) equidistant square inclusions in a matrix, random orientation of inclusions (matrix and inclusions two different isotropic linearly elastic materials) -> should be isotropic if mesh size h->0 and number of orientations->infinity

(4) equidistant square inclusions in a matrix, random orientation of inclusions (random isotropic linearly elastic materials) -> should be isotropic if mesh size h->0 and number of orientations->infinity

Also to consider: In 3d, there are 21 elastic constants for a linearly elastic material. In a mathematical 2d scenario, it should be 6. This suggests that, in example (3), it is not strictly necessary to have random orientation. Rather, a few (equally spaced) discrete orientations might be good enough. If that's the case, how many does one need? I'm betting on 6, not sure though. Alternatively, one could replace the square inclusions in (2) by regular polygons and ask how many vertices the polygon needs for isotropy. Again, I'm betting on 6.

Related: Are there crystal structures with such a high degree of symmetry, that they are elastically isotropic? For dielectric properties, a cubic crystal is good enough already. But the dielectric tensor is rank 2 and the elastic one rank 4. So you'll need much more symmetry in the crystal; and considering that a crystal can have 6-fold rotational symmetry at most that might be impossible.

What I'm just noticing: Hexagonal crystals are elastically isotropic perpendicular to the hexagonal axis. So, my bet on 6 might be good. And it might explain your observation that the triangular elements are relatively isotropic (though maybe not perfectly).

I hope that gives you some input. If you have definitive answers to any of the questions, I'm curious.

Regards,
Sebastian

Am 11.07.20 um 00:12 schrieb David F:

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/759917d0-227a-456e-bfac-00eab3ed71cdo%40googlegroups.com.

[David F]

Hi Sebastian,

first of all, sorry for my late reply. Thank you very much for comment, it certainly raises some very interesting points. I think the only thing which is still left to be explained is the fact that a triangular grid, which yields the right isotropic result, becomes anisotropic just by rearranging the elastic properties to form squares-like clusters. I think this case does not correspond to any of your 4 points since I do this without h->inf. In fact, each of those clusters are just 2 FE.

Best,

David.

To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.