[tiat]

hi, dear all

I'm a student doing my project with stress analysis using OpenFOAM, and I would like to ask about a question related to the fvVectorMatrix UEqn in the SolidDisplacementFoam.

Is it possible to obtain a fvScalarMatrix for a component of U from UEqn, like UEqn[0], UEqn[1], UEqn[2] (i tried, but failed..) but i guess there can be some way to access the information for the single component of U?

The reason why i would like to obtain the fvScalarMatrix for the Ux, Uy, Uz, is because i am trying to use the block matrix solver instead of the segregated method, so all the components could be solved implicitly...

i appreciate any of your kind suggestions

Tian

[bigphil]

Hi Tian,

It is good to hear that other people are looking at the block solver for stress analysis. I have started to look at it myself and derive the coupled coefficients, it starts to get very long. I plan to get back looking at this in the coming months.

To answer your question, you can obtain all the fvVectorMatrix coefficients using the fvMatrix and lduMatrix class functions e.g. diag() are the diagonal coefficients, upper() and lower() are the off-diagonal coefficients, and source() gives you the right-hand-side source vector.

Out of interest, which terms do you plan on making implicit and how do you plan to implicitly discretise them?

Best regards,
Philip

[tiat]

Hi, Philip

thanks a lot for the reply! i'm so happy that you have the same interest as me about using block matrix for the stress analysis.

i have been reading through the BlockScalarTransportFoam in openfoam-1.6-ext, and i plan to use that idea for my case.

For the momentum equation, if i understood it correctly, the traditional way of decomposition by Jasak can give the coefficients in the diagonal terms of the block.

fvm::laplacian(2*mu+lambda, U)==fvc::div(sigmaExp), sigmaExp=mu*fvc::grad(U).T()+lamda*tr(fvc::grad(U) )-(mu+lambda)*fvc::grad(U).

therefore i was planning to construct a BlockLduMatrix as the following:

1. fvVectorMatrix UEqn(fvm::laplacian(2*mu+lambda,U)) to update the diagonal terms in the block (both in the diagonal blocks and the upper lower blocks)
2. then i need some other fvMatrices for updating the upper and lower terms (six in total) in the block (again in the digonal blocks and upper lower blocks). The problem comes out how to derive the coeffcients from ∂^2(u_y)/(∂x∂y), ∂^2(u_z)/(∂x∂z)in the ux equation, and similarly ∂^2(u_x)/(∂y∂x), ∂^2(u_z)/(∂y∂z) in uy equation, and ∂^2(u_x)/(∂z∂x), ∂^2(u_y)/(∂z∂y)in uz equation. Since we only have the fvm::laplacian() discretization in Openfoam, no fvm::grad() available yet, then it becomes quite difficult to get the coefficients i mentioned above...(sorry for the bad format of partial derivatives, i hope you can understand my messy expressions )

how do you think? it would be great to listen to your ideas

Tian

[bigphil]

Hi Tian,

Yep I believe you are on the right track.

The steady state momentum equation for solids neglecting body forces is:
div(sigma) = 0

To solve this using the OpenFOAM segregated approach, we write it as:
div(sigma_implicit) + div(sigma_explicit) = 0,
where
div(sigma_implicit) = fvm::laplacian(2mu+lambda, U)
and
div(sigma_explicit) = div( mu*gradU.T() + lambda*tr(gradU)*I - (mu + lambda)*gradU )

Now, the block-coupled solver potentially allows the div(sigma_explicit) terms to be included implicitly.

The first explicit term can carry a large contribution (especially during simulations with significant bending), so we need a "fvm::laplacianTranspose" operator to discretise it:
div(mu*gradU.T()) = fvm::laplacianTranspose(mu, gradU)

Unfortunately, this operator has not been written yet, and this is something I have started to look at. Essentially what I have started to do was discretise the div(mu*gradU.T()) term and then end up with the coefficients for the block matrix. The derivation becomes quite lengthy, and neighbours-neighbours term must be treated explicitly (because OpenFOAM only allows neighbours be implicit). Also care must be taken with the boundary conditions, treating them explicitly is the most straight-forward solution.

At the moment I am trying to finish my thesis, but in a few months I will hopefully have a bit of time to see if I can finish what I started and get it working.

Best regards,
Philip

[tiat]

Hi Philip,

thanks a loooot for your kind reply
it's so good to know that you are planing to write the laplaciantranspose operator. it would be extremely helpful for the implicit way of solving stress problem. i'm looking forward it!
btw, does it take a long time to write such an implicit operator (like the laplaciantranspose you are planning)? i'm also interested in a fvm::grad() operator since my project is tring to have the pressure coupling in the momentum equation as well. probably i could start to think about such an implicit gradient operator at least.

regards,
Tian

[tiat]

hi again Philip,

i would like to ask you one more question. if the laplacianTranspose operator can do the implicit discretization of div(mu*grad(U).T()), i guess there might be still a problem left for the block matrix, since we need the cofficients of u_y and u_z separately in u_x equation (similar for the u_y, u_z equation) and they are just mixed in the div(mu*grad(U).T()), have you thought about some way of decompose the div(mu*grad(U).T())?

i hope my question is not too silly and thanks for your time in advance

regards,
Tian

[bigphil]

Hi Tian,

btw, does it take a long time to write such an implicit operator (like the laplaciantranspose you are planning)?

Developing an implicit operator requires a good understanding of the finite volume method. Then operator must be tested and debugged which can be the slowest part. If you do plan to try derive such an implicit operator, I would recommend you become comfortable with the current fvm:: operators like fvm::laplacian. Hrv's slides describe it nicely (Finite Volume Discretisation with Polyhedral Cell Support ).

i'm also interested in a fvm::grad() operator since my project is tring to have the pressure coupling in the momentum equation as well. probably i could start to think about such an implicit gradient operator at least.

You could also check out developments by Mangani et al. at the recent OpenFOAM workshop. I think they may have already achieved an fvm::grad. Page 149 in the Book of Abstracts, 7th OpenFOAM Workshop. L. Mangani , M. Darwish2 and M. Buchmayr. A fully integrated pressure based coupled solver using block GAMG acceleration.

if the laplacianTranspose operator can do the implicit discretization of div(mu*grad(U).T()), i guess there might be still a problem left for the block matrix, since we need the cofficients of u_y and u_z separately in u_x equation (similar for the u_y, u_z equation) and they are just mixed in the div(mu*grad(U).T()), have you thought about some way of decompose the div(mu*grad(U).T())?

After successful derivation of the block coupled coefficient, you will be left with coefficient tensor a_p and coefficient tensor a_n such that:
a_p u_p + a_n u_n = b is the block coupled system.
So the derived coefficients can then just be inserted into the block system and solved. I am not quite sure what you mean when you say we need the coefficients separately… the coefficients tensors are composed of the individual term coefficients.

Best regards,
Philip

[tiat]

Hi, Philip,

sorry for my confusing question. let me try to explain it in a clear way:

To construct the a_p or a_n tensor cofficients (call it a) in the block matrix, if we use the origianl fvVectorMatrix UEqn(fvm::laplacian(2mu+lambda,U)), we obtain the a_11, a_22, and a_33, then if the laplacianTranspose works, we could have fvVectorMatrix UEqn2(fvm::laplacianTranspose(mu, U)), which possibly gives the other six elements, is my understanding correctly so far?

but one things confused me, laplacian(2mu+lambda, U) gives very clear component-wise structure (a_11*U_x, a_22*U_y, a_33*U_z), while the laplacianTranspose(mu, U) is somehow giving (a_11*U_x+a_12*U_y+a_13*U_z, a_21*U_x+a_22*U_y+a_23*U_z, a_31*U_x+a_32*U_y+a_33*U_z), (if my understanding is correct), but how could we seprarate out the a_12 and a_13, a_21 and a_23, a_31 and a_32 from the above structure... that's the really difficult part i could not understand. am I going a wrong way of learning the block matrix structure? could you give me some helps?

regards,
Tian

[bigphil]

Originally Posted by tiat

we could have fvVectorMatrix UEqn2(fvm::laplacianTranspose(mu, U)), which possibly gives the other six elements, is my understanding correctly so far?

fvm::laplacianTranspose would return a blockFvVectorMatrix (or blockLduMatrix) not a fvVectorMatrix, because fvVectorMatrix cannot hold the coupled coefficients.

Originally Posted by tiat

but one things confused me, laplacian(2mu+lambda, U) gives very clear component-wise structure (a_11*U_x, a_22*U_y, a_33*U_z), while the laplacianTranspose(mu, U) is somehow giving (a_11*U_x+a_12*U_y+a_13*U_z, a_21*U_x+a_22*U_y+a_23*U_z, a_31*U_x+a_32*U_y+a_33*U_z), (if my understanding is correct), but how could we seprarate out the a_12 and a_13, a_21 and a_23, a_31 and a_32 from the above structure... that's the really difficult part i could not understand. am I going a wrong way of learning the block matrix structure? could you give me some helps?

Essentially, you are correct. The vector momentum equation is made up of three scalar equation - the X, Y and Z momentum equations.
Each scalar momentum equation is function of (Ux, Uy and Uz), this is why the system is coupled.
For the segregated approach, in essence we assumes that the X momentum equation is only a function of Ux (we use old guesses of Uy and Uz). And then we do the same with the Y and Z equations.
But using the block solver, the X equation can have implicit Uy and Uz in it, and their coefficients give us the block coupled coefficients.

Philip

[tiat]

hi, Philip

now it becomes much clear for me. really thanks!

btw, could you please inform me any information when u have progessed with the laplacianTranspose operator in the future? i would really appreciate it! if this could be ok, here is my email, tiat@byg.dtu.dk.

best luck with your thesis

regards,
Tian

[pbohorquez]

I wondered whether fvm::laplacianTranspose operator is available in foam-extend-3.1? Is there a foamer who has already implemented it? tetFem incorporates it. This operator will be very helpful to use pUCoupledFoam for non-Newtonian fluids ...

[bigphil]

Hi Patricio,

I made quite a bit of progress on the implementation of fvm::laplacianTranspose and fvm::laplacianTrace (see OpenFOAM Workshop presentation here); but as yet it is not fully finished/stable/usable. I really hope to get back to this soon and finish it off; I will then be happy to share it.

Best regards,
Philip

[pbohorquez]

Thanks for your reply Philip.
Very good presentation. The new functionalities that you are developing together with Hrv and co-workers will be really useful. I tried an implementation of these operators by migrating the code from an old 1.4.1-dev version but it is not easy.

BTW, could we evaluate the gradient at the boundaries with second-order accuracy too? I found out a jump in the boundaries values of (for instance) the strain rate with respect to the interior cell. I believe that it is due to a first order evaluation of the gradient at the boundary faces instead of the second-order accuracy of the interior cells.

[bigphil]

Originally Posted by pbohorquez

BTW, could we evaluate the gradient at the boundaries with second-order accuracy too? I found out a jump in the boundaries values of (for instance) the strain rate with respect to the interior cell. I believe that it is due to a first order evaluation of the gradient at the boundary faces instead of the second-order accuracy of the interior cells.

Fully implicit block-coulpled strictly second-order boundary conditions are possible but may required a bit of effort.
For the segregated case I certainly have noticed that non-orthogonal correction on boundaries makes a huge difference to the accuracy at the boundaries for solid mechanics problems.
By default OpenFOAM disables boundary non-orthogonal corrections, this results in the solution only being first-order accurate for meshes with non-orthogonal faces at the boundary.
Maybe this could be contributing to your problems in non-newtonian fluids?
This can easily be fixed using custom boundary conditions, for example such as fixedDisplacement (in $FOAM_SRC/solidMechanics/fvPatchFields) instead of fixedValue and the equivalent for fixedGradient boundaries.


For the coupled approach, I am attempting to implicitly include these same corrections.

Philip

[Luca Cornolti]

Dear Philip,

I saw in your recent presentation "Updates on a Block Coupled Solver for
Linear Elasticity" that you implemented the three discretization operators of the laplacian term discussed in this topic.

Are you planning to release them in the next foam-extend version?

[bigphil]

Originally Posted by Luca Cornolti

Dear Philip,

I saw in your recent presentation "Updates on a Block Coupled Solver for
Linear Elasticity" that you implemented the three discretization operators of the laplacian term discussed in this topic.

Are you planning to release them in the next foam-extend version?

Hi Luca,

We recently had a paper accepted on this, see here:
http://www.sciencedirect.com/science...45794916306046

Yes, I plan to share the code; I am in the process of tidying things up.

Philip

Edit: the code from this paper is available in solids4foam, for example, see https://github.com/solids4foam/solid.../narrowTmember.

[fmuni]

Hello Philip,

I read your paper

We recently had a paper accepted on this, see here:
http://www.sciencedirect.com/science...45794916306046

and I found it very interesting.
Specifically, I am interested in the vertex-based ldu addressing and how it is implemented in foam-extend rather than the details of laplacianTranspose and laplacianTrace.
Such new addressing (and the implicit bonds) can enable the implementation of a variety of fvm discretisation operators in OpenFOAM (basically all operators that involve face tangent components)
and perhaps set the code free from the "abuse" of fixed-point iterations in linear systems.
Unfortunately, the paper does not explain in details how this is done.
I have been working with OpenFOAM for quite a while but i switched to foam-extend just recently, so I am not very familiar with BlockLduMatrices yet.
Perhaps it is a trivial matter.

Therefore, i would kindly ask if you are still planning to make the implementation public (ideally in the short period) so that you can share this formidable improvement with the community.
I suppose even a forgotten repo with almost broken code and hard coded stuff would be very helpful...surely better than nothing!

Thank you and best regards,
Federico

[randolph]

Originally Posted by tiat

hi, dear all

I'm a student doing my project with stress analysis using OpenFOAM, and I would like to ask about a question related to the fvVectorMatrix UEqn in the SolidDisplacementFoam.

Is it possible to obtain a fvScalarMatrix for a component of U from UEqn, like UEqn[0], UEqn[1], UEqn[2] (i tried, but failed..) but i guess there can be some way to access the information for the single component of U?

The reason why i would like to obtain the fvScalarMatrix for the Ux, Uy, Uz, is because i am trying to use the block matrix solver instead of the segregated method, so all the components could be solved implicitly...

i appreciate any of your kind suggestions

Tian

Tian,

I know this is quite old post, but did you manage to extract the fvScalarMatrix for each velocity component?

Thanks,
Rdf

[bigphil]

Originally Posted by fmuni

Hello Philip,

I read your paper

and I found it very interesting.
Specifically, I am interested in the vertex-based ldu addressing and how it is implemented in foam-extend rather than the details of laplacianTranspose and laplacianTrace.
Such new addressing (and the implicit bonds) can enable the implementation of a variety of fvm discretisation operators in OpenFOAM (basically all operators that involve face tangent components)
and perhaps set the code free from the "abuse" of fixed-point iterations in linear systems.
Unfortunately, the paper does not explain in details how this is done.
I have been working with OpenFOAM for quite a while but i switched to foam-extend just recently, so I am not very familiar with BlockLduMatrices yet.
Perhaps it is a trivial matter.

Therefore, i would kindly ask if you are still planning to make the implementation public (ideally in the short period) so that you can share this formidable improvement with the community.
I suppose even a forgotten repo with almost broken code and hard coded stuff would be very helpful...surely better than nothing!

Thank you and best regards,
Federico

Hi Federico,

Apologies for missing your post. In case you (or others) are still interested, the extended LDU addressing is available in solids4foam at https://github.com/solids4foam/solid.../solidPolyMesh. To be honest, it was a real pain implementing this; my current plans aim for extended addressing aim to avoid the OpenFOAM LDU addressing entirely, and instead write OpenFOAM discretisation operators that produce easier-to-use linear systems that can solved with, for example, PETSc. The vertexCentredLinGeomSolid solid model in solids4foam uses this approach and it seems to work fine.
This approach makes higher-order finite volume implementations seem very doable, especially since PETSc uses global indexing into the matrix.

Philip

[bigphil]

Originally Posted by randolph

Tian,

I know this is quite old post, but did you manage to extract the fvScalarMatrix for each velocity component?

Thanks,
Rdf

Hi Rdf,

It should be possible to do this by looking at the "solveSegregated" function in the fvMatrix class (see fvMatrixSolve.C).

The "solveSegregated" function forms the linear system for each direction before solving it.

fvMatrix stores the matrix without the boundary condition contributions because the boundary condition contribution may differ for each direction.
Before solving each direction, fvMatrix makes a copy of the matrix diagonal and modifies it with the boundary conditions for that direction before solving it; you could use similar code to prepare the matrix for each direction. For example, see how the addBoundaryDiag function is used inside the solveSegregated function:

   for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
    {
        if (validComponents[cmpt] == -1) continue;

	// copy field and source                                                                                                                                

        scalarField psiCmpt(psi.primitiveField().component(cmpt));
        addBoundaryDiag(diag(), cmpt);

        scalarField sourceCmpt(source.component(cmpt));

        // ...

Philip