[tschenkel]

Hi,


I have moved my code to OF1912 (easier for students and have implemented another model that was difficult to port to foam-extend).

Now I realised that one of the solids4Foam BCs I use (fixedDisplacementZeroShear) hasn't been ported yet. Are you currently working on that?


I could give it a try, but while my local branch does compile the BC now, it does not behave nicely, so I must have missed a different field definition. My OF developing expertise is not on the "implement from scratch" level. I'm more a "copy and adapt" hacker, in order to get a new model equation in, so maybe out of my depth with that.


Do you maybe have a hint to what the differences in BC definition are between fe and (ESI)-OF? On the fluid side of things I'd be more comfortable, but since solids4Foam is in development software I'm struggling to see the wood for the trees.

[bigphil]

Hi Torsten,

I have now ported fixedDisplacementZeroShear to compile with OpenFOAM-v1912; there were only a couple of minor changes. It compiles but I have not double-checked its behaviour. I pushed the commit to the development branch.

On a side, solid tutorials seem to run much slower in OpenFOAM-v1912 than v1812; I thought this was related to the dataTemplates fix but it seems to be something else (probably related to storing residuals). Have you noticed the same behaviour?

Philip

[Daniel_Khazaei]

Originally Posted by bigphil

Hi Torsten,

I have now ported fixedDisplacementZeroShear to compile with OpenFOAM-v1912; there were only a couple of minor changes. It compiles but I have not double-checked its behaviour. I pushed the commit to the development branch.

On a side, solid tutorials seem to run much slower in OpenFOAM-v1912 than v1812; I thought this was related to the dataTemplates fix but it seems to be something else (probably related to storing residuals). Have you noticed the same behaviour?

Philip

Well I have reported this awhile ago by email, it seems that the implicit loop of the solid models requires more iterations to converge when running the cases on OpenFOAM distributions. If I remember correctly, I first noticed the difference when working on HronTurekFsi3 test case.

[tschenkel]

Originally Posted by bigphil

Hi Torsten,

I have now ported fixedDisplacementZeroShear to compile with OpenFOAM-v1912; there were only a couple of minor changes. It compiles but I have not double-checked its behaviour. I pushed the commit to the development branch.

On a side, solid tutorials seem to run much slower in OpenFOAM-v1912 than v1812; I thought this was related to the dataTemplates fix but it seems to be something else (probably related to storing residuals). Have you noticed the same behaviour?

Philip

I haven't run comparisons between 1812 and 1912 - and also no solids only cases (I know I should).

The one thing that I did notice is that RBF is a lot slower than AMI/GCI on the interpolation.

[tschenkel]

Edit: The problem below was only an instability, I had accidentally copied the wrong timestep over. It seems to work now.

Originally Posted by bigphil

Hi Torsten,

I have now ported fixedDisplacementZeroShear to compile with OpenFOAM-v1912; there were only a couple of minor changes. It compiles but I have not double-checked its behaviour. I pushed the commit to the development branch.
Philip

I have just pulled the updated files and your changes are exactly the same ones I made.


So I'll need to run a few more tests, but the first tests show a weird behaviour that wasn't there in foam-extend. See attached image. It seems as if there is only deformation on the block boundaries (4 block blockMesh on the solid). Case is essentially the 3dTube case in full 3d (D=20mm, T=2mm, L=100mm) with a linear elastic solid E=1e6, nu=0.3, rho=1000. Newtonian luid nu=4e-3, rho=1000. Inlet pressure BC jumps to 5 (5000Pa) at t=0. Based on:

Lee, S.H., Choi, H.G., Yool, J.Y., 2012. Finite element simulation of blood flow in a flexible carotid artery bifurcation. J Mech Sci Technol 26, 1355–1361. https://doi.org/10.1007/s12206-012-0331-9



On a side note: with fixedDisplacementZeroShear on D, what would be the appropriate BC for pointD and pointDD?
At the moment I run those at "calculated", because non of the others made sense.


zeroShearProblem.jpg

[tschenkel]

I just checked the error log on the cluster (I didn't check that before since the case started ok), but there is an error on the v1912 run:



Duplicate entry skewCorrected in runtime selection table snGradScheme │Reconstructing point fields


This was moved from newSkewCorrected to skewCorrected recently. Could it be that that clashes with another scheme of the same name in OF1912?

[tschenkel]

Originally Posted by bigphil

Hi Torsten,

I have now ported fixedDisplacementZeroShear to compile with OpenFOAM-v1912; there were only a couple of minor changes. It compiles but I have not double-checked its behaviour. I pushed the commit to the development branch.

Philip

Your changes were exactly as mine, and the problems I had were stability problems caused by me not being careful while porting my case to OF.

It does work now and the BC behaves as it did in FE41.

However, there are some differences between FE and OF (v1812 at the moment, v1912 is still compiling on the cluster) and also between the results of solids4Foam and the benchmark results I'm using (may have to compare those, but it's three different solvers against solids4Foam).

The case is from (adapted from the 3dTube tutorial):

Lee, S.H., Choi, H.G., Yool, J.Y., 2012. Finite element simulation of blood flow in a flexible carotid artery bifurcation. J Mech Sci Technol 26, 1355–1361. https://doi.org/10.1007/s12206-012-0331-9

Meshes are blockMesh (butterfly / O-Grid for the fluid, 4 quadrants aligned with the fluid mesh for the solid), with cell counts (length refinement ratio is 1:2):

Solid fine: 69120
Solid coarse: 8640

Fluid fine: 201600
Fluid coarse: 25200

FSI interpolation: FE - GCI, OF, RBF (non-matching meshes), direct (matching meshes)

Solid: rho = 1000, E=1e6, nu=0.3, linear Elastic
Fluid: icoFoam, rho = 1000, nu = 4e-3, Newtonian, laminar


deltaT = 1e-5

IQNLS



Figure 3 in Lee's paper shows the deformation at different times:


lee2012_Fig3.png



And these are the results from solids4Foam (refinement labels are fluid/solid, i.e. coarse/fine is coarse fluid and fine solid mesh - fine fluid mesh simulations still running):


Lee_Validation.jpg


As you can see the shape of the deformation is similar, but the amplitude is about 25% lower. It also drops below the z=0 deformation with time, where the published data shows the main wave travelling on without loosing in amplitude. There is also a rather large difference between FE and OF.

I had hoped this would be a straightforward benchmarking exercise that I could share with you, but I seem to be missing something.

[bigphil]

Did you check the effect if time-step size?

I wonder do all the methods approach the same solution as the mesh and time-step sizes are reduced.

This case was examined in the paper so it should work, at least in fe40.

Philip

[tschenkel]

Originally Posted by bigphil

Did you check the effect if time-step size?

I wonder do all the methods approach the same solution as the mesh and time-step sizes are reduced.

This case was examined in the paper so it should work, at least in fe40.

Philip

Thanks for the paper. I noticed that the number of cells for the 3dTube case is very high for both fluid and solid in that paper. So it seems that the diffusive nature of the FVM does indeed require a lot more cells than FEA would require elements (or so my local FEA expert tells me).

Coupling stability also seems to be very sensitive to the time-step. It works for 1e-5, but going to either 1e-4 or 5e-5 becomes unstable.

My case (Lee et al.) seems to be less benign than yours. Differences are:

- both ends are not fixed, but allowed to move in plane: BCs D: fixedDisplacementZeroShear, poindD[D]: calculated

- diameter 20mm, wall thickness 2mm, E=1e6, nu=0.3, rho=1000, fluid: nu=4e-3, rho=1000 (so density ratio a little less favourable than yours).

- pressure jumps to 5000Pa at t=0, and then stays there

The PWV seems to come out nicely, haven't done a full evaluation of that yet, but the deformation seems to be a lot more sensitive to diffusion (both temporal and spatial?).

I have tested the case on OF and CalculiX, coupled with preCICE. That gives me the right deformation, but wave speed is a lot lower.


I'll send your paper to my MSc student ;-)

[bigphil]

Hi again Torsten,

Originally Posted by tschenkel

Thanks for the paper. I noticed that the number of cells for the 3dTube case is very high for both fluid and solid in that paper. So it seems that the diffusive nature of the FVM does indeed require a lot more cells than FEA would require elements (or so my local FEA expert tells me).

The element number is high in the paper to ensure small spatial discretisation error, however, I politely do not agree with the statement "diffusive nature of the FVM does indeed require a lot more cells than FEA would require elements": this is not my experience
There are two main contributions to diffusion error:

• spatial discretisation: there are many variants of FE and FV discretisation, each with their own accuracy and order-of-accuracy. For linear local variations (e.g. linear shape functions), FV approaches can often display 2nd order accuracy (on average) in displacement and stress/strain, with FE 2nd in displacement and 1st in stress/strain (of course, there are many exceptions). For a concrete example of comparing FE/FV accuracy and order-of-accuracy, see this paper,
  https://www.sciencedirect.com/scienc...45794916306046; in particular the bending beam problem (example 1) which shows two FV formulations (both in solids4foam), one which is more accurate than Abaqus reduced integration linear elements, and one which is less accurate.
• temporal discretisation: for problems where inertia is important, you will need to use a 2nd order temporal scheme to reduce energy diffusion or a 1st order method with "small" time-steps. In solids4foam/OpenFOAM, we have 1st order Euler which requires small time-steps to minimise diffusion (it is actually not very efficient for wave problems); there is also the 2nd order "backward" scheme, which can be OK but also tends to be quite unstable for large time-steps. If a 2nd order Newmark-beta scheme was implemented (common in FE and works quite well), then you could expect very similar behaviour to common FE procedures.

There is also stabilisation, which can affect accuracy (as it tends to be a diffusion term); some FV formulations require it (e.g. Rhie-Chow) like some FE schemes (hourglass control).

So if in your case, your are finding that you have premature wave dissipation then I would guess it is probably related to the time-scheme; ideally Newmark-beta would be implemented in solids4foam (it has been on my list for a while!) but failing that use a smaller time-step with Euler or backward.

Originally Posted by tschenkel

Coupling stability also seems to be very sensitive to the time-step. It works for 1e-5, but going to either 1e-4 or 5e-5 becomes unstable.

Yes, this is true. I believe this is partly related to the definition of the residual, which currently depends only on what happens in a time-step and not on the accumulated error (it would be interesting to examine other residual definitions). Also, it may be related to the implementation of the coupling schemes (e.g. Aitkens/IQN-ILS implementations can vary); I suggest switching between Aitkens and IQN-ILS to see if it helps.

Originally Posted by tschenkel

My case (Lee et al.) seems to be less benign than yours. Differences are:
- both ends are not fixed, but allowed to move in plane: BCs D: fixedDisplacementZeroShear, poindD[D]: calculated

- diameter 20mm, wall thickness 2mm, E=1e6, nu=0.3, rho=1000, fluid: nu=4e-3, rho=1000 (so density ratio a little less favourable than yours).

- pressure jumps to 5000Pa at t=0, and then stays there

The PWV seems to come out nicely, haven't done a full evaluation of that yet, but the deformation seems to be a lot more sensitive to diffusion (both temporal and spatial?).

Yes it does seem less benign

If both ends are "fixedDisplacementZeroShear", you may have an issue with rigid body translation, with the pipe "floating away"; this would not be good for convergence.

Regarding diffusion, see comments above.

Originally Posted by tschenkel

I have tested the case on OF and CalculiX, coupled with preCICE. That gives me the right deformation, but wave speed is a lot lower.

I'll send your paper to my MSc student ;-)

Interesting, I wonder why the wave speed is lower with OF/CalculiX, do you have an idea? I guess a mesh/time-step sensitivity analysis would also be useful.

Philip

[sita]

Hi all,

I'm currently working on a simulation in which I apply a linearly increasing pressure to the interior of a geometry, similar to the cylindricalPressureVesselBijelonja tutorial. This works perfectly fine as long as I'm using some non-linear elastic model, like solidModel nonLinearGeometryUpdatedLagrangian in solidProperties, and type neoHookeanElastic in mechanicalProperties.

However, when I switch to a linear elastic model, like solidModel unsLinearGeometry in solidProperties, and type linearElastic in mechanicalProperties, sigmaHyd is giving me all sorts of errors (the solver seems to expect scalar values in the sigmaHyd boundary conditions, and grad(sigmaHyd) isn't a valid volVectorField).

What would be the correct way of applying a linearly increasing pressure to the interior of a geometry if a linear elastic model is used?

Many thanks!
Sita

[bigphil]

Hi Sita,

I have just tried this myself and it seems to work as expected for both linearGeometry and unsLinearGeometry.

For applying a linearly increasing surface pressure, this should be applied in the solidTraction boundary condition for that patch using a pressureSeries; is this what you are doing? For example, it is not correct to apply a linearly increasing boundary condition to the sigmaHyd field as sigmaHyd only represents the hydrostatic component of the stress, whereas solidTraction will enforce the total surface pressure to the requested value.

Also sigmaHyd is only needed if solvePressureEqn is yes in mechanicalProperties (it defaults to no if not found).
Philip

[sita]

Hi Philip,

Thanks a lot, you're fast!

I used the cylindricalPressureVesselBijelonja tutorial as an example for my simulation, so yes, I'm indeed applying the time dependent pressure as a pressure series in the solidTraction boundary condition.

I just noticed that in this tutorial solvePressureEqn is set to yes, so I'll try running my case with solvePressureEqn set to no, see if that helps.

Just out of curiosity: when would you typically need to set solvePressureEqn to yes? And why am I getting these errors only in case of a linear elastic material?

Thanks,
Sita

[sita]

Oh, too bad, turns out I didn't copy this "solvePressureEqn yes" bit from the tutorial case, so that's not the cause of the errors I'm getting.

Leaving out the sigmaHyd boundary condition did help in that it eliminated the error messages, but without it my geometry just doesn't deform at all.

Think I'll have to experiment a bit more. Apart from the cylindricalPressureVesselBijelonja tutorial, are there other (simpler) tutorials that I could check?

[bigphil]

Hi sita,

You could look at the following tutorial which uses a time-varying surface pressure (pressureSeries): solids4foam-release/tutorials/solids/linearElasticity/pressurisedCylinder.

As regards when to set solvePressureEqn to yes, this is typically requires (or improves accuracy and convergence speed) when your material is incompressible or close to incompressible e.g.a Poisson's ratio close to 0.5. It is also useful if you have large plastic deformation and your plasticity model assumes incompressible deformation.

Philip

[sita]

Hi Philip,

Thanks, that looks like exactly the tutorial I was looking for. Unfortunately it's not among the tutorials I have. Do you think I should update my S4F version (my current version was installed last July)?

Thanks for the additional explanation too; slowly but surely I'm starting to actually understand what I'm doing

Cheers,
Sita

P.S. Not sure whether different versions exist, but I'm using the Git-repository from Bitbucket, and I'm currently on the master branch (up to date)



P.P.S. The tutorials that I do have in solids/linearElasticity are:
biMatPlate oscillatingPlate rigidPunch
biMatPlateSkew patchTest slidingFrictionBall
cantilever2d plateHole tetBeam
ellipticPlate plateWithFixedInternalCells uniaxialTensionRod
multiMaterialBallBrick punch waveBar
narrowTmember rigidCylinderContactBrick

[bigphil]

Apologies, the pressurisedCylinder case was not pushed yet (I need to check why). Either way, I will attach it here.

[sita]

Originally Posted by bigphil

Apologies, the pressurisedCylinder case was not pushed yet (I need to check why). Either way, I will attach it here.

Ah, that explains. Great, thanks a mill!

[sita]

Hi Philip,

Sorry, I think I didn't quite understand what was happening earlier, when I wrote that after removing sigmaHyd from the boundary conditions the geometry didn't deform. Am I correct in assuming that, for linear elastic simulations, S4F doesn't compute/store the deformed mesh?

In the simulation I'm trying to do, the Young's modulus is low, so deformation can be considerable, even if the material is linear elastic (i.e. linear elastic material, non-linear deformation). Is there a way to have S4F compute/store the deformed mesh, like it does for non-linear materials?

Many thanks,
Sita


P.S. At the moment I'm using solidModel unsNonLinearGeometryTotalLagrangian; combined with a StVenantKirchhoff type material

[bigphil]

Originally Posted by sita

Sorry, I think I didn't quite understand what was happening earlier, when I wrote that after removing sigmaHyd from the boundary conditions the geometry didn't deform. Am I correct in assuming that, for linear elastic simulations, S4F doesn't compute/store the deformed mesh?

Yes you are correct.

For linear geometry (aka small strain) solid models and total Lagrangian nonlinear geometry (finite/large strain formulation) solid models, the deformed mesh is not stored or written. However, you can see it in Paraview by using the "Warp by vector" option, selecting the total displacement field (D or pointD) and setting the scale factor to 1.

Originally Posted by sita

In the simulation I'm trying to do, the Young's modulus is low, so deformation can be considerable, even if the material is linear elastic (i.e. linear elastic material, non-linear deformation). Is there a way to have S4F compute/store the deformed mesh, like it does for non-linear materials?

Can I ask why you need to store the deformed mesh? Maybe "warp by vector" is all you need?

Originally Posted by sita

P.S. At the moment I'm using solidModel unsNonLinearGeometryTotalLagrangian; combined with a StVenantKirchhoff type material

That's fine, however, the StVenantKirchhoff law is not really intended for truly large strains, it is really just a simple extension of a Hookean (linear elastic) law to deal with finite rotations. Instead I suggest using "neoHookeanElastic" as it gives physically realistic results for larger strains, and it will be more stable too.