Modelling Anistropic materials for electrostatic simulations

[phillip mobley]

Hello all,

This is a follow up discussion to the Jean-Paul answer to the question "Is the approach for the electrostatic Bi-linear form correct?". The question (and answer) can be found at the following link:

https://groups.google.com/forum/#!topic/dealii/i8P4JTwm7kQ

In short, I am modeling the maxwell equations for the electric field and voltage scalar field. The equations that I am using are displayed below:

div(E) = rho / epsilon where epsilon = epsilon_{0} * epsilon_{r} and rho is the charge density of the material.

-grad(V) = E

Using Dr. Bangerth's recommendation, I am solving for the scalar voltage field first then taking the gradient of my solution using the DataPostprocessorVector class. This has worked extremely well in my test program. For more information on that, see this post: https://groups.google.com/forum/#!topic/dealii/XIiPyMh0Jz4

However, now that I am actually coding the solver of the simulation, I will need to expand on my test simulation to include modeling anisotropic materials.

When the material is anisotropic, the epsilon value (the permittivity) of the material is represented by a tensor. To make things slightly simplified, I am only running 2D simulations.

I am attempting to determine the best method on modeling these types of materials. One approach that I have considered is to still solve for the voltage scalar field. If I go with this approach, then I will end up simulating a Possion equation where f = rho / epsilon and epsilon is a tensor. My concern for this direction would be that the Laplacian operator results in a scalar value. So I am not sure how I would handle the tensor on the RHS. Unless Deal.II has some sort of provision for this.

I have also been kicking around the idea of solving for the displacement vector D(i)= epsilon(i)(j) * E by substituting this equation into one of the equations above. Or at the very least, to use this equation as a constitutive relation to the equations above.

A third approach that I haven't quite explored very much is solving for the polarization of the material. But I am not sure if this is a practical approach since I could unnecessarily complicate the problem.

I wanted to post a discussion on this form to discuss what the best direction I should take to model anisotropic materials in Deal.ii? I have been looking at 3 different approaches and I would like to discuss which one of these 3 directions is the better. Or if there might be others that I have not considered yet.

[Wolfgang Bangerth]

Philip,

> In short, I am modeling the maxwell equations for the electric field and
> voltage scalar field. The equations that I am using are displayed below:
>

> div(*E*) = rho / epsilon where epsilon = epsilon_{0} * epsilon_{r} and

> rho is the charge density of the material.
>

> -grad(V) = *E*
> *
> *

> Using Dr. Bangerth's recommendation, I am solving for the scalar voltage
> field first then taking the gradient of my solution using the
> DataPostprocessorVector class. This has worked extremely well in my test
> program. For more information on that, see this post:
> https://groups.google.com/forum/#!topic/dealii/XIiPyMh0Jz4

> <https://groups.google.com/forum/#%21topic/dealii/XIiPyMh0Jz4>

>
> However, now that I am actually coding the solver of the simulation, I
> will need to expand on my test simulation to include modeling
> anisotropic materials.
>
> When the material is anisotropic, the epsilon value (the permittivity)
> of the material is represented by a tensor. To make things slightly
> simplified, I am only running 2D simulations.
>
> I am attempting to determine the best method on modeling these types of
> materials. One approach that I have considered is to still solve for the
> voltage scalar field. If I go with this approach, then I will end up
> simulating a Possion equation where f = rho / epsilon and epsilon is a
> tensor. My concern for this direction would be that the Laplacian
> operator results in a scalar value. So I am not sure how I would handle
> the tensor on the RHS. Unless Deal.II has some sort of provision for this.

Correct -- rho/epsilon doesn't make any sense.

The way to formulate this is to look at the mixed Laplace first. There,
if you start with the primal formulation
-div (K grad p) = f
and write it as a first order system, you arrive at the equation

K^{-1} \vec u - grad p = 0
div \vec u = -f

as is explained in step-20. If you have an anisotropic material
behavior, then K is not just a scalar, but in fact a matrix, and K^{-1}
is its inverse. For physical reasons, K will have to be a symmetric and
positive definite matrix, and consequently the same holds for K^{-1}.

Now how you implement something like this: This is actually already done
in step-20, where we describe K (or, rather, its inverse) as a dim x dim
tensor:
https://dealii.org/developer/doxygen/deal.II/step_20.html#Theinversepermeabilitytensor

This is then used in the matrix assembly here:
https://dealii.org/developer/doxygen/deal.II/step_20.html#MixedLaplaceProblemassemble_system

Best
WB

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

[phillip mobley]

Hello Dr. Bangerth,

Thank you for the reply to my discussion! I greatly I appreciate the help. I did not realize that the flow equations are very similiar if not the same to the electrostatic equations. I am still going through the material and formulating a plan on how to code the solver. There are a few other components that I would like to add to the data out put if possible.

From what I have read so far, the pressure would be the voltage scalar field and the velocity would be the electric field. The tutorial also mentions a flux. In electrostatics, this is the Electric Displacement Field.

The equations for finding these are exactly in step-20. But the symbols used are different.

I do have some questions regarding the material found in step-20 and with finite element as a whole. I apologize if these are beginner questions but there are many topics in finite elements that I am still learning.

1) Is Raviet-Thomas elements specific for mixed Laplace equations? I understand that step-20 is a tutorial and for learning. Other then being a good learning opportunity, is there any specific reason why RT elements were chosen for step-20?

1b) From another perspective, If I am solving for the mixed Laplace equation, do I need to setup the system to use the RT elements?

2) The weak derivation for step-20 has Dirichlet boundary values incorporated into it. I believe the sentence is "Note how in this formulation, Dirichlet boundary values of the original problem are incorporated in the weak form." I am guessing that this is because of the p = g on d(omega) term. Later down the road, I might need to implement different B.C. such as a periodic, anti-periodic, or a mixed coefficient B.C. If this is the case, would I need to re-derive the weak form? Or can I utilize the one found in step-20? Or, would it be better for me to derive a more general weak form that can be used with any B.C.?

3) From what I have read so far, I get the impression that the R.T elements are also incorporated into the weak form. Basically this sentence here: "Both xh and wh are from the space Xh=RT(k)×DQ(k), where RT(k) is itself a space"? is this thought correct? If so, if I wanted to change the element type later, would I need to re-derive the weak form?

4) The weak enforcement of the pressure boundary conditions, does this need to be factored in even if the B.C are not Dirichlet?

[Jean-Paul Pelteret]

Dear Philip,

Let me offer an alternative to what Wolfgang has suggested. I don’t think that it is necessary to use a mixed formulation to introduce material anisotropy into your problem. If you go back to the derivation in our previous discussion, there was this line:

(6) div [\epsilon_{0} \epsilon_{r} e] = \rho^{f} 

For tensor valued coefficient you could restate this in one of two ways. Either

(6a) div [ \epsilon_{r} e] = \rho^{f} / \epsilon_{0}

where \epsilon_{r} is a tensor of relative electric permittivity coefficients, or

(6b) div [ \epsilon e] = \rho^{f}

where \epsilon is a tensor of electric permittivity coefficients.

The latter, (6b) is analogous to linear elasticity. In fact, we have a code-gallery example (here’s a link to the theory pdf) that deals with anisotropic linear elasticity. I think that you would benefit greatly by reading the literature on piezoelectric materials, which are regularly modelled with a linear constitutive law. 

Touching on one of your other questions, I would say that nowadays many people would tend to use the electric scalar potential approach (solving for d with e as the independent field) to modelling these materials under electrostatic conditions, rather than the electric vector potential formulation (solving for e with d as the independent field). I’m not saying that one is more correct than the other - I’m just reporting what I perceive the trend in the literature to be (I’d be happy to have someone offer a different opinion that me on this). There is literally a mountain of literature dating back to the 1980’s on this topic, so I leave it up to you to read into this more.

Best,

Jean-Paul

--
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.
For more options, visit https://groups.google.com/d/optout.

[Wolfgang Bangerth]

Philip,

> The equations for finding these are exactly in step-20. But the symbols used
> are different.

Yes -- the beautiful universality of mathematics :-)

> 1) Is Raviet-Thomas elements specific for mixed Laplace equations? I
> understand that step-20 is a tutorial and for learning. Other then being a
> good learning opportunity, is there any specific reason why RT elements were
> chosen for step-20?

From the introduction of step-20:

"It is a well-known fact stated in almost every book on finite element theory
that if one chooses discrete finite element spaces for the approximation of
u,p inappropriately, then the resulting discrete saddle-point problem is
instable and the discrete solution will not converge to the exact solution.

To overcome this, a number of different finite element pairs for u,p have been
developed that lead to a stable discrete problem. One such pair is to use the
Raviart-Thomas spaces RT(k) for the velocity u and discontinuous elements of
class DQ(k) for the pressure p. For details about these spaces, we refer in
particular to the book on mixed finite element methods by Brezzi and Fortin,
but many other books on the theory of finite elements, for example the classic
book by Brenner and Scott, also state the relevant results."

> 1b) From another perspective, If I am solving for the mixed Laplace equation,
> do I need to setup the system to use the RT elements?

Yes. Or some of the other Hdiv elements.

> 2) The weak derivation for step-20 has Dirichlet boundary values incorporated
> into it. I believe the sentence is "Note how in this formulation, Dirichlet
> boundary values of the original problem are incorporated in the weak form." I
> am guessing that this is because of the p = g on d(omega) term. Later down the
> road, I might need to implement different B.C. such as a periodic,
> anti-periodic, or a mixed coefficient B.C. If this is the case, would I need
> to re-derive the weak form? Or can I utilize the one found in step-20? Or,
> would it be better for me to derive a more general weak form that can be used
> with any B.C.?

In general, whenever you integrate by parts, you get a boundary term and you
will have to think what to do with it. That's the place where you can apply
boundary values into the weak formulation. In the case of the mixed Laplace,
that just happens to be where Dirichlet boundary conditions go, just like with
the usual (non-mixed) formulation of the Laplace equation, the boundary term
allows you to take care of Neumann boundary conditions.

> 3) From what I have read so far, I get the impression that the R.T elements
> are also incorporated into the weak form. Basically this sentence here: "Both

> xhand whare from the space Xh=RT(k)×DQ(k), where RT(k)is itself a space"? is

> this thought correct? If so, if I wanted to change the element type later,
> would I need to re-derive the weak form?

No. You derive the weak formulation of the PDE, i.e., the continuous model. It
makes no reference to the (later step of) discretization.

> 4) The weak enforcement of the pressure boundary conditions, does this need to
> be factored in even if the B.C are not Dirichlet?

If you have Neumann boundary conditions, then you'll have to see what to do
with it in the context of the weak form you happen to have. I think I have a
whole video lecture on boundary conditions :-)

Best
W.

[phillip mobley]

Hello Jean-Paul,

Thank you for the links and resources. I am going through them now. I am not very familiar with linear elastic. I was able to better understand the step-20 and how it applies to electrostatics because I was able to relate the equations to their electrostatic equivalents.

I have been reading through you paper a couple of times now but I od not see how it relates to anisotropic materials. I could be missing something or misunderstanding some parts. From what I can tell from the paper, the strong form of linear elasticity is:

∇· σ + b = 0 on Ω.

where σ is the stress tensor which is modeled by σ = σ m + σ f. The isotropic linear constitutive law is given by: σm = C m : ε. Which is basically Hooke's law. Now, hooke's law does have an analogue for the polarization of a dielectric material when the material is exposed to an electric field. I am not very sure if this is the direction that I need to go in.

I have been looking up anisotropic materials for linear elasticity. I might be missing something but I am having a hard time seeing how this relates to the electrostatic case? Maybe I need to spend a little bit more time on this topic. I do see that there is a similarity between div [ \epsilon e] = \rho^{f} and ∇· σ + b = 0 where σ = \epsilon e and b = - \rho^{f}.

[phillip mobley]

Hello Dr. Bangerth,

Thank you for your response and clearing up those questions. I am continually going through steps 8, 20, and 21.

I ma currently exploring 2 directions, Jean-Paul linear elastic direction and the Mixed-Laplace direction. From what I can tell, modeling the anisotropic materials will not be a problem in deal.ii. Which is a nice thing.

[Wolfgang Bangerth]

On 11/6/18 8:10 PM, phillip mobley wrote:
>
> I ma currently exploring 2 directions, Jean-Paul linear elastic direction and
> the Mixed-Laplace direction. From what I can tell, modeling the anisotropic
> materials will not be a problem in deal.ii. Which is a nice thing.

Yes, definitely not a problem!

[Wolfgang Bangerth]

> I have been reading through you paper a couple of times now but I od not see
> how it relates to anisotropic materials. I could be missing something or
> misunderstanding some parts. From what I can tell from the paper, the strong
> form of linear elasticity is:
>
> ∇· σ + b = 0 on Ω.
>
> where σ is the stress tensor which is modeled by σ = σ m + σ f. The isotropic
> linear constitutive law is given by: σm = C m : ε. Which is basically Hooke's
> law. Now, hooke's law does have an analogue for the polarization of a
> dielectric material when the material is exposed to an electric field. I am
> not very sure if this is the direction that I need to go in.
>
>
> I have been looking up anisotropic materials for linear elasticity. I might be
> missing something but I am having a hard time seeing how this relates to the
> electrostatic case? Maybe I need to spend a little bit more time on this

> topic. I do see that there is a similarity between div *[* *\epsilon* *e*] =
> \rho^{f} and ∇· σ + b = 0 where σ = *\epsilon* *e *and b = - \rho^{f}.

You have
sigma = C eps(u)
and
-div sigma = b
which together yields
-div (C eps(u)) = b
and is the equivalent of
-div (K grad u) = f

The elastic material is anisotropic if the stress-strain tensor C represents
an anisotropic material. In the elasticity case, it isn't just a dxd matrix,
but a dxdxdxd (rank-4) tensor. So it looks more complicated, but it isn't really.

[phillip mobley]

Hello Dr. Bangerth,

Thank you stating this!

I was able to go through the paper from Jean-Paul more thoroughly a few times. I do have a few questions regarding the paper.

1) On equation 23, where do the subscripts ijkl come from? I understand that subscripts ij are from equation 5. But I do not see where the subscripts k and l factor into the problem.

2) Why is the shape function for the problem in the form as eqn. 23? Is this the basic definition of the shape function in deal.ii? Or is this specific to the problem?

3) How dos equation 27 reduce down from this:

to this:

Is this because of the Kronecker delta?

4) Also, I am not familiar with the abbreviation eps. Which math function does this relate to?

5) Since  -div (C eps(u)) = b is equivalent to  -div (K grad u) = , I am almost thinking I should do the Mixed Laplace approach for the problem? But, after reading the paper by Jean-Paul a few times, it would seem that the  -div (C eps(u)) = b approach does not require Mixed Laplace.

[phillip mobley]

I almost forgot:

6) I am starting to get the feeling that the difference between the isotropic and the anisotropic material is the stiffness matrix. This might be a long shot but would it be possible for me to utilize the same weak form as in the paper but change the code on the stiffness matrix? Or should I re-derive the weak form with anisotropic materials in mind?

6a) As a sub note on this one, there would be a concern since the weak form in the paper has a provision for "history variables". These variables also appear before the weak form and are apart of the derivation. Now since I am working on an electrostatic problem, I would not have these history variables and thus, I am unable to utilize this weak form.

[Wolfgang Bangerth]

Philip,

> I was able to go through the paper from Jean-Paul more thoroughly a few
> times. I do have a few questions regarding the paper.

So many questions :-)

Most of your questions -- like why the stress-strain tensor in
elasticity has four indices, and whether or not to use mixed elements --
are well answered in the literature. For the former, take a book about
solid mechanics. For the latter take a good book about the theory of
finite elements, but I think I also discuss this in one of the lectures
about how to choose the element for a given equation.

I, or J-P, or others could easily answer all of your questions here. But
my impression is that you are lacking some foundational knowledge about
common PDEs and their finite element discretization. Every answer we
will give you will raise more questions. The only way to get ahead of
this curve is to read.

Of course, if you have questions *specific to deal.II*, then we're more
than happy to answer them!

Best
Wolfgang

[phillip mobley]

Hello Dr. Bangerth,

Thank you very much for your feedback. I do appreciate you (and everyone else) for taking their time in answering the question.

Yes, it is true that there are many lackings in the theory of finite element. For the most part, I am self taught and doing my best to learn the material at hand. With that, I do have a final question, would you be able to recommend some books or authors that I can start looking into and reading?

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/ubRbrzYne4s/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+unsubscribe@googlegroups.com.

[Wolfgang Bangerth]

> Yes, it is true that there are many lackings in the theory of finite
> element. For the most part, I am self taught and doing my best to learn
> the material at hand. With that, I do have a final question, would you
> be able to recommend some books or authors that I can start looking into
> and reading?

I learned finite elements from a more formal, mathematical perspective
and so that's the kinds of books that I usually look things up in. The
ones on my shelf are Brenner & Scott, and Braess. Braess is often a bit
more intuitive, but it's still a math book at the end of the day. I also
have a German one by Grossmann & Roos that I like because it is more
written with applications in mind. This may not be useful to you because
of the language, of course.

There are many finite element books written from the engineering
perspective. I don't know them well. You may want to go to the library
and see whether they have anything by Zienkiewicz, Hughes, Oden,
Wriggers, or others.

Best
W.

[phillip mobley]

Great,

Thank you Dr. Bangerth, I will start to look from some of the book. Thank you for all of your help.

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/ubRbrzYne4s/unsubscribe.

To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.