[Brayanashel]

Philip,

How can stress residual be added beside displacement residual to see it by pyFoamPlotWatcher?

[iblu]

hello
a quick search throughout the OpenFoam threads show some people as me wanting to use OpenFoam to simulate erosion (and particularly erosion rate, rather than "simply" movement of particles in a flow (gaseous, in my case), but basically no many suggestions as to how to proceed.
I guess perhaps OpenFOAM does not have yet the right tools. Am I correct?
Is there something in the pipeline with regards to flow-solid interactions?
I understand ANSYS FLUENT could be used, but there's no opensource module out there (perhaps Impact? )

I am fairly new to CFD but we have a machine to test resistance to erosion by solid particles in an air jet and have been asked to try and model it. 2g/min particles of alumina 50micron in a 5l/min air flow exiting a 1.5mm diameter nozzle and impacting materials at 10mm distance. Supposedly at a speed of 30m/sec (up to 70m/sec). Happy to cooperate if anybody is taking up the challenge

I

[Brayanashel]

How can results of stressFemFoam be seen in ParaView?

[sloti]

Hi everyone,

I was wondering if and/or how it would be possible to simulate the rotation of a solid around a given point with the solidMechanics library.
Imagine for example, you want to simulate a pair of balances. So you have a beam pivoted at some point and two forces acting at the tips, such that the beam can rotate but not translate in order to reach some equilibrium state.
Another application with respect to fluid structure interaction is some wind wheel that starts rotating around a given axis due to fluid flow around it.

Thank you for any ideas!
Regards, Anna

[bigphil]

Originally Posted by sloti

Hi everyone,

I was wondering if and/or how it would be possible to simulate the rotation of a solid around a given point with the solidMechanics library.
Imagine for example, you want to simulate a pair of balances. So you have a beam pivoted at some point and two forces acting at the tips, such that the beam can rotate but not translate in order to reach some equilibrium state.
Another application with respect to fluid structure interaction is some wind wheel that starts rotating around a given axis due to fluid flow around it.

Thank you for any ideas!
Regards, Anna

Hi Anna,

As regards rotation, you should look at these tutorials:

elasticNonLinULSolidFoam/rotateSphereUL
elasticNonLinTLSolidFoam/rotateSphereTL

As regards large rotations within the FSI solver, you would probably need to use GGIs within the fluid mesh.

Best regards,
Philip

[bigphil]

Originally Posted by Brayanashel

How can results of stressFemFoam be seen in ParaView?

Hmnnn… stressFemFoam generates tetPointFields (scalar and vector) and the foam ParaView reader does not recognise these. This problem was discussed previously here but no solution.
Also foamToVTK does not recognise these fields either.

So you could update the ParaView foam reader and/or foamToVTK, or probably the easiest way would be to write a small OpenFOAM utility which reads in these tetPointFields and then transfers their contents to a pointScalarField for example, which ParaView can read.

Philip

[bigphil]

Originally Posted by Brayanashel

Philip,

How can stress residual be added beside displacement residual to see it by pyFoamPlotWatcher?

I am not familiar with pyFoamPlotWatcher, but I presume you would just need to adjust how the residual is written within the solver so as it conforms to how pyFoamPlotWatcher expects the residual to be written.

Philip

[codder]

@ Dr. Cardiff -

I believe there exists a bug in the arbitraryCrack solid model. I have tested 2D trigonal meshes with the elasticAcpSolidFoam solver using meshes generated in both gmsh and ICEM (having initially suspected the error to be associated with a mesh import utility).

But I am now thinking the issue is:

$FOAM_SRC/solidModels/arbitraryCrack/faceCracker/detachFaceCracker.C

As part of this issue, I see under paraview warp-by-vector an unphysical separation of faces not a member of (but adjacent to) the crack patch; conjointly, elements with no opportunity to detach nodes, because no traction will develop on attached faces.

By thumbnail:

#1 - uncracked mesh - crackable boundary patch specified along the five triangular notches
#2 - cracked mesh with labeled displaced problem face (time = 54)
#3 - crack patch in slightly 3-D view, similarly labelled (time = 54)

For example, compilation using "detachFaceCracker_orig.C" worsens the issue.

Best, Eric

[Brayanashel]

Hi,

How can a triangular distributed load be applied on a beam in OF, please?

[misklach]

Hello everyone,
I'm using elasticNonLinULSolidFoam to simulate large displacements of a 2D membrane under a constant load, fixed at the extremities.

With a fixed time-step, I can get grid converged results on modes and amplitude of oscillations, but modifing the time-step I get a different dynamics.

As you can see from the animation (https://drive.google.com/file/d/0B1d...it?usp=sharing), the solution converges for the the final displacement, but the amplitude and the modes of vibration seem highly dependent from the time-step (decreasing the time-step things go even worse).

here the details of the animation
mesh 600x2
deltaT : 0.00125 (blue), 0.0025 (green), 0.005 (red)


Is it a physical or numerical issue? Can anyone explain this behaviour?
Could you link me some reference?

Thank you very much
Giampaolo

[bigphil]

Originally Posted by misklach

Hello everyone,
I'm using elasticNonLinULSolidFoam to simulate large displacements of a 2D membrane under a constant load, fixed at the extremities.

With a fixed time-step, I can get grid converged results on modes and amplitude of oscillations, but modifing the time-step I get a different dynamics.

As you can see from the animation (https://drive.google.com/file/d/0B1d...it?usp=sharing), the solution converges for the the final displacement, but the amplitude and the modes of vibration seem highly dependent from the time-step (decreasing the time-step things go even worse).

here the details of the animation
mesh 600x2
deltaT : 0.00125 (blue), 0.0025 (green), 0.005 (red)


Is it a physical or numerical issue? Can anyone explain this behaviour?
Could you link me some reference?

Hi Giampaolo,

Nice animation.

OK, this time-step effect is a numerical problem.

If you refine the time-step much more (i.e. 0.0005, 0.0001, etc.), do the results converge?

What time-discretisation scheme are you using (e.g. Euler)?
Actually, from a previous post we realised that the 1st Order Euler is currently the only implemented time scheme (backward can be chosen but it just uses Euler for the elasticNonLinULSolidFoam solver by default).

You could actually try the 2nd Order backward time scheme, by editing this line in the solver:

fvm::d2dt2(rho,DU)

and changing it to:

rho*fvm::d2dt2(DU)

and then in fvSchemes, set:

d2dt2(DU) backward;

as dtdt2(rho,DU) is not implemented for backward but dtdt2(DU) is implemented.

Philip

[misklach]

Hi Philip, thank you very much for your quick reply

unfortunately with a smaller time-step the solution doesn't converge (at some point the displacement is not even symmetrical)
https://drive.google.com/file/d/0B1d...it?usp=sharing
dt= 0.00002 (red), 0.00125 (blue)

I'll try if with the backward scheme the dynamic solution converges.

An other quick question: I checked the differences between elasticNonLinULSolidFoam and icoFsiElasticNonLinULSolidFoam, because I remember that in the second one I could use the backward scheme ...basically the two solid solvers should be the same, right? Is there any particular reason why the momentum equation is solved in a different way (d2dt2 does not appear) in the solid solver of icoFsiElasticNonLinULSolidFoam?

[bigphil]

Originally Posted by misklach

Hi Philip, thank you very much for your quick reply

unfortunately with a smaller time-step the solution doesn't converge (at some point the displacement is not even symmetrical)
https://drive.google.com/file/d/0B1d...it?usp=sharing
dt= 0.00002 (red), 0.00125 (blue)

I'll try if with the backward scheme the dynamic solution converges.

OK.

Another thing: the solver does not explicitly include any damping so the membrane should keep oscillating forever; the reason it slows and stops is due to the numerical damping in the 1st order Euler approximation i.e. diffusive error.

Also, what solution tolerance are you using in fvSolution for DU (not the linear solver tolerance)? as this may have an effect given the large number of time-steps: errors may build up.

Originally Posted by misklach

An other quick question: I checked the differences between elasticNonLinULSolidFoam and icoFsiElasticNonLinULSolidFoam, because I remember that in the second one I could use the backward scheme ...basically the two solid solvers should be the same, right? Is there any particular reason why the momentum equation is solved in a different way (d2dt2 does not appear) in the solid solver of icoFsiElasticNonLinULSolidFoam?

The momentum equation is equivalent in both solvers; the d2dt2(DU) term can be written in terms of the difference of ddt(DU) and ddt(DU_0).

Philip

[misklach]

With the backward scheme the membrane keeps oscillating (as expected I guess) ...but then, how can I calculate the steady displacement, given a constant load?

https://drive.google.com/file/d/0B1d...it?usp=sharing
0.005 (green), 0.00125 (red), 0.0003125 (blue)

Also, the solution seems to be in the same range only until a certain time, is it possibly the effect of the large number of time-steps?

In all the simulations I have in the solidMechanics subdictionary
nCorrectors 10000;
DU 1e-07;

Shall I decrease DU tolerance in order to have kind of a converged solution for longer?

[bigphil]

Originally Posted by misklach

With the backward scheme the membrane keeps oscillating (as expected I guess) ...but then, how can I calculate the steady displacement, given a constant load?

OK, that makes sense.

To perform a steady-state simulation, set the d2dt2 time scheme to steadyState in the fvSchemes dictionary:

    //d2dt2(DU)    Euler; // 1st order
    //d2dt2(DU)    backward; // 2nd order
    d2dt2(DU)    steadyState;

steadyState will set the temporal term to zero in the momentum equation; then it is only necessary to perform one time-step (of any size) to get the steady-state solution.

Originally Posted by misklach

Also, the solution seems to be in the same range only until a certain time, is it possibly the effect of the large number of time-steps?

Yes, it is possible: errors can build-up when there are lots of time-steps; but increasing the solution tolerance should help stop them build up.

Originally Posted by misklach

In all the simulations I have in the solidMechanics subdictionary
nCorrectors 10000;
DU 1e-07;

Shall I decrease DU tolerance in order to have kind of a converged solution for longer?

Yep, try 1e-8, 1e-9, etc. and see how it affects the solution.

Philip

[misklach]

I run different cases with fixed time-step and different convergenceTolerance values (1e-9, 1e-12)and apparently there weren't any difference ...too weird.
Then I checked the log files and basically for each case the rel res (which is the convergenceTolerance, right?) jumped from e-7 to 0 making the loop terminate. That means that no matter the value of DU I assigned in solidMechanics all the cases behaved the same.

	Time 0.275, Corrector 182, Solving for DU using DICPCG, res = 1.29921441212e-08, rel res = 2.73569592574e-06, inner iters 2
	Time 0.275, Corrector 183, Solving for DU using DICPCG, res = 2.0893153815e-08, rel res = 2.7124471853e-06, inner iters 2
	Time 0.275, Corrector 184, Solving for DU using DICPCG, res = 2.03314639413e-08, rel res = 2.31764770932e-06, inner iters 1
	Time 0.275, Corrector 185, Solving for DU using DICPCG, res = 1.84189167588e-08, rel res = 7.36534949676e-07, inner iters 1
	Time 0.275, Corrector 186, Solving for DU using DICPCG, res = 7.07846147695e-09, rel res = 0, inner iters 0

Do you know where can I set the precision? Is "writePrecision" in controlDict? I tried to increase that value but without any success

An other question, the pressure at the boundary condition is actual pressure (N/m^2) or pressure/rho?

[bigphil]

Originally Posted by misklach

I run different cases with fixed time-step and different convergenceTolerance values (1e-9, 1e-12)and apparently there weren't any difference ...too weird.
Then I checked the log files and basically for each case the rel res (which is the convergenceTolerance, right?) jumped from e-7 to 0 making the loop terminate. That means that no matter the value of DU I assigned in solidMechanics all the cases behaved the same.

	Time 0.275, Corrector 182, Solving for DU using DICPCG, res = 1.29921441212e-08, rel res = 2.73569592574e-06, inner iters 2
	Time 0.275, Corrector 183, Solving for DU using DICPCG, res = 2.0893153815e-08, rel res = 2.7124471853e-06, inner iters 2
	Time 0.275, Corrector 184, Solving for DU using DICPCG, res = 2.03314639413e-08, rel res = 2.31764770932e-06, inner iters 1
	Time 0.275, Corrector 185, Solving for DU using DICPCG, res = 1.84189167588e-08, rel res = 7.36534949676e-07, inner iters 1
	Time 0.275, Corrector 186, Solving for DU using DICPCG, res = 7.07846147695e-09, rel res = 0, inner iters 0

What is the tolerance for the linear solver in fvSolution (ie for PCG)? The linear solver tolerance must be tighter than the DU tolerance; I typically set the PCG solver tolerance to 1e-12 and relTol as 0.1 (it will never reach 1e-12 due to relTol and DU tolerance but I set it tight to make sure DU tolerance is met).

Originally Posted by misklach

Do you know where can I set the precision? Is "writePrecision" in controlDict? I tried to increase that value but without any success

The numbers inside the solver are stored to double precision, this writePrecision just affects the precision to which they are written to file, so it should not affects your results/convergence.

Originally Posted by misklach

An other question, the pressure at the boundary condition is actual pressure (N/m^2) or pressure/rho?

It is actual pressure (N/m^2 == Pa).

Philip

[misklach]

Thank you very much Philip for your support.

Ok I fixed the issue with the tolerance (the tolerance for the linear solver was bigger than the convergenceTolerance ...shame on me!) ...but apparently it doesn't affect very much the result.

For a time step of 0.005 at time 10 the displacement is still the same. A zoom-in at the middle point shows that the difference is really negligible.
DU convergenceTolerances are 1e-7(blue), 1e-9(red), 1e-12(green)
tolerance2.jpg

tolerance1.jpg

Decreasing the timestep still makes the solution diverge at some point
https://drive.google.com/file/d/0B1d...it?usp=sharing
deltaT = 5e-3(blu), 3.125e-4(red), 7.8125e-5(green)
DU convergenceTolerance = 1e-9

[Brayanashel]

Hi,

In tutorial of elasticSolidFoam, how can analyticalPlateHole be run for different cases just by once compilation, please? In other words, by changing the dimensions of plate we don't have to recompile analyticalPlateHole.C.

[bigphil]

Originally Posted by Brayanashel

Hi,

In tutorial of elasticSolidFoam, how can analyticalPlateHole be run for different cases just by once compilation, please? In other words, by changing the dimensions of plate we don't have to recompile analyticalPlateHole.C.

Hi,

The geometry and loading is hard-coded into the plateHoleSolution function which is implemented in the analyticalPlateHole utility and also within the $FOAM_SRC/solidModels/fvPatchFields/analyticalPlateHoleTraction boundary condition (also in the plateHoleSolution function).
It should be clear to see where the hole radius and applied traction is defined (it should be at the top of the function).
You would need to change these to what ever you want then recompile the utility (with wmake) and the solidModels library (with wmake libso).

Philip