Rotation-free Kirchhoff solver

[iBatistic]

@philipcardiff, I tried solving the Kirchhoff plate (rotation-free) equation for a supported plate and encountered some odd behavior. Instead of achieving fourth-order convergence, I am only getting second-order convergence (but with a lower error value). Interestingly, in a method of manufactured solutions (MMS) scenario involving a Hookean solid, I do achieve the expected fourth-order convergence. The key difference is that here we have two Laplace equations. I thought it would be worthwhile to share this with you to get any insights or suggestions you might have.

[philipcardiff]

Interesting.

Is this the L1 norm error? i.e. the average of the absolute error values.
If it were the max (L infinity) error, then maybe this is a boundary condition issue; if it is L1, then I am not sure.

I will look through the code and see if I notice anything

[philipcardiff]

I wonder if it is related to the source terms.

In laplacianFoam.py, you have

    laplacian = fvm.construct(T, 'Laplacian', gammaValue)
    bodyForce = fvm.construct(T, 'bodyForce')

I guess bodyForce is explicit? And I guess you use quadrature for the volume integral.

In kirchhoffPlate.py, you have an additional function, which you do not use in laplacianFoam.py:

sourceW = fvm.construct(M, 'Su')

I guess this is an implicit volume integral. How do you deal with this? You could assume the mid-point rule for the implicit component and then explicitly correct the calculation in the source with the quadrature approach. Or you could use the quadrature (extended stencil) approach implicitly.
Since this is the only function, which is different, I suspect this may be the cause of the 2nd order behaviour, i.e. if this is only second order then the entire solution will be second order.

Also, what does momentumMatrix.relax(relaxation=1) do?

[iBatistic]

I wonder if it is related to the source terms.

Yes, I believe it is!

In the code, bodyForce refers to pressure, which is constant across all cells, and I’ve added it to the right-hand side:

I’m calling it bodyForce because I’m reading the pressure from a file rather than defining it as a field.

sourceW is incorrect. I’m using M obtained from solving the first equation and multiplying it by the cell volume (which corresponds to the midpoint rule, as you mentioned). What I actually need is quadrature approach for M to get the correct value. Hmmm, I think that to get values of M at the cell Gauss points, we need to treat them using the same procedures as for the ones on the faces. I need to think about this.

Also, what does momentumMatrix.relax(relaxation=1) do?

As for momentumMatrix.relax, that’s implicit under-relaxation. I needed it for the pressure-displacement solver. In this case, it has no effect becouse under-relaxation is set to 1.

[iBatistic]

@philipcardiff, good news; I’m now getting fourth-order convergence! (all edges clamped plate)

In the plot, the empty circles represent the previous results (midpoint integration for the source term in the deflection equation). The filled circles correspond to the source term integration using quadrature points over the cell.

[philipcardiff]

Very cool!

How long did each mesh take for 2nd vs 4th order? Was the 4th order much slower?
It is amazing that the coarsest 4th-order error is orders of magnitude smaller than the finest 2nd-order mesh.

[iBatistic]

@philipcardiff, the finest mesh with linear interpolation achieves the same accuracy as quadratic interpolation on a coarser mesh. At this accuracy level, linear interpolation (on the finest mesh) takes 29.85 seconds, while quadratic interpolation (on the coarser mesh) requires 3.32 seconds. Cubic interpolation (fourth-order) completes in 4.35 seconds on the coarsest grid. These timings look promising.

EDIT I was using 7 Gauss points per face in all cases. When i put 1 Gauss point per face for linear interpolation (what we use in standard discretisation) I'm getting 22 s which is still notable difference.

[philipcardiff]

Very promising, especially since the 4th-order coarse solution is much more accurate than the 2nd-order finest solution.