Register Blogs Members List Search Today's Posts Mark Forums Read

 January 10, 2012, 11:31 Manipulating source fvc::grad(vsf1+vsf2) #1 Member   Andreas Ruopp Join Date: Aug 2009 Location: Stuttgart / Germany Posts: 31 Rep Power: 9 Hello, I need to manipulate the source term (fvc::grad(vsf1+vsf2)) of the momentum equation. The physical behaviour should be, that every identified internal face should be treated as a zero gradient condition: Code: ```int icounterinternalfaces = 0; forAll(phi,phii) { label P = mesh.faceOwner()[phii]; label N = mesh.faceNeighbour()[phii]; if(vsf1[P]<(vsf2[N]-vsf2[P])) { if(vsf1[N]

 January 26, 2012, 20:31 #2 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 9 Try something like this in your code. Code: ```volScalarField sumVf = vf1+vf2; surfaceScalarField sumVff = fvc::interpolate(sumVf); volVectorField gradVsf = fvc::grad(sumVff); const unallocLabelList& owner = mesh.owner(); const unallocLabelList& neighbour = mesh.neighbour(); const vectorField& Sf = mesh.Sf(); const scalarField& V = mesh.V(); scalarField& sumVfi = sumVf.internalField(); scalarField& issf = sumVff.internalField(); vectorField& igGrad = gradVsf.internalField(); forAll(owner, facei) { label P = owner[facei]; label N = neighbour[facei]; if (vf1[P] < (vf2[N]-vf2[P])) & (vf1[N] < vf3[N]) { igGrad[P] += Sf[facei]*(sumVfi[P] - issf[facei]) / V[P]; igGrad[N] += Sf[facei]*(sumVfi[N] - issf[facei]) / V[N]; } }``` This is based on the Gauss grad scheme. For every selected face it removes the original face contribution and adds the zero-gradient contribution. I'm not guaranteeing it will work but it's probably the easiest approach. mm.abdollahzadeh likes this.

 February 3, 2012, 12:46 #4 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 9 correctBoundaryConditions(vsf, gGrad) modifies the value of the gradient on all normal (non-coupled) boundaries so that the surface-normal component of gradient matches the surface-normal of the boundary condition. This correction has no influence on the internal field values of the gradient so, unless the boundary values are important to you, you can safely ignore this. Conveniently this function is static so you can call it yourself to correct the surface-normal component afterwards. Code: ```volScalarField sumVf = vf1+vf2; surfaceScalarField sumVff = fvc::interpolate(sumVf); volVectorField gradVsf = fvc::grad(sumVff); const unallocLabelList& owner = mesh.owner(); const unallocLabelList& neighbour = mesh.neighbour(); const vectorField& Sf = mesh.Sf(); const scalarField& V = mesh.V(); //- Apply correction to internalField scalarField& sumVfi = sumVf.internalField(); scalarField& issf = sumVff.internalField(); vectorField& igGrad = gradVsf.internalField(); forAll(owner, facei) { label P = owner[facei]; label N = neighbour[facei]; if (vf1[P] < (vf2[N]-vf2[P])) & (vf1[N] < vf3[N]) { igGrad[P] += Sf[facei]*(sumVfi[P] - issf[facei]) / V[P]; igGrad[N] += Sf[facei]*(sumVfi[N] - issf[facei]) / V[N]; } } //- Now apply correction to all coupled boundaries forAll(mesh.boundary(), patchi) { const unallocLabelList& pFaceCells = mesh.boundary()[patchi].faceCells(); const vectorField& pSf = mesh.Sf().boundaryField()[patchi]; const fvsPatchField& pssf = sumVff.boundaryField()[patchi]; if (sumVf.boundaryField()[patchi].coupled()) { const scalarField vf1P = vf1.boundaryField()[patchi].patchInternalField(); const scalarField vf1N = vf1.boundaryField()[patchi].patchNeighbourField(); const scalarField vf2P = vf2.boundaryField()[patchi].patchInternalField(); const scalarField vf2N = vf2.boundaryField()[patchi].patchNeighbourField(); // const scalarField vf3P = vf3.boundaryField()[patchi].patchInternalField(); const scalarField vf3N = vf3.boundaryField()[patchi].patchNeighbourField(); forAll(mesh.boundary()[patchi], facei) { if (vf1P[facei] < (vf2N[facei]-vf2P[facei])) & (vf1N[facei] < vf3N[facei]) { label celli = pFaceCells[facei]; igGrad[celli]] += pSf[facei]*(sumVfi[celli] - pssf[facei]) / V[celli]; } } } } fv::gaussGrad::correctBoundaryConditions(sumVf, gradVsf);``` Again, I haven't checked that this works but the concept is there.

 February 7, 2012, 17:00 summarising your last post #5 Member   Andreas Ruopp Join Date: Aug 2009 Location: Stuttgart / Germany Posts: 31 Rep Power: 9 Hello cliffoi, thank you for your very fast reply. Please correct me, if I summaries something wrong: 1) First we go through the internal field, calculating the gradients with: volVectorField gradVsf = fvc::grad(sumVff); Then we manipulate the internal field with the first loop. Ok. In file gaussGrad.C, after dividing the cell gradient with the volume, then comes: Code: `gGrad.correctBoundaryConditions();` 2) This, we have to manipulate again with your suggestion: Code: ```//- Now apply correction to all coupled boundaries forAll(mesh.boundary(), patchi) { const unallocLabelList& pFaceCells = mesh.boundary()[patchi].faceCells(); const vectorField& pSf = mesh.Sf().boundaryField()[patchi]; const fvsPatchField& pssf = sumVff.boundaryField()[patchi]; if (sumVf.boundaryField()[patchi].coupled()) { const scalarField vf1P = vf1.boundaryField()[patchi].patchInternalField(); const scalarField vf1N = vf1.boundaryField()[patchi].patchNeighbourField(); const scalarField vf2P = vf2.boundaryField()[patchi].patchInternalField(); const scalarField vf2N = vf2.boundaryField()[patchi].patchNeighbourField(); // const scalarField vf3P = vf3.boundaryField()[patchi].patchInternalField(); const scalarField vf3N = vf3.boundaryField()[patchi].patchNeighbourField(); forAll(mesh.boundary()[patchi], facei) { if (vf1P[facei] < (vf2N[facei]-vf2P[facei])) & (vf1N[facei] < vf3N[facei]) { label celli = pFaceCells[facei]; igGrad[celli]] += pSf[facei]*(sumVfi[celli] - pssf[facei]) / V[celli]; } } } }``` 3) Then we have to do finally: correctBoundaryConditions(vsf, gGrad); Which is similar to your call: Code: `fv::gaussGrad::correctBoundaryConditions(sumVf, gradVsf);` Is it correct? Thank you very much. You are a great help! Best regards Andy Last edited by andyru; February 7, 2012 at 17:24. Reason: Better explanation

 February 7, 2012, 17:30 #6 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 9 That's pretty much it. Let me know if it actually works...

 February 9, 2012, 15:14 I made progess #7 Member   Andreas Ruopp Join Date: Aug 2009 Location: Stuttgart / Germany Posts: 31 Rep Power: 9 Hello Cliffoi, the implementation of the source worked and I could manipulate the source term of the momentum equation. Here I did one correction to your suggestion: Code: ```... igGrad[P] += Sf[facei]*(sumVfi[P] - issf[facei]) / V[P]; igGrad[N] += Sf[facei]*(sumVfi[N] + issf[facei]) / V[N]; ...``` with using + instead of - for the neighbour cell gradient value, hence the surface normal vector in gaussGrad is pointing in opposite direction (from the neighbour cell point of view) and, I think, must use then opposite sign. Then I got almost the expected result in my testcase. However, I manipulated the source term of the momentum equation before the piso-loop. So, I have a flux-field, which is not divergence free, but approximately satisfies momentum (because this grad term has the same function as "fvc::grad(p)" in my incompressible flow). Now I have to do the same corrections in the piso-loop [having a correction (?) of the grad(p)/grad(vsf1+vsf2) values as I see in Ferziger & Peric] because I have fluctuations in the flux-values after the time-step ... So first, I must understand the piso-loop. I think, some new interesting questions will arise from my side soon. Sorry about that in advance... Best regards, Andy

 Tags fvc, gauss, grad, manipulate, source

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post GGerber OpenFOAM Installation 54 April 24, 2015 16:02 Greg Givogue OpenFOAM Programming & Development 3 March 4, 2011 18:18 allenzhao OpenFOAM Installation 127 January 30, 2009 20:08 hjasak OpenFOAM Post-Processing 69 April 24, 2008 01:24 jens_klostermann OpenFOAM Bugs 11 June 28, 2007 17:51

All times are GMT -4. The time now is 07:59.