[Barry]

Hi (Structural) Foamers,

Thought this post might help some of you to understand the importance of time step in a structural solver using solids4foam toolbox.

The frequency of cantilever beam was validated using the finite volume method by the authors Slone et al https://www.sciencedirect.com/scienc...07904X02000604. Please see fig 3 in this paper. The natural frequency is approx 20 s.

The same case was run using linearGeometryTotalDisplacement model from solids4foam toolbox. The left part of the attached image shows the displacement variation with time for the force inputs as shown in right for different time steps.

As given in the paper, the force was linearly ramped upto 10 s and then released. However, the displacement from simulation does not capture the natural frequency when time step, dt = 1 s. There is too much damping arising from first order Euler scheme. With finer time steps, this numerical damping disappeared to a large extent and the natural frequency is as expected. Unlike the Slone et al's article, the peak displacement is greater than 0.1 m and it gradually reduced with increasing time. Thanks to Philip Cardiff (aka bigphil) for suggesting to perform time step checks.

[Barry]

I forgot to mention in the above message. At time step = 0.1 s, the results for backward scheme and Euler scheme are same though. As backward scheme is second order and Euler is first order, why is displacement response not sensitive to the type of time scheme?

[bigphil]

Originally Posted by Barry

I forgot to mention in the above message. At time step = 0.1 s, the results for backward scheme and Euler scheme are same though. As backward scheme is second order and Euler is first order, why is displacement response not sensitive to the type of time scheme?

As the time-step is reduced, all consistent methods will approach the same answer. The difference is that second-order methods will approach the time-step-independent answer at a quicker rate the first-order methods (well, actually only in the asymptotic region i.e. when getting close to the answer).

Also, for a given time-step size, second-order methods are often more accurate in absolute terms than first-order methods, but this does not have to be the case.

In your case, I suggest you perform the same time-step analysis using both Euler and backward and compare them, and you should see that backward approaches the final solution at a higher rate.

Philip

EDIT: by the way, first order Euler is generally known to require very small toe-steps to be accurate in comparison to Newmark or trapezoidal second order schemes, which can produce the same accuracy with much larger time-step sizes.

[chaosCFD]

Originally Posted by Barry

Hi (Structural) Foamers,

Thought this post might help some of you to understand the importance of time step in a structural solver using solids4foam toolbox.

The frequency of cantilever beam was validated using the finite volume method by the authors Slone et al https://www.sciencedirect.com/scienc...07904X02000604. Please see fig 3 in this paper. The natural frequency is approx 20 s.

The same case was run using linearGeometryTotalDisplacement model from solids4foam toolbox. The left part of the attached image shows the displacement variation with time for the force inputs as shown in right for different time steps.

As given in the paper, the force was linearly ramped upto 10 s and then released. However, the displacement from simulation does not capture the natural frequency when time step, dt = 1 s. There is too much damping arising from first order Euler scheme. With finer time steps, this numerical damping disappeared to a large extent and the natural frequency is as expected. Unlike the Slone et al's article, the peak displacement is greater than 0.1 m and it gradually reduced with increasing time. Thanks to Philip Cardiff (aka bigphil) for suggesting to perform time step checks.

Hello Barry, thanks for sharing! May I ask a question? How is the calculation time of your case for the cantilever? Is the calculation slow? Recently, I also conducted a free vibration simulation of a cantilever using solids4Foam, but the solving procedure is very slow. A 4s simulation lasts nearly 6 hours. The total cell number of the solid model is only 2000, and the timeStep is set as 0.01. How about your case?

[bigphil]

Originally Posted by chaosCFD

Hello Barry, thanks for sharing! May I ask a question? How is the calculation time of your case for the cantilever? Is the calculation slow? Recently, I also conducted a free vibration simulation of a cantilever using solids4Foam, but the solving procedure is very slow. A 4s simulation lasts nearly 6 hours. The total cell number of the solid model is only 2000, and the timeStep is set as 0.01. How about your case?

The segregated solution algorithm is slow for high-aspect-ratio geometry, particularly for slender cantilevers. You may find coupled solvers to be a better option in your case, e.g. I showed a coupled solver to be much faster than a segregated solver in cases like this in this paper. In solids4foam, you can try coupledUnsLinGeomLinearElasticSolid and vertexCentredLinGeomSolid.

[chaosCFD]

Originally Posted by bigphil

The segregated solution algorithm is slow for high-aspect-ratio geometry, particularly for slender cantilevers. You may find coupled solvers to be a better option in your case, e.g. I showed a coupled solver to be much faster than a segregated solver in cases like this in this paper. In solids4foam, you can try coupledUnsLinGeomLinearElasticSolid and vertexCentredLinGeomSolid.

Thanks for your suggestion, Prof. Cardiff! Yeah, I also find the coupled solver is much faster, but what is the mechanism behind it? Why the segregated method is slow for high-aspect-ratio structure and why the coupled is faster? Yesterday, I ran a case of which a cantilever subjected by a sudden load with segregated method. I set a very small time step, about 1/5000 of the cantilever's nature period, and the solution was much faster, because the iteration steps became fewer in every calculation time. By the way, I couldn't find the paper with the url. Is this the paper? A block-coupled finite volume methodology for linear elasticity and unstructured meshs.

[bigphil]

Originally Posted by chaosCFD

Thanks for your suggestion, Prof. Cardiff! Yeah, I also find the coupled solver is much faster, but what is the mechanism behind it? Why the segregated method is slow for high-aspect-ratio structure and why the coupled is faster? Yesterday, I ran a case of which a cantilever subjected by a sudden load with segregated method. I set a very small time step, about 1/5000 of the cantilever's nature period, and the solution was much faster, because the iteration steps became fewer in every calculation time. By the way, I couldn't find the paper with the url. Is this the paper? A block-coupled finite volume methodology for linear elasticity and unstructured meshs.

Sorry, I have now fixed the URL to the paper; here it is again: P. Cardiff, Ž. Tuković, H. Jasak, A. Ivanković, A block-coupled Finite Volume methodology for linear elasticity and unstructured meshes, Computers & Structures, 175, 2016, https://doi.org/10.1016/j.compstruc.2016.07.004.


In reply to your questions

• Why is the segregated solver slow (relative to a coupled solver) for these cases? The segregated solver uses an approximation of the divergence of stress term in the momentum equation for the implicit term (a Laplacian term is used); then explicit (deferred corrections) are added so that once the outer iterations converge then the solution satisfies the original momentum equation (with the divergence of stress). This approach is efficient for "blocky" geometry with low aspect ratios; however, for bending of slender cantilevers, explicit terms are large relative to the implicit term resulting in slow convergence of the outer iterations. A coupled method treats all (or most) of the divergence of stress term implicitly, so it does not have this problem. Another solution could be to use a geometric multigrid algorithm, as done by Demirdzic, Muzaferija, Peric and co-workers.
  
• Why is the segregated solver much faster when solving the case transiently using small time steps?
  When using small time steps in transient problems, the inertia term (rho*dv/dt) is large relative to the stress term's divergence, so convergence is fast. In these cases, the coupled method could be slower because it is more expensive to build and solve the coupled linear system than the segregated system. For really small time steps, fully explicit methods become competitive as the inertia term fully dominates the solution.

[chaosCFD]

Originally Posted by bigphil

Sorry, I have now fixed the URL to the paper; here it is again: P. Cardiff, Ž. Tuković, H. Jasak, A. Ivanković, A block-coupled Finite Volume methodology for linear elasticity and unstructured meshes, Computers & Structures, 175, 2016, https://doi.org/10.1016/j.compstruc.2016.07.004.


In reply to your questions

• Why is the segregated solver slow (relative to a coupled solver) for these cases? The segregated solver uses an approximation of the divergence of stress term in the momentum equation for the implicit term (a Laplacian term is used); then explicit (deferred corrections) are added so that once the outer iterations converge then the solution satisfies the original momentum equation (with the divergence of stress). This approach is efficient for "blocky" geometry with low aspect ratios; however, for bending of slender cantilevers, explicit terms are large relative to the implicit term resulting in slow convergence of the outer iterations. A coupled method treats all (or most) of the divergence of stress term implicitly, so it does not have this problem. Another solution could be to use a geometric multigrid algorithm, as done by Demirdzic, Muzaferija, Peric and co-workers.
  
• Why is the segregated solver much faster when solving the case transiently using small time steps?
  When using small time steps in transient problems, the inertia term (rho*dv/dt) is large relative to the stress term's divergence, so convergence is fast. In these cases, the coupled method could be slower because it is more expensive to build and solve the coupled linear system than the segregated system. For really small time steps, fully explicit methods become competitive as the inertia term fully dominates the solution.

I think I got it. Thanks so much!

[Steve1999]

Originally Posted by bigphil

As the time-step is reduced, all consistent methods will approach the same answer. The difference is that second-order methods will approach the time-step-independent answer at a quicker rate the first-order methods (well, actually only in the asymptotic region i.e. when getting close to the answer).

Also, for a given time-step size, second-order methods are often more accurate in absolute terms than first-order methods, but this does not have to be the case.

In your case, I suggest you perform the same time-step analysis using both Euler and backward and compare them, and you should see that backward approaches the final solution at a higher rate.

Philip

EDIT: by the way, first order Euler is generally known to require very small toe-steps to be accurate in comparison to Newmark or trapezoidal second order schemes, which can produce the same accuracy with much larger time-step sizes.

Hi Philip,

You mentioned that using the first-order Euler method or the second-order backward method can yield results similar to Newmark-beta, though they require different time-steps. However, I notice that in the field of structural dynamics (actually within FE), Newmark-beta and other similar methods, such as HHT-alpha, are commonly used. Are there significant differences between these two approaches (Euler and Newmark)? I rarely come across comparisons between them.

Additionally, I’m curious why Newmark-beta or HHT-alpha hasn't been implemented in solids4foam for solid dynamics problems. What are the reasons for this? I saw in one of your other posts that you mentioned, “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.” Have you managed to implement this in solids4foam, particularly in the context of the cell-centered finite volume framework? I'm trying to incorporate it and am considering an approach similar to backwardD2dt2Scheme.C in solids4foam. Do you think this is reasonable?

I look forward to hearing your thoughts.
Best,
Steve

[bigphil]

Originally Posted by Steve1999

Hi Philip,
You mentioned that using the first-order Euler method or the second-order backward method can yield results similar to Newmark-beta, though they require different time-steps. However, I notice that in the field of structural dynamics (actually within FE), Newmark-beta and other similar methods, such as HHT-alpha, are commonly used. Are there significant differences between these two approaches (Euler and Newmark)? I rarely come across comparisons between them.

The HHT-alpha scheme is a generalisation of the Newmark-beta scheme (i.e., you can get Newmark-beta behaviour with a particular choice of parameters), where the benefit of the HHT-alpha (as per my understanding) is that you have greater control of the frequencies which dissipate energy. For example, you can damp high-frequency numerical noise without significantly affecting physical lower frequencies. In contrast, Newmark-beta tends to dissipate energy across the spectrum of frequencies.

For a certain choice of parameters, HHT-alpha/Newmark-beta give the trapezoidal rule (same as Crank-Nicolson), which is energy conserving (consequently, it may be unstable).

HHT-alpha/Newmark-beta and backwards schemes are typically 2nd order, i.e., the numerical temporal error reduces at a second order rate as the time step size is reduced, e.g. halving the time step size will reduce the temporal error by a factor of 4, assuming you are in the asymptotic regime (close to the true solution).
In contrast, Euler is 1st order accurate, e.g. halving the time step size will reduce the temporal error by a factor of 2. In practice, the 1st order Euler method is highly dissipative for structural problems for reasonable time step sizes. Of course, the temporal error can be made arbitrarily small but this will likely require prohibitively small time steps for the 1st order Euler scheme.

I compared 1st Euler, with backwards and Newmark-beta (with parameters corresponding to the trapezoidal rule) in https://www.researchgate.net/publica...cs_in_OpenFOAM; see slide 34 (temporal discretisation) and slide 53 (comparison of predictions for a 3-D vibrating cantilever). In this case, it is based on a vertex-centred formulation (irrelevant, but FYI).

In terms of comparing Newmark-beta and backwards, they will be similar for dynamic problems, albeit backwards can no parameters to control the amount of dissipation.

Originally Posted by Steve1999

Additionally, I’m curious why Newmark-beta or HHT-alpha hasn't been implemented in solids4foam for solid dynamics problems. What are the reasons for this? I saw in one of your other posts that you mentioned, “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.” Have you managed to implement this in solids4foam, particularly in the context of the cell-centered finite volume framework? I'm trying to incorporate it and am considering an approach similar to backwardD2dt2Scheme.C in solids4foam. Do you think this is reasonable?

As you will see above, there are some implementations of Newmark-beta in solids4foam and related codes, but not yet for the segregated cell-centred solvers. There is no fundamental reason apart from finding time to do it. One inconvenience of the Newmark-beta scheme from an OpenFOAM perspective is that it is expressed in terms of the update of the displacement, velocity and acceleration (typically, in OpenFOAM, we only need the d2dt2 discretisation). So it means we would need to ensure the displacement, velocity and acceleration are all updated using the same consistent schemes. This is fine, but it likely means adding extra checks and code into the solver/solid-model.


Philip

[Steve1999]

Hi Philip,

Thank you for your reply. It seems that there isn’t much difference between the backward and Newmark methods, but I’m working on implementing it in a cell-centered solver and running comparisons.

In addition, as I try to set up my plane stress example, I reviewed `correct(vol sigma)` in `linearElastic.C` from s4f v1.1. I noticed that in this code, the calculation of `epsilon.zz()` at the n+1 iteration step uses two components of \sigma , which are clearly from the n iteration step. This approach seems different from that in FEM. Could this calculation method impact the results?

Furthermore, when I compared this with `linearElastic.C` in s4f v2.1, I noticed that plane stress is addressed by using different formulas for \lambda, but in `correct(vol sigma, epsilon)`, K_ * tr (\epsilon) is applied. As I understand it, if a user specifies E and \nu, then K_ should also vary for plane stress, just like \lambda, since K_ = 2/3* \mu + \lambda.

Do you think this understanding is correct? I look forward to your insights.
Best,
Steve