From b1a1597c3258d03043a9caf8d602da0950c9d55a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Chassagne?= Date: Fri, 12 Jan 2024 17:28:42 +0100 Subject: [PATCH 1/4] Fix Taua computation and allow granular phase momentum budget --- solver/createFields.H | 177 +++++++++++++++++++++++++++++++++++++++++- solver/writeOutput.H | 57 +++++++++++++- 2 files changed, 229 insertions(+), 5 deletions(-) diff --git a/solver/createFields.H b/solver/createFields.H index c33ebe9e5..44b3305e3 100644 --- a/solver/createFields.H +++ b/solver/createFields.H @@ -491,7 +491,10 @@ Switch writeTau ( forceProperties.getOrDefault("writeTau", false) ); - +Switch writeMomentumBudget +( + forceProperties.getOrDefault("writeMomentumBudget", false) +); Info<< "Initializing with specified pressure gradient:" << gradPOSC.value() << endl << endl; @@ -876,6 +879,178 @@ volScalarField p p_rbgh + rhob*gh ); + +volVectorField SolidPressureContr +( + IOobject + ( + "SolidPressureContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(p_rbgh) +); + +volVectorField divTauaContr +( + IOobject + ( + "divTauaContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + + +volVectorField ViscStressContr +( + IOobject + ( + "ViscStressContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + +volVectorField FrictContr +( + IOobject + ( + "FrictContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + +volVectorField phiRaContr +( + IOobject + ( + "phiRaContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + +volVectorField phiRa2Contr +( + IOobject + ( + "phiRa2Contr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + +volVectorField divRcaContr +( + IOobject + ( + "divRcaContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + +volVectorField RcaContr +( + IOobject + ( + "RcaContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + +volVectorField BulkContr +( + IOobject + ( + "BulkContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + +volVectorField gravityContr +( + IOobject + ( + "gravityContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + +volVectorField DragContr +( + IOobject + ( + "DragContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + +volVectorField SuspContr +( + IOobject + ( + "SuspContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + +volVectorField ExtContr +( + IOobject + ( + "ExtContr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::grad(pff) +); + + label pRefCell = 0; scalar pRefValue = 0.0; setRefCell diff --git a/solver/writeOutput.H b/solver/writeOutput.H index 731b9e714..ad8a33b4a 100644 --- a/solver/writeOutput.H +++ b/solver/writeOutput.H @@ -25,13 +25,62 @@ if (favreAveraging) } if (runTime.outputTime()) { - if (writeTau == 1) - { - Taua = rhoa*(alpha*nuEffa+nuFra)*dev(fvc::grad(Ua)); - Taub = rhob*(scalar(1.0)-alpha)*nuEffb*dev(fvc::grad(Ub)); + if (writeTau == 1){ + K.write(); + Taua = fvc::reconstruct(fvc::interpolate(rhoa*(alpha*nuEffa+nuFra))*mesh.magSf()*fvc::snGrad(Ua)) + + rhoa*(alpha*nuEffa+nuFra)*dev(gradUaT); + Taub = fvc::reconstruct(fvc::interpolate(rhob*(scalar(1.0)-alpha)*nuEffb)*mesh.magSf()*fvc::snGrad(Ub)) + + rhob*(scalar(1.0)-alpha)*nuEffb*dev(gradUbT); Taua.write(); Taub.write(); } + if (writeMomentumBudget == 1){ + + volTensorField Rca( + "Rca", + ((2.0/3.0)*I)*((nuEffa+nuFra/(alpha+alphaSmall))*tr(gradUaT) + + spherSigmaSGSa) + - (nuEffa+nuFra/(alpha+alphaSmall))*gradUaT + ); + + volTensorField Rca2( + "Rca2", + ((2.0/3.0)*I)*((alpha*nuEffa+nuFra)*tr(gradUaT) + + alpha*spherSigmaSGSa) + - (alpha*nuEffa+nuFra)*gradUaT + ); + volVectorField Ra( + "Ra", + -nuEffa*fvc::grad(alpha)/(alpha+alphaSmall) + ); + + surfaceScalarField phiRa( + -fvc::interpolate(nuEffa)*mesh.magSf()*(fvc::snGrad(alpha)) + /fvc::interpolate(alpha+ alphaSmall) + ); + + SolidPressureContr=-fvc::grad(pa+pff); + ViscStressContr=rhoa*alpha*fvc::laplacian(nuEffa, Ua); + FrictContr=rhoa*fvc::laplacian(nuFra, Ua); + phiRaContr=-rhoa*alpha*fvc::div(phiRa, Ua, "div(phiRa,Ua)"); + phiRa2Contr=rhoa*alpha*fvc::Sp(fvc::div(phiRa), Ua); + divRcaContr=-rhoa*alpha*fvc::div(Rca); + RcaContr=-rhoa*alpha/(alpha+alphaSmall)*(fvc::grad(alpha) & Rca); + BulkContr=fvc::grad(lambdaUa*tr(gradUaT)); + divTauaContr=ViscStressContr+FrictContr+phiRaContr+phiRa2Contr+divRcaContr+RcaContr+BulkContr; + + gravityContr=alpha/(alpha+alphaSmall)*(rhoa-rhob)*g; + DragContr=alpha*beta*K*(Ub-Ua); + ExtContr=rhoa*alpha*gradPOSC/(tilt*rhob +(1.0-tilt)*rhoa); + SuspContr=-SUS*K*beta*turbulenceb->nut() + *((SUS_I*iso-SUS_A*aniso) & fvc::grad(alpha)); + + divTauaContr.write(); + gravityContr.write(); + DragContr.write(); + ExtContr.write(); + SuspContr.write(); + } if (granularRheology.on()) { pa.write(); From d28dd402c6b9fb87b6aaeedf4bde0208beef153e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Chassagne?= Date: Fri, 12 Jan 2024 17:40:24 +0100 Subject: [PATCH 2/4] Small modif in writeTau before pull request --- solver/writeOutput.H | 1 - 1 file changed, 1 deletion(-) diff --git a/solver/writeOutput.H b/solver/writeOutput.H index ad8a33b4a..cfe207610 100644 --- a/solver/writeOutput.H +++ b/solver/writeOutput.H @@ -26,7 +26,6 @@ if (favreAveraging) if (runTime.outputTime()) { if (writeTau == 1){ - K.write(); Taua = fvc::reconstruct(fvc::interpolate(rhoa*(alpha*nuEffa+nuFra))*mesh.magSf()*fvc::snGrad(Ua)) + rhoa*(alpha*nuEffa+nuFra)*dev(gradUaT); Taub = fvc::reconstruct(fvc::interpolate(rhob*(scalar(1.0)-alpha)*nuEffb)*mesh.magSf()*fvc::snGrad(Ub)) + From 9ff5d524c9dd1d153bf1a03462a1df5b8a570cd4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Chassagne?= Date: Fri, 12 Jan 2024 17:52:00 +0100 Subject: [PATCH 3/4] Fix checkstyle --- solver/writeOutput.H | 39 +++++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/solver/writeOutput.H b/solver/writeOutput.H index cfe207610..eb0eb6056 100644 --- a/solver/writeOutput.H +++ b/solver/writeOutput.H @@ -25,15 +25,21 @@ if (favreAveraging) } if (runTime.outputTime()) { - if (writeTau == 1){ - Taua = fvc::reconstruct(fvc::interpolate(rhoa*(alpha*nuEffa+nuFra))*mesh.magSf()*fvc::snGrad(Ua)) + - rhoa*(alpha*nuEffa+nuFra)*dev(gradUaT); - Taub = fvc::reconstruct(fvc::interpolate(rhob*(scalar(1.0)-alpha)*nuEffb)*mesh.magSf()*fvc::snGrad(Ub)) + - rhob*(scalar(1.0)-alpha)*nuEffb*dev(gradUbT); + if (writeTau == 1) + { + Taua = fvc::reconstruct( + fvc::interpolate(rhoa*(alpha*nuEffa+nuFra)) * + mesh.magSf()*fvc::snGrad(Ua)) + + rhoa*(alpha*nuEffa+nuFra)*dev(gradUaT); + Taub = fvc::reconstruct( + fvc::interpolate(rhob*(scalar(1.0)-alpha)*nuEffb) * + mesh.magSf()*fvc::snGrad(Ub)) + + rhob*(scalar(1.0)-alpha)*nuEffb*dev(gradUbT); Taua.write(); Taub.write(); } - if (writeMomentumBudget == 1){ + if (writeMomentumBudget == 1) + { volTensorField Rca( "Rca", @@ -57,16 +63,17 @@ if (runTime.outputTime()) -fvc::interpolate(nuEffa)*mesh.magSf()*(fvc::snGrad(alpha)) /fvc::interpolate(alpha+ alphaSmall) ); - - SolidPressureContr=-fvc::grad(pa+pff); - ViscStressContr=rhoa*alpha*fvc::laplacian(nuEffa, Ua); - FrictContr=rhoa*fvc::laplacian(nuFra, Ua); - phiRaContr=-rhoa*alpha*fvc::div(phiRa, Ua, "div(phiRa,Ua)"); - phiRa2Contr=rhoa*alpha*fvc::Sp(fvc::div(phiRa), Ua); - divRcaContr=-rhoa*alpha*fvc::div(Rca); - RcaContr=-rhoa*alpha/(alpha+alphaSmall)*(fvc::grad(alpha) & Rca); - BulkContr=fvc::grad(lambdaUa*tr(gradUaT)); - divTauaContr=ViscStressContr+FrictContr+phiRaContr+phiRa2Contr+divRcaContr+RcaContr+BulkContr; + + SolidPressureContr = -fvc::grad(pa+pff); + ViscStressContr = rhoa*alpha*fvc::laplacian(nuEffa, Ua); + FrictContr = rhoa*fvc::laplacian(nuFra, Ua); + phiRaContr = -rhoa*alpha*fvc::div(phiRa, Ua, "div(phiRa,Ua)"); + phiRa2Contr = rhoa*alpha*fvc::Sp(fvc::div(phiRa), Ua); + divRcaContr = -rhoa*alpha*fvc::div(Rca); + RcaContr = -rhoa*alpha/(alpha+alphaSmall)*(fvc::grad(alpha) & Rca); + BulkContr = fvc::grad(lambdaUa*tr(gradUaT)); + divTauaContr = ViscStressContr + FrictContr + phiRaContr + + phiRa2Contr + divRcaContr + RcaContr + BulkContr; gravityContr=alpha/(alpha+alphaSmall)*(rhoa-rhob)*g; DragContr=alpha*beta*K*(Ub-Ua); From 290e79c646b6bca782e5a5199447eab07e2ca98e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Chassagne?= Date: Fri, 12 Jan 2024 17:53:51 +0100 Subject: [PATCH 4/4] Fix chackstyle --- solver/writeOutput.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/solver/writeOutput.H b/solver/writeOutput.H index eb0eb6056..039548b10 100644 --- a/solver/writeOutput.H +++ b/solver/writeOutput.H @@ -72,9 +72,9 @@ if (runTime.outputTime()) divRcaContr = -rhoa*alpha*fvc::div(Rca); RcaContr = -rhoa*alpha/(alpha+alphaSmall)*(fvc::grad(alpha) & Rca); BulkContr = fvc::grad(lambdaUa*tr(gradUaT)); - divTauaContr = ViscStressContr + FrictContr + phiRaContr + + divTauaContr = ViscStressContr + FrictContr + phiRaContr + phiRa2Contr + divRcaContr + RcaContr + BulkContr; - + gravityContr=alpha/(alpha+alphaSmall)*(rhoa-rhob)*g; DragContr=alpha*beta*K*(Ub-Ua); ExtContr=rhoa*alpha*gradPOSC/(tilt*rhob +(1.0-tilt)*rhoa);