[tschenkel]

Hi,

Yes the new solvers in solids4foam seem to have lower numerical diffusion than the ones in the FluidSolidInteraction toolkit (bazaar). This is good, since I get a better match on the QNET-FSI validation case, which is a more realistic version of the Hron-Turek benchmark, based on experimental results. I get within 10% of frequency on a setup similar to the Hron-Turek benchmark tutorial, where the old FSI toolkit was struggling to get in the same regime.

Unfortunately the new toolkit is a lot less stable (well, obviously, if it has lower diffusion), so the HronTurek benchmark is really tricky to get working.

a) Underrelaxation helps, not only on D (i did not actually need to do that), but in the FSI underrelaxation (I went down by more of an order of magnitude), which is only used for the first three FSI loops, with the IQN method, then IQN takes over. But if the first loops are too far out, it won't converge.

I set max FSI loops to 50, but it should converge within around 20-30, the 50 is just for the unruly timesteps.

You have to check the FSI residuals, so I do a

grep "FSI residual" solids4foam.logfile | tail -n 500 > FSI-residuals

and plot that residual file in gnuplot:

set logscale y

set xlabel "Iteration"
set ylabel "R"
set grid

plot "< tail -n 500 FSI-residuals | cut -d' ' -f7 | tr -d ','" title 'FSI Res' with lines ls 2

This one is nice:

FSI_convergence.png

And this one is not (that's two or three time steps leading up to the problem you describe):

FSI_convergence_borked.png

BTW, there may be an issue with the coupling, I do get weird artefacts in the structure of my QNET case:

weird_artefacts.jpg

Pressures are vertex values, sigma is cell based.

Material model is unsNonLinearGeometryTotalLagrangian, StVenantKirchhoff (rho-1360, E 16e6, nu 0.48).

The weird artefacts seem to be correlated with the mesh size on the fluid side. I can imagine that these could trigger instabilities in the solid solver.

[tschenkel]

Originally Posted by sita

Hi Mikko,

Yeah, I noticed those tutorials, couldn't figure out the difference either though. I still don't quite get how my flexible plate was able to bend straight through the lower wall of my domain, but I'm glad everything is working as it should now.

Thanks again,
Sita

I guess that was due to the fluid mesh not deforming, but the solid solver working on the undeformed mesh. So with the mesh movement inhibited, you basically have a 1-way coupling. So the solid will calculate a steady state deformation based on the current (time dependent) flow. But the solid motion will not influence the vortex shedding (which is still there, since the Re=200 case will have an oscillatory mode even with a fixed blade).

[sita]

Hi Torsten,

Thanks for the tip! Reducing the FSI under-relaxation by a factor of 10 didn't help, even in combination with under-relaxation of D, Deqn, DD, and DDeqn. But I'll play around with the settings a bit further, see if I can get things running.

Regards,
Sita


EDIT: I'm at an FSI relaxation factor of 0.001 now, with under-relaxation of D, Deqn, DD, and DDeqn all at 0.5, but still no luck. I've tried both Aitken and IQNILS, but my run keeps crashing. I didn't think this would be such a hard case, Hron Turek FSI 1 (i.e. steady flow) with gravity in negative y (I wanted to test large deformation). Anything else I can try?

[tschenkel]

Originally Posted by sita

Hi Torsten,

Thanks for the tip! Reducing the FSI under-relaxation by a factor of 10 didn't help, even in combination with under-relaxation of D, Deqn, DD, and DDeqn. But I'll play around with the settings a bit further, see if I can get things running.

Regards,
Sita

BTW, forgot to mention, the errors (or rather warnings about enforcing linear model) you get are no problem, if they go away in the final FSI loops. As long as the final loop is using the correct equations, the results should be fine.

[bigphil]

Originally Posted by sita

EDIT: I'm at an FSI relaxation factor of 0.001 now, with under-relaxation of D, Deqn, DD, and DDeqn all at 0.5, but still no luck. I've tried both Aitken and IQNILS, but my run keeps crashing. I didn't think this would be such a hard case, Hron Turek FSI 1 (i.e. steady flow) with gravity in negative y (I wanted to test large deformation). Anything else I can try?

One comment on using equation relaxation for the solid solvers: I have found the solution to be very sensitive to the value used; generally values less than 0.99 should not be used as the solver tends to prematurely converge (to the wrong answer). Using field relaxation does not have this same sensitivity.

[sita]

Hi Philip,

Thanks, that's good to now. For now I'm afraid premature convergence isn't a problem though; it looks like the solver uses the maximum number of iterations at each time step, without properly converging.

Is there anything else I should try to get my case running properly?

Many thanks,
Sita

[sita]

Hi all,

I'm still struggling to get my case running (see previous posts), it keeps warning me that Material linearity was enforced for stability, and then crashes after a while.

I've been trying to increase the number of iterations, using nIterations in unsNonLinearGeometryTotalLagrangianCoeffs, but it looks like that's not being picked up. What's the proper way to set the (maximum) number of iterations for unsNonLinearGeometryTotalLagrangian?

Also, looking through the code for unsNonLinearGeometryTotalLagrangian, I noticed a parameter stabilisePressure, which is false by default. What's the function of stabilisePressure?

Thanks in advance,
Sita


EDIT: After some more testing, I noticed that the coupling start time makes a difference. If I set it to t = 0 (i.e. coupled = yes), even the original Hron Turek FSI 1 case (steady deformation, no gravity) crashes fairly soon. This isn't the case in the old fsiFoam, does anyone know what's the difference?


Setting the coupling start time to t = 0.1, the case runs alright for g = -2 m/s2, but still not for higher, more earth-like values. Any tips, anyone?

[tschenkel]

Hi,

I'm looking into the HronTurek benchmark solids side of things and have discovered a strange damping.

Setup is as per the tutorial case, but solid only (so CSM2 from the paper, but as a transient case):

NeoHookean Material
E = 5.6e6
nu = 0.4
rho = 1000

unsNonLinearTotalLagrangian


Case description:

Beam of length 0.6m and thickness (height) of 0.02m is fixed at the left side

At t=0 the structure is released in a constant gravity field of (0 -9.81 0).

The simulation runs nicely, but the oscillation is damped.

If I understand correctly there should be no damping in this setup:

- no fluid surrounding the structure, so no external damping
- neoHookeanElastic or StVenantKirchhoffElastic material should be perfect elastic, so no internal damping

So the only damping I can see can come from discretisation errors. I tried a finer mesh, which does improve things very slightly, and going from second order (linear) to third order (cubic):

34c34,35
<     default            Gauss linear;
---
>  //   default            Gauss linear;
>     default            Gauss cubicCorrected;
40,41c41,42
<     laplacian(DD,D)    Gauss linear newSkewCorrected 1;
<     laplacian(DDD,DD)  Gauss linear newSkewCorrected 1;
---
>     laplacian(DD,D)    Gauss cubic newSkewCorrected 1;
>     laplacian(DDD,DD)  Gauss cubic newSkewCorrected 1;

But there is only a slight improvement with refinement and none with order.

Am I missing something?

damped_oscillation.png

[bigphil]

This damping comes from the time discretisation (d2dt2 and ddt).

Euler is only 1st order accurate. You should perform a time-step size sensitivity analysis, and may also try "backward" which is 2nd order accurate (but introduces a lot of dispersion errors).

Either way, Euler with a very small time-step will show minimal time discretisation errors.

Philip

EDIT: FE methods often use trapezoidal or Newmark methods or some combination: these could be implemented

[tschenkel]

Originally Posted by bigphil

This damping comes from the time discretisation (d2dt2 and ddt).

Euler is only 1st order accurate. You should perform a time-step size sensitivity analysis, and may also try "backward" which is 2nd order accurate (but introduces a lot of dispersion errors).

Either way, Euler with a very small time-step will show minimal time discretisation errors.

Philip

EDIT: FE methods often use trapezoidal or Newmark methods or some combination: these could be implemented

Thanks a lot, I am on backward already, so had thought d2t2 would be alright. Time step analysis was next on the list.

As a fluid dynamicist, I have to ask a possibly stupid question: What is the solids equivalent to the Courant number?

EDIT:

It does require quite small time steps, but does work (reduced timestep from 1e-3 to 1e-4s). There is still a little bit of damoing, but it's a lot better:

small_deltat.png

[bigphil]

Originally Posted by tschenkel

As a fluid dynamicist, I have to ask a possibly stupid question: What is the solids equivalent to the Courant number?

Courant number refers to the speed at which information passes through the domain i.e. Courant of 1 means one cell per time step. The Courant number in a solid is based on the speed of sound in the solid (diffusion waves rather than convection waves). For an elastic solid, the speed of sound is approximately sqrt(E/rho) where E is the Young's modulus and rho is the density: actually this is the lower "unconstrained" speed of sound; replace E with (2*mu + lambda) for plane-strain and 3-D, where mu and lambda are Lame parameters.

The speed of sound in solids is typically higher than fluids so the wave speed is high; this means to stick to a Courant number of less than 1 requires very small time steps (especially as you refine the mesh).

However, the stability of implicit solution methodologies (like in unsNonLinGeomTotalLag) does not depend on the Courant number. In contrast, explicit solution methodologies (like in explicitLinGeomSolid) do require the Courant to be less than 1 for stability.

For your time sensitivity analysis, there is probably no need to approach the Courant number, although if you do use such a small time-step then the diffusion error due to the time-step size should be very small.

Philip

[tschenkel]

Originally Posted by bigphil

Courant number refers to the speed at which information passes through the domain i.e. Courant of 1 means one cell per time step. The Courant number in a solid is based on the speed of sound in the solid (diffusion waves rather than convection waves). For an elastic solid, the speed of sound is approximately sqrt(E/rho) where E is the Young's modulus and rho is the density: actually this is the lower "unconstrained" speed of sound; replace E with (2*mu + lambda) for plane-strain and 3-D, where mu and lambda are Lame parameters.

The speed of sound in solids is typically higher than fluids so the wave speed is high; this means to stick to a Courant number of less than 1 requires very small time steps (especially as you refine the mesh).

However, the stability of implicit solution methodologies (like in unsNonLinGeomTotalLag) does not depend on the Courant number. In contrast, explicit solution methodologies (like in explicitLinGeomSolid) do require the Courant to be less than 1 for stability.

For your time sensitivity analysis, there is probably no need to approach the Courant number, although if you do use such a small time-step then the diffusion error due to the time-step size should be very small.

Philip

Thanks for confirmation of the wave propagation speed as the characteristic speed for this case. I'm using implicit methods, so no need for Co=1, but I usually use it as an asymptote for errors.

EDIT:

Quick time step estimate:

SOS = sqrt(E/rho) = sqrt(5.6e-6 Pa / 1000 kg/m^3) = 75 m/s

for Co = 1:

deltat = deltax/V = 0.0033m/75m/s = 4.4e-5s

so we are close with the 1e-4s. Will try 2e-5 and see if that gives close to zero damping.

[tschenkel]

Originally Posted by tschenkel

Quick time step estimate:

SOS = sqrt(E/rho) = sqrt(5.6e-6 Pa / 1000 kg/m^3) = 75 m/s

for Co = 1:

deltat = deltax/V = 0.0033m/75m/s = 4.4e-5s

so we are close with the 1e-4s. Will try 2e-5 and see if that gives close to zero damping.

Update:

There is still a fair amount of damping for a time step of 2e-5 s.

Co<1.png

A lot less than for 1e-3 or 1e-4, but still there.

[bigphil]

Hi Torsten,

The damping error should reduce based on the order of accuracy of the time discretisation: for 1st order Euler, halving the time step will half the error; for 2nd order schemes, halving the time step will reduce the error by a factor of 4.

See the attached image comparing different time schemes for an oscillating beam using a cell-centre finite volume approach (from this thesis).

Philip

[sita]

Hi everyone,

Good, it looks like I'm slowly getting there with my Hron-Turek FSI 1 case with gravity. Delaying the onset of coupling to 0.01 s seemed to help (can't delay it too much though, as around t = 0.24 s the flexible plate would fall out of my domain).

Still, the simulation crashes after just over 0.7 s, and the results start looking funny after some 0.3 s. The tip displacement starts oscillating, as can be seen in the attached plot. Also, at t = 0.4 s waves appear in the plate shape, and the fluid mesh deformation doesn't seem to keep up with the solid mesh deformation (see attached screenshots, one is coloured by p (fluid) and sigma_eq (solid)). With increasing time, the plate shape returns to something looking more physically realistic, and the fluid and solid mesh deformations seem to catch up. But just when I thought things might be getting back to normal and stable, both the interface displacement and the velocity magnitude blow up and the simulation crashes.

This looks like I should be paying more attention to my simulation settings. Is there anything in particular that I should be changing to get stable, physically realistic results?

Many thanks,
Sita


P.S. Adding relaxation only seems to make things worse (both FSI relaxation and D/DD relaxation)

[sita]

Brief update: with RBFMeshMotionSolver, the simulation crashes even more spectacularly, at t = 0.402 s this time (see attached screenshot).

I'm getting a bit desperate, why is this case so much harder in solids4Foam than in fsiFoam (foam-extend-3.1)? Any tips anyone?

Thanks!
Sita

[bigphil]

Hi Sita,

Can you attach a simple example case to reproduce the problem?

Philip

[sita]

Hi Philip,

Sure, here you go. I haven't changed an awful lot since last time I uploaded it (post #496), but clearly the changes that I did make worked out rather disastrously.

Some additional notes:
- With a lower gravity value the case runs just fine
- With a linear elastic material model the case is stable as well, but as then I had to use linear deformation too, the tip displacement in x-direction was negligible, which clearly is not physically correct here.

Thanks in advance for having a look,
Sita

[Daniel_Khazaei]

Originally Posted by sita

Hi Philip,

Sure, here you go. I haven't changed an awful lot since last time I uploaded it (post #496), but clearly the changes that I did make worked out rather disastrously.

Thanks in advance for having a look,
Sita

Hi Sita,

I have looked through the attached case and made three minor changes:

1- file: fsiProperties:

• disabled couplingStartTime and set coupled to yes.

2- file: 0/U

• inlet: change the transition period of transitionalParabolicVelocity to a non-zero value, e.g. 0.2

The case is currently running but it's too soon to judge the outcome.

Regards,
D. Khazaei

[sita]

Hi Daniel,

Thanks a lot for your suggestions! I initially set coupled to yes, and disabled couplingStartTime, just like you did now, but I didn't yet try playing with this transition period of the inlet velocity. With the transition period kept at 0, even the case without gravity (i.e. the original Hron-Turek FSI 1 benchmark) wouldn't run, so I'm definitely going to try playing with this transition period.

Cheers,
Sita