Harmonic Balance Method for nonlinear solid mechanics

[Sami]

Hi all,

I am new to sparselizard but currently I am looking for FEM solution that can yield steady state response to periodic loads for large displacements and rotations (geometric nonlinearity). I read on a forum that sparselizard has an implementation of Harmonic Balance. I am curious if that is accurate information and if it is possible to recreate frequency vs amplitude response for a MEMS duffing oscillator type equation (see fig.) https://en.wikipedia.org/wiki/Duffing_equation.

If yes, then where would be the best place to start wrt creating such FEM models in the sparselizard tool.

Cheers,

Sami

[alexandr...@gmail.com]

Hi Sami,

That's totally true, sparselizard natively supports harmonic balance, and is probably one of the very few FEM softwares (open source or commercial) to support that.

Furthermore the geometric nonlinearity (small-strain but large displacements/rotation) equations are predefined (see the corresponding flavour of predefinedelasticity in the doc or see sparselizard/examples/nonlinear-truss-elasticity-2d/ for an implementation)

and heavily validated on MEMS (vibrating micromembranes, my previous job as a MEMS design engineer for 2 years).

Truth is I had hoped to be able to add exactly such an example online before christmas, but did not have time

--> see attached main.cpp to start with anyways (not validated yet, no guarantees, just a draft I quickly wrote a few weeks ago)

You can easily take any prestress into account in the predefinedelasticity (see doc).

Let me know how that goes,

Alex

[Alexandre Halbach]

Hi Sami,

Well it is how the error message says: it is not allowed to compute the norm of a multiharmonic field as such (I could add it though).

But the workaround is straightforward: compute the norm yourself based on the norm of each harmonic.

The norm of each harmonic, e.g

 harmonic say 5, can be obtained as:

norm(u.harmonic(5))

I will try to find time to allow this for MH fields in the future byt this week is qyite busy.

Alex

On Mon, 8 Mar 2021, 19:36 Sami, <saransh.s...@gmail.com> wrote:

Hi Alex,

I tried running a code to do a frequency sweep as discussed last time. However, I am getting an error with respect to the norm calculation. The idea is as follows :- set the fundamental frequency of the system and a drive frequency which varies in the region of the fundamental frequency to get the softening/hardening response. But the norm calculations throw an error. I've attached the code and the error here. I tried 2 different ways of applying loads (both in the code and referenced here) and I get the same error.

elasticity += integral(load, array1x3(0,0,-2*sin(2*M_PI*fd*t()))*tf(u.harmonic(1)));

    // elasticity += integral(load, array1x3(0,0,-2)*tf(u.harmonic(2)));

Reference for nonlinear loop :

while (relchange > 1.0e-6)

        {

            solve(elasticity);

            double um = norm(u.harmonic(1)).max(vol,2)[0];

            /* Want to calculate norm of total displacement, not just a harmonic term, error thrown that norm can't be computed for multiharmonic field u */

            //double um = norm(u).max(vol,2)[0];

            relchange = std::abs(um-umax)/um;

            umax = um;

            std::cout<<"Max U2: " << um << " ( relative change " << relchange << ")" << std::endl;

        }

Let me know your thoughts on what I am doing incorrectly.

Cheers,

Saransh

--
You received this message because you are subscribed to the Google Groups "sparselizardforum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sparselizardfo...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sparselizardforum/08e130dc-71c8-4c43-b065-c8b8eef4a4f6n%40googlegroups.com.

[Sami]

Hey Alex,

I have the example working but for higher DOFs for some reason the memory exceeds (have 188GB RAM on the server) and the process is killed. I've attached the code for reference and also the snapshot of memory consumption for each nonlinear iteration. Is it a code bug?
As I understand, total # of DOFs being solved are 5553 nodes x 3 x 7 (harmonic coeffs) = 116613 DOFs so not that big an example and the first iteration is fine.

Cheers,

[alexandr...@gmail.com]

Hi Sami,

1. "First iteration is fine": yes it takes at least one nonlinear iteration or more to have all relevant harmonics non-zero so that definitely affects the RAM used and is not surprising

2. "5553 nodes x 3 x 7 (harmonic coeffs) = 116613 DOFs so not that big an example": that's only correct for u.setorder(reg, 1), but you do not use order 1 interpolation but order 3, which has massively more dofs, you are probably closer to 10 million than 100k dofs --> check the exact amount of dofs with:

int numdofs = yourformulationobject.countdofs();

Sparselizard does not use Lagrange shape functions so above order 1 the number of dofs is more than the number of nodes. Shape functions are not associated only to nodes but also to edges and faces and volumes.

--> Maybe use order 2 for u instead of 3 (but not order 1, it sucks for mechanics)

Alex

[Sami]

Hi Alex,

Thanks for the suggestion. With regards to damping, is there a way to use the tangent stiffness matrix in computing the damping (while using rayleigh damping alpha*mass + beta*stiffness)?

Cheers,

[alexandr...@gmail.com]

Hi,

I am not sure to understand the question: if you want to use Rayleigh damping you can simply (for a formulation 'elastic' with a mass M and stiffness K matrix):

mat K = elastic.K();

mat M = elastic.m();

double alpha = 1e-5, beta = 0.1;

mat C = alpha*K + beta*M;

-> C is your damping matrix obtainmed with rayleigh damping.

OR

See https://github.com/halbux/sparselizard/blob/master/examples/eigenvalues-damped-elasticity-membrane-3d/main.cpp line 41 to add terms as usual to the formulation if you do not want to access the C matrix yourself or need to do more advanced stuff.

Alex

[Sami]

Hi Alex,

I referred to the code blob you've linked here while defining damping. My question is if the stiffness proportional term is updated over different nonlinear iterations as the stiffness matrix is updated? Just to understand if the damping I am defining is either linear C = alpha*K(linear, if computed only once) + beta*M or nonlinear C = alpha*K(nonlinear, if computed over every iteration) + beta*M.

Other than that, I've completed the validation example for nonlinear frequency sweep. Let me know if you want to add it to your examples folder. 
Cheers,

[alexandr...@gmail.com]

Hi Sami,

If you use the way using the formulation terms as pointed out in the example then at every time you .generate your formulation the damping matrix will be updated (is that answering your question?).

Looks great! Nice work. Yes I would love to. Once you are all done (is the damping already included?) would it be possible for you to do a code cleaning round (if needed), to put the reference to the paper used for validation at the top and to tell me which characteristic you used for the validation and to how accurate the fit was?

Alex

[Sami]

Hi Alex,

Yes that answers the question. I assumed as much, just wanted to confirm.

I'll need to do some code cleanup and % error in terms of curve fit which I'll do in a week or so. I'll package it in with the literature when I send it your way. Yes, it does include damping (only mass proportional term though). I had a couple of follow up questions with respect to modeling the damping (in context of FSI) :-

1. I saw the FSI example you've included in the examples folder and was curious if I could use multiharmonic analysis in the fluid domain too (to model full fluid damping for my project, if needed)
2. Amongst your solver implementations I saw you also have iterative solvers (e.g. gmres with preconditioning). Are they validated to be used in multiharmonic analysis? I didn't find any examples using iterative solvers so I'm curious where you are with its implementation. If I go the full FSI route, I'll most probably have to switch to iterative solvers since direct solvers take a lot of RAM.

Cheers,
Sami

[alexandr...@gmail.com]

Hi Sami,

Great. Could you as well print out the quantity you validated with as well as the exact value expected, so that I can run this example as well with all other tests for future code verification.

1. Yes. I do not see any limitation for fluids (I actually would like to use MH to simulate  the von Karman vortex example one day). For MH on a fluid geometry deformed by a mechanical displacement you will however need a little more work to bring back calculations on the undeformed mesh (there is an example for that). Future developments should also make this fully automatic, but not yet.

2. Yes, works correctly (like anything else in the documentation should). I do not use it myself because iterative solvers are too often disapointing in their lack of robustness

[Sami]

Hi Alex,

Did you get a chance to take a look at the material I sent over? 

I have a question regarding load ramping while doing a nonlinear iteration. I've written a loop which increments a scaling variable however while invoking the solve function, it doesn't seem to update this. To give a gist of the code structure ->

1. Create formulation "elasticity" and add mass, stiffness and damping contributions using += integral(...)
2. Define the load contribution as : elasticity += integral(domain, 3, array1x3(0,0, scale*(load))*tf(u.harmonic(2))) (initial value of scale is 0)
3. Create nested loop as follows :
        a. Outermost loop -> Change the frequency using setfundamentalfrequency()
        b. Middle loop -> Change the load by updating the value of scale (for e.g. 0.1, 0.2, 0.3 ... 1.0)
        c. Innermost loop -> perform nonlinear iteration where solve(elasticity) is invoked, field u is updated accordingly

The idea is to solve for multiple loads in the same simulation (thus reducing overall number of nonlinear iterations)
Now what happens is the middle loop doesn't do what it's supposed to as it's still solving for the initial definition of scale, i.e. when scale is 0.0 and hence the load contribution is 0. I checked the solve function which invokes the generation function so my understanding is everytime generate function is called, it will create the formulation with updated contributions. Is that correct or am I missing something?

Cheers,
Sami

[alexandr...@gmail.com]

Hi,

When defining an expression with scale (I assume scale is a double) the double value is evaluated and no link to the scale variable is made, so changing scale after adding the formulation term has no effect.

For what you want to do you can ise a parameter: parameters are reevaluated as the are at the evaluation moment.

E.g.: 

parameter scale;

scale|wholedomain = 0.2; // init value

...

yourformulation += integral(... scale...);

// Change scale at any moment:

scale|wholedomain =0.3;

Alex

[Sami]

Hi Alex,

I am trying to understand the nonlinear solver within sparselizard. In this regard, I was wondering if there is a way to compute Jacobian matrix/tangent stiffness matrix (different terms used in different literature) in sparselizard? I would like to update my solution field based on a relaxation factor as follows :-

elasticity.generate();

vec b               = elasticity.b();

mat A               = elasticity.A();

relres             = (b - A*solu).norm()/b.norm();

Kt                   = ......compute Kt..........

solu                = solu + relaxfact*solve(Kt, b - Kt*solu);

u.setdata(selectall(), solu);

where u is updated as follows u = u_old + relaxationfactor*(du_old), du_old = solve(Kt, b - Kt*u_old). In sparselizard terms, Kt would be dA/du_old

[Jeremy Theler]

Hi Sami

Just a side comment: it is usually not recommended to "manually" code a solution scheme when one already has all the PETSc artillery. Did you already check the line search and trus regions algorithms? I don't know the details of your problem, but I bet tweaking PETSc parameters will give better results for sure:

https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESNEWTONLS.html

https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESNEWTONTR.html

[Sami]

Hi Jeremy,

I'll have a look. I am not sure if Sparselizard has SNES implemented, just KSP as far as I know. So might be simpler to play with manually implementing a solution scheme as a first step. Thanks for the info.

[Jeremy Theler]

Hi Sami

Good point! I assumed SNES was there, but now I grepped for "snes" and did not find anything in the source tree.

So we should wait for Alex’ answer then.

Regards

--

jeremy theler

www.seamplex.com

--
You received this message because you are subscribed to the Google Groups "sparselizardforum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sparselizardfo...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/sparselizardforum/1a9618b2-3e28-435c-a168-2dd5b0983fc3n%40googlegroups.com.

[alexandr...@gmail.com]

Hi Sami,

What you are using now is probably a fixed point nonlinear iteration: you reuse the field as it is from the previous iteration and you start over again until convergence. That's totally fine and I would stick to that unless you have a really slow nonlinear convergence.

If you want to add the relaxation factor to your current nonlinear iteration the easiest is to keep the solu (let's call it soluprev) of the previous iteration and when you have your current iteration solu just to:

solu = relaxfactor * solu + (1.0-relaxfactor)*soluprev;

setdata(u);

...

If you really want to use a Newton iteration that's absolutely fine (I provide the predefinedelasticity and predefinednavierstokes as a Newton iteration but it's hidden to the user). For that:

See the post of last week from okazaki3410 on the forum, I provided 2 examples to compute the tangent stiffness matrix for a magnetic saturation problem. If you know the weak formulation of your tangent stiffness matrix then you just add the formulation block for that and the rest uses things that you know, in a nonlinear loop. If you do not have a weak form of your tangent stiffness matrix then it's more difficult of course but I bet you can calculate it by hand in your case.

In the future I will provide that saturation example with a Newton iteration in a clean way online but it's not there yet.

Using petsc as Jeremy mentions is of course one way but a Newton iteration is quite simple to program and one sees exactly what's happening.

Alex

[alexandr...@gmail.com]

Damn: I meant solu is your nonlinear solution 'vec' object:

solu = relaxfactor * solu + (1.0-relaxfactor)*soluprev;

setdata(solu);

[Jeremy Theler]

On Mon, 2021-05-03 at 05:42 -0700, alexandr...@gmail.com wrote:

Using petsc as Jeremy mentions is of course one way but a Newton iteration is quite simple to program and one sees exactly what's happening.

That's my point. Plain-old Newton is way too simple. PETSc gives a plethora of state-of-the-art methods which beat any other handcrafted approach. Actually the first sentence of all PETSc tutorials is: "please use all the solvers we provide, not just the linear solvers and code your own SNES & TS."

--

jeremy theler

www.seamplex.com

[Alexandre Halbach]

I am sure their solvers are great of course but for all the cases until now a fixed point or Newton iteration is more than enough. I could probably go more brutal in the FSI example or in the bending pillar with geometric nonlinearity but that was until now never needed and there are other things to develop first :)

[Sami]

Ok Thanks, I'll have a look. So far fixed point hasn't worked well. Even with load ramping, the solution is oscillating!

[alexandr...@gmail.com]

Did you make it more nonlinear compared to before when it was working?

BTW: is it ok if I mention your name on the example I put online?

Alex

[Sami]

The one before was a different example, a smaller model so to speak. The one I am modelling right now is closer to the actual MEMS mirror design we fabricate which is rather complex with large number of DOFs and actually quite nonlinear. The idea is to reduce the nonlinearity in design phase of our next gen mirrors.

Yeah sure. I believe you have the details from my company account.

[Alexandre Halbach]

Let me know how the relaxation has helped. It's the easiest immediate potential solution.

I just ran your cpamped clamped example, I could speedup a bit the generation in the future: it should be very fast but with the elasticity with geometric nonlinearity it becomes quite heavy of course in 3D.

Alex

--
You received this message because you are subscribed to the Google Groups "sparselizardforum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sparselizardfo...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/sparselizardforum/5c13b17d-9070-4580-a172-d63291f60d34n%40googlegroups.com.

[Alexandre Halbach]

In any case the DDM framework is in the finalizing phase. That should allow to speedup the calculation nicely in a month or so.

Alex

[alexandr...@gmail.com]

Sami / Jeremy:

Sparselizard matrix storage has changed (and will not change anymore) 3 weeks ago. The Dirichlet constraints are now not stored in the mat object any more but instead A is now in the form:

A = [ Aa   Ad]

       [  0       1]

i.e. only As and Ad are stored (see documentation at mat::getainds).

If you want to use SNES and do anything by yourself in petsc you can still use .getpetsc() for vec objects but now you need to do .getapetsc and .getdpetsc to get Aa and Ad petsc objects (see doc).

You might want to use functions (soon in the doc):

vec bwithoutdirichlet = A.eliminate(b);

Vec bpetsc = bwithoutdirichlet.getpetsc();

Mat Apetsc = A.getapetsc();

so that you do not care about the d part anymore and just have to solve the Aa*xa = ba part.

[Sami]

Hi Alex,

I am currently debugging the snes implementation and discovered that there might be an issue with ksp instead. When I use fixed point iteration without snes, the convergence is quick.

But when I switch to SNES and try the same approach, i.e. retrieve ksp from snes and set the pc etc. I get the convergence of linear solve to be quite bad. Have you ever encountered this problem before?

I've attached the sl.cpp file for your reference. Look at the function solvenonlinear.

[alexandr...@gmail.com]

Hi,

When you say there might be an issue with the ksp: do you mean the ksp as created in sparselizard?

The convergence in 2 iterations for the truss nonlinear simulation seems to fast to be correct, unless you have a very small load and are almost linear? Did you check if the solution was correct?

For the snes part I have never used it so maybe petsc team will be best to help for general questions.

Alex

[Sami]

Yes, I reduced the load for faster convergence. I started out with your code verification to make sure everything is correct. 

Just for curiosity, when I call formul.generate() and then .A() and .b(), does it flush the matrices/vectors from the memory while copying. If so, I might have an idea. I am calling .generate() function in both the main solvenonlinear function and myresidual function, used to calculate the residual. It might be that the 2 are not related anymore. 

[Jeremy Theler]

The usual way to try to catch problems in SNES is to ask PETSc to compute a numerical jacobian of the residual with finite differences with -snes_test_jacobian -snes_test_jacobian_view 

See https://www.mcs.anl.gov/petsc/documentation/faq.html#newton

--

jeremy theler

www.seamplex.com

[alexandr...@gmail.com]

Hi Sami,

Yes, when you .A() or .b() then the A and b live their own life and everything related is flushed from the formulation. If you want to reuse the previously generated fragments you can use .A(true) or .b(true) (see doc)

Alex

[Jaka Pribosek]

Hi all,

I am a first time user of Sparselizard and really impressed by its capabilities.

Over last couple of days, I am trying to run Sami's example provided  here and got somewhat different result from the one used in validation example.  The response (blue curve) did not reach the right amplitude nor foldover frequency (here 1035rad/s) on the scan up.

I tried increasing the number of harmonics to 7  and while it helped to somewhat increase the fold over frequency  (1070rad/s ), the response got too slanted, fell over too quickly and did not reach the expected amplitude (should reach value beyond 0.02m). However, scan-down works fine with only very little difference  between 5 or 7 harmonics and seems to be consistent with the validation example  (amplitude jump at 1025rad/s).

I am running these examples with coarser mesh of the both provided in the example and 5 rad/s steps. Could my problem stem from the fact that mesh is too coarse? I tried running it with finer mesh but terminated the simulation after three or so hours as it did not find solution for elasticity on first iteration.

Any help or hint in the right direciton would be much appreciated.

 

sreda, 26. maj 2021 ob 12:32:22 UTC+2 je oseba alexandr...@gmail.com napisala:

[alexandr...@gmail.com]

Hi Jaka,

Thanks for your user feedback and question. (Just for info all future examples will be provided at https://github.com/halbux/sparselizard-users and the backbone example is there as well).

@Sami is of course the best person to answer this but I can try myself:

The value that the curve reaches is clearly impacted by the attenuation factor used on line 55:

// Mass damping only, C = alpha*M

double alpha = 3.0;

You can surely get any amplitude by changing the attenuation. I think Sami was still working on that. I don't know if they provide a proportional damping set of parameters in the papers but having the right loss factors is very needed for an accurate simulation.

Regarding the finer mesh: it surely affects the results but I understand 3D mechanics with geometric nonlinearity and 7 harmonics can take a while to solve since it leads to a quite large and dense problem (possibly all harmonics are coupled together). There

are ways to speed that up though.

Hope this helps,

Alex

[Sami]

Hi Jaka,

Could you attach the modifications you made (the .cpp file) here? With both coarse mesh and fine mesh, you should get similar results, time taken obviously depends on your system specs. The solution also depends on what convergence criteria you provide in the nonlinear iteration loop (along with any physical parameters such as damping). Just out of curiosity, where are you measuring the amplitude? It should be measured at the center of cross-section area of the beam, in the middle of the beam.

[Jaka Pribosek]

Thank you both for your inputs and prompt response!

Sami, your question about the amplitude got me thinking, I am not plotting the right ampltiude here. Instead, I am plotting maximum of the second harmonic term ( norm(u.harmonic(2)).max(vol,5)[0];) which is wrong. I noted from your previous discussion, that taking the norm of the multiharmonic field is not (yet) possible in sparselizard. Is there other convenient way to do it ?  Sami, how did you measure the amplitude?

ponedeljek, 14. junij 2021 ob 13:25:54 UTC+2 je oseba Sami napisala:

[alexandr...@gmail.com]

Hi Jaka,

It is not yet possible with a one liner but you can easily write a 10 liner function for that since it all boils down to taking the norm of a sum of sines and cosines at frequencies multiples of the fundamental f0. And every harmonic can be accessed by .harmonic()

Alex

[Jaka Pribosek]

Thank you for the hint. Now it works!

torek, 15. junij 2021 ob 08:05:19 UTC+2 je oseba alexandr...@gmail.com napisala: