
[Sponsors] 
July 25, 2012, 05:15 
fvVectorMatrix & blockMatrixSolver

#1 
New Member
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5 
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 

July 27, 2012, 14:02 

#2 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 572
Rep Power: 19 
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 offdiagonal coefficients, and source() gives you the righthandside 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 

July 29, 2012, 03:57 

#3 
New Member
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5 
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 openfoam1.6ext, 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 

July 30, 2012, 16:49 

#4 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 572
Rep Power: 19 
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 blockcoupled 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 neighboursneighbours 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 straightforward 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 

July 31, 2012, 02:29 

#5 
New Member
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5 
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 Last edited by tiat; July 31, 2012 at 04:13. 

July 31, 2012, 04:10 

#6 
New Member
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5 
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 

July 31, 2012, 05:29 

#7  
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 572
Rep Power: 19 
Hi Tian,
Quote:
Quote:
Quote:
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 

July 31, 2012, 06:05 

#8 
New Member
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5 
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 componentwise 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 

July 31, 2012, 10:36 

#9  
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 572
Rep Power: 19 
Quote:
Quote:
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 

July 31, 2012, 10:55 

#10 
New Member
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5 
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 

April 24, 2015, 05:02 

#11 
Member
Patricio Bohorquez
Join Date: Mar 2009
Location: Jaén, Spain
Posts: 94
Rep Power: 8 
I wondered whether fvm::laplacianTranspose operator is available in foamextend3.1? Is there a foamer who has already implemented it? tetFem incorporates it. This operator will be very helpful to use pUCoupledFoam for nonNewtonian fluids ...


April 27, 2015, 22:27 

#12 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 572
Rep Power: 19 
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 

April 28, 2015, 09:12 

#13 
Member
Patricio Bohorquez
Join Date: Mar 2009
Location: Jaén, Spain
Posts: 94
Rep Power: 8 
Thanks for your reply Philip.
Very good presentation. The new functionalities that you are developing together with Hrv and coworkers will be really useful. I tried an implementation of these operators by migrating the code from an old 1.4.1dev version but it is not easy. BTW, could we evaluate the gradient at the boundaries with secondorder 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 secondorder accuracy of the interior cells. 

May 4, 2015, 16:17 

#14  
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 572
Rep Power: 19 
Quote:
For the segregated case I certainly have noticed that nonorthogonal correction on boundaries makes a huge difference to the accuracy at the boundaries for solid mechanics problems. By default OpenFOAM disables boundary nonorthogonal corrections, this results in the solution only being firstorder accurate for meshes with nonorthogonal faces at the boundary. Maybe this could be contributing to your problems in nonnewtonian 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 

Tags 
fvm::laplaciantranspose, pucoupledfoam 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
What does mean Operation Source in fvVectorMatrix ( .H())?  EmadTandis  OpenFOAM Programming & Development  2  April 29, 2011 06:03 
FvVectorMatrix questions  stephan  OpenFOAM Running, Solving & CFD  6  May 4, 2007 06:23 
Matrix  shuo  OpenFOAM Running, Solving & CFD  4  October 23, 2006 03:09 
FvVectorMatrix  cfdengineering  OpenFOAM Running, Solving & CFD  0  December 1, 2005 10:51 