Transient 2D piezoelectric problem

Hello All,

I am currently working on trying to implement a transient piezoelectric system in 2D. I have seen static and harmonic cases done and for the most part follow the formulation and implementation. However, in the case of the transient problem I am struggling due to the existence of the second order time derivative. Would anybody be able to recommend anything for solving this issue? I was thinking about implementing the Newmark method, however, was unsure how this would translate to the FreeFem code.

In the end I am hoping to set up a small animation but that comes after figuring this out! I am new to FreeFem and relatively new to FEM in general so it is very possible I am missing something basic but any help or guidance would be greatly appreciated.


I’m looking at some electrochemical systems such as polishing hoping to get to
blood coagulation soon. I saw a couple of interesting ways to handle
time - in one case the authors said something IIRC about correct to
second order but I can’t find ti offhand. I looked at just doing FEM along
the time axis which is possible in FF if you only have a 2D model.
I use the final time value as the initial condition for the next iteration
( see my earlier posts ).

You may want to look at that literature the equations are often similar
across disciplines but the names are entirely different :slight_smile:

1 Like

I guess acoustics may have some interesting literature, for example[1]. There is
an issue called pre-echo in audio compression not sure if that is related
to this but generally sampling rate issues may have some similarities.

[1]Dispersion-reducing finite elements for transient acoustics, Yue , Bin and Guddati , Murthy N.; The Journal of the Acoustical Society of America. 2005 2005

Thank you, I’ll be sure to take a look!

piezosquare.edp (3.0 KB)

(sorry I realised I uploaded some incomplete code before, this should be the full version)

Hello again,

So I attempted trying to implement this problem using an easier method, however, when I look at the output it seems that it isn’t behaving quite as expected. I have set the potential at the applied plate to follow a sin curve from 0V to 100V and back before staying at close to 0V. While the simulation does show that the material is expanding when the potential is applied, it does not seem to be relaxing back closer to the original shape after the potential is near 0. I would have expected the elastic term to take over at this point and essentially pull the material back to it’s original shape (save for slight deformation left due to hysteresis) but there is essential no change to the displacement, only getting greater but never less and then maintaining the deformation completely with the (almost) removed potential.

I was wondering if I perhaps set up the variational form incorrectly or did something wrong using the time stepping method that I did.

If anyone has any suggestions on fixing this or using a different method then please let me know. Perhaps I am misunderstanding something really simple.

Thanks again for reading!

I cut your time step in half and the deformation doubled :slight_smile:
You need to look at the formulation again . If you are getitng an incmrental
result the time increment can appear in the numerator.

piezosquare.edp (2.8 KB)

I’m sorry, I don’t think I understand what you mean. When I decrease the time step I still have the same ux displacement?? I certainly feel like something is wrong but can’t quite figure it out, are you saying that the scheme I am using to try and deal with the second order time derivative is formulated incorrectly? I admit I am quite new to FEM in general and am still trying to work out the kinks in my own understanding alongside learning how to implement it in FreeFem so I am always open to learning where I have gone wrong!

I did, however, realise why my output looked like it was only increasing though, I was changing the mesh via:

mesh Th2 = movemesh(Th, [x+(ux)*scl, y+(uy)*scl]);

but have now changed it to account for the relative change in displacement from the previous point, it seems to update now according to the value of ux displayed:

mesh Th2 = movemesh(Th, [x+(ux-uxold)*scl, y+(uy-uyold)*scl]);


It also looks like you updated uxold and uxold2 in the wrong order
if I undestand what you are trying to do but I have not tried to
reconstruct any particular equations but just looking for things
that don’t look right. Certainly your result should not
depend significantly on numerical parameters in most cases.
Checking the dimensions of the terms may help with simple
stuff like that. It may just be the time step is too big to produce
an accucrate result however.

Normally with most
materials iyou expect it to compress
not expand due o the attract between charged surfaces.

I guess if you can post the source equations or a link that would make
it easier to examine.