[Jose Manuel Olmos]

Hello everyone,

I'm trying to simulate the flow of a solution (containing a solute) inside a pipe. I'm using the icoFoam solver with some modifications on the code for including the transport of the solute (just after calculating the distribution of velocities):

[I]fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
}

// --- PISO loop
while (piso.correct())
{
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU);

// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();

if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}
}
#include "continuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
}
// --- Mass transport

dictionary soluteSolveDict = mesh.solutionDict().subDict("solvers").subDict("C_ solutes");
forAll(soluteNameList, i)
{
fvScalarMatrix MassTransport
(
fvm::ddt(soluteList[i]) + fvm::div(phi, soluteList[i], "div(phi,C_solutes)")
- fvm::laplacian(soluteList[i].D(), soluteList, "laplacian(D,C_solutes)")
);
MassTransport.solve(soluteSolveDict);
}

// --- Write run time

The solver is running and the distribution of concentrations is calculated. However, it only takes part from the inlet of the pipe.. I mean, the transport of the momentum and the concentration of the solute is not "synchronized. You can see this problem in the following images, that show the distribution of velocities and concentrations at the same time:

https://ibb.co/bQLCf2v
https://ibb.co/7XSj25K

Before running the simulation, I created a new field ("concentration") with the following boundary conditions:

dimensions [0 -3 0 0 1 0 0];

internalField uniform 0;

boundaryField
{
CAD_inlet
{
type fixedValue;
value uniform 0.01; //this is the concentration of the solute at
//the inlet
}

CAD_outlet
{
type zeroGradient;
}

CAD_wall
{
type zeroGradient;
}
}

I don't understand why the solution reach the outlet of the pipe (the velocity is non-null at this zone) but the solute is only at the first of the pipe. I would be very grateful if anyone could tell me if I'm taking any mistake on the modification of the solver or in the establishment of the boundary conditions.

Thanks in advance,
Jose Manuel

[GregNordin]

This code is similar to How to add transport of N solutes to icoFoam, but there is one difference in the Mass Transport section. The line

- fvm::laplacian(soluteList[i].D(), soluteList, "laplacian(D,C_solutes)")

in the code above is

- fvm::laplacian(soluteList[i].D(), soluteList[i], "laplacian(D,C_solutes)")

in the linked code, where the difference is in the 2nd argument to fvm::laplacian. I don't know if this accounts for the problem you are seeing, but I thought it might be worth pointing out.

On another note, can you share your full solver code? I can't get the code I linked to above to compile because it is for an old version of OF and I am using 2112 in a docker container, which doesn't have the "basicMultiComponentSolution.H" header file used with the old version of OF. I am trying to simulate microfluidic device geometries that generate various solute concentration gradients.

[Jose Manuel Olmos]

Hi Greg,


You're right, my code is based on the code found in How to add transport of N solutes to icoFoam. The difference in the arguments fvm::laplacian was only a mistake I had when I created the post.


I also asked the same issue in another section of the forum, you can find it here:


Transport of solutes


People said me that the code was ok, my problem was only the duration of my simulation. After that, the solute reaches the end of the pipe.


Here you have the full code solver:

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.

OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.

Application
diffusionIcoFoam

Description
Transient solver for incompressible, laminar flow of Newtonian fluids.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "pisoControl.H"
#include "solute.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"

pisoControl piso(mesh);

#include "createFields.H"
#include "initContinuityErrs.H"


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



Info<< "\nStarting time loop\n" << endl;

while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;

#include "CourantNo.H"

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
}

// --- PISO loop

while (piso.correct())


{
volScalarField rAU(1.0/UEqn.A());

volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);

adjustPhi(phiHbyA, U, p);

// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU);

// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())

{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();

if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();

}
}

#include "continuityErrs.H"

U = HbyA - rAU*fvc::grad(p);

U.correctBoundaryConditions();
}

// --- Mass transport

dictionary soluteSolveDict = mesh.solutionDict().subDict("solvers").subDict("C_ solutes");
forAll(soluteNameList, i)
{
fvScalarMatrix MassTransport
(
fvm::ddt(soluteList[i]) + fvm::div(phi, soluteList[i], "div(phi,C_solutes)")
- fvm::laplacian(soluteList[i].D(), soluteList[i], "laplacian(D,C_solutes)")
);
MassTransport.solve(soluteSolveDict);
}

// --- Write run time

runTime.write();


Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Info<< "End\n" << endl;

return 0;
}


// ************************************************** *********************** //



Also, you need the additional files ("createFields.H", "solute.H", "solute.C" and "speciesTable.H" ) that you can find in How to add transport of N solutes to icoFoam.

[GregNordin]

Thanks, Jose! After fixing one typo, everything runs and the results appear reasonable.

The typo is in this line:

dictionary soluteSolveDict = mesh.solutionDict().subDict("solvers").subDict("C_ solutes");

where at the end there is a space between "C_" and "solutes" that has to be removed.

One other item to note is that in the cavity case that is in the original How to add transport of N solutes to icoFoam, I changed the C_solutes solver in "system/fvSolution" to "PBiCGStab".

Thanks again for your help!

By the way, have you found this solver useful for your use cases? I plan to do some microfluidic device simulation where we deliberately set up a concentration gradient, and I want to evaluate different geometries that we can create with our high resolution 3D printer that we have built specifically for microfluidic device fabrication.