|
[Sponsors] |
September 27, 2019, 22:52 |
The velocity shows checkerboard feature
|
#1 |
Member
|
Hi Foamers,
I create a new steady state solver based on the standard buoyantSimpleFoam to study natural convection. Part of solid-fluid zones are modelled with porous media. Thus I add some features from the standard solver rhoPorousSimpleFoam. The major equation is coded as, 1. UEqn.H Code:
// Construct the Momentum equation MRF.correctBoundaryVelocity(U); tmp<fvVectorMatrix> tUEqn ( fvm::div(phi, U) + MRF.DDt(rho, U) + turbulence->divDevRhoReff(U) == fvOptions(rho, U) ); fvVectorMatrix& UEqn = tUEqn.ref(); UEqn.relax(); // Include the porous media resistance and solve the momentum equation // either implicit in the tensorial resistance or transport using by // including the spherical part of the resistance in the momentum diagonal tmp<volScalarField> trAU; tmp<volTensorField> trTU; if (pressureImplicitPorosity) { /**** tmp<volTensorField> tTU = tensor(I)*UEqn.A(); pZones.addResistance(UEqn, tTU.ref()); trTU = inv(tTU()); trTU.ref().rename("rAU"); fvOptions.constrain(UEqn); volVectorField gradp(fvc::grad(p)); for (int UCorr=0; UCorr<nUCorr; UCorr++) { U = trTU() & (UEqn.H() - gradp); } U.correctBoundaryConditions(); fvOptions.correct(U); ****/ } else { pZones.addResistance(UEqn); fvOptions.constrain(UEqn); // solve(UEqn == -fvc::grad(p)); // fvOptions.correct(U); if (simple.momentumPredictor()) // modify coeffs in rAU { solve ( UEqn == fvc::reconstruct ( ( // 2019-09-16 // - ghf*fvc::snGrad(rhok) - ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh) )*mesh.magSf() ) ); fvOptions.correct(U); } trAU = 1.0/UEqn.A(); trAU.ref().rename("rAU"); } Code:
{ // 1. p --> p_rgh // 2. comment implicit part // const volScalarField& psi = thermo.psi(); rho = thermo.rho(); rho.relax(); // Info<< "Pre rho max/min : " << max(rho).value() << " " << min(rho).value() // << endl; tmp<volVectorField> tHbyA; if (pressureImplicitPorosity) { // tHbyA = constrainHbyA(trTU()&UEqn.H(), U, p); } else { tHbyA = constrainHbyA(trAU()*UEqn.H(), U, p_rgh); } volVectorField& HbyA = tHbyA.ref(); tUEqn.clear(); surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*trAU())); // 2019-09-16 // surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rhok)*mesh.magSf()); surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); bool closedVolume = false; surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA)); MRF.makeRelative(fvc::interpolate(rho), phiHbyA); closedVolume = adjustPhi(phiHbyA, U, p_rgh); phiHbyA += phig; // Update the pressure BCs to ensure flux consistency constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF); while (simple.correctNonOrthogonal()) { tmp<fvScalarMatrix> tp_rghEqn; if (pressureImplicitPorosity) { // tpEqn = // ( // fvm::laplacian(rho*trTU(), p) // + fvOptions(psi, p, rho.name()) // == // fvc::div(phiHbyA) // ); } else { tp_rghEqn = ( fvm::laplacian(rhorAUf, p_rgh) // + fvOptions(psi, p, rho.name()) == fvc::div(phiHbyA) ); } fvScalarMatrix& p_rghEqn = tp_rghEqn.ref(); p_rghEqn.setReference ( pRefCell, getRefCellValue(p_rgh, pRefCell) ); p_rghEqn.solve(); if (simple.finalNonOrthogonalIter()) { // Calculate the conservative fluxes phi = phiHbyA - p_rghEqn.flux(); // Explicitly relax pressure for momentum corrector p_rgh.relax(); // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure U = HbyA + trAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf); U.correctBoundaryConditions(); fvOptions.correct(U); } } #include "continuityErrs.H" // 2019-09-16 // p = p_rgh + rhok*gh; p = p_rgh + rho*gh; /**** if (pressureImplicitPorosity) { U = HbyA - (trTU() & fvc::grad(p)); } else { U = HbyA - trAU()*fvc::grad(p); } U.correctBoundaryConditions(); fvOptions.correct(U); ****/ // pressureControl.limit(p); // For closed-volume cases adjust the pressure and density levels // to obey overall mass continuity /**** if (!thermo.incompressible() && closedVolume) { p += (initialMass - fvc::domainIntegrate(psi*p)) /fvc::domainIntegrate(psi); p_rgh = p - rho*gh; } ****/ rho = thermo.rho(); // depends on p and T rho.relax(); Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; // update display var // gradP_rgh = fvc::grad(p_rgh); // gradRho = fvc::grad(rho); // driveSource = -gradP_rgh - gh_display*gradRho; } Code:
{ volScalarField& he = thermo.he(); /**** dimensionedScalar Cp1 ( "liquidHeatCapacity", dimEnergy/dimMass/dimTemperature, thermo.Cp() ); volScalarField Cp = thermo.Cp(); ****/ gammaPhi = gammaf*phi; volScalarField alphaMixEff ( "alphaMixEff", gamma*turbulence->alphaEff() + (1.0-gamma)*lambda2/Cp1 ); // h1 -> fluid // h2 -> solid // h2 = Cp2 *T = Cp2 * h1/Cp1 // laplacian(alpha1*gamma, h1) + laplacian(alpha2*(1-gamma), h2) //=laplacian(alpha1*gamma, h1) + laplacian(alpha2*Cp2/Cp1*(1-gamma), h1) // alpha = conductivity/heat capacity [kg/m/s] // Cp [J/kg/K] // lambda: conductivity [kg m /s^3 /K] // alpha = nu / prandtl // the following he refers to enthalpy of fluid fvScalarMatrix EEqn ( fvm::div(gammaPhi, he) + ( he.name() == "e" ? fvc::div(gammaPhi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho)) : fvc::div(gammaPhi, volScalarField("K", 0.5*magSqr(U))) ) - fvm::laplacian(alphaMixEff, he) == rho*(U&g)*gamma + radiation->Sh(thermo, he) + fvOptions(rho, he) ); fvScalarMatrix EEqnFvO ( - fvOptions(rho, he) ); const scalarField& V = mesh.V(); const scalarField& eqnSource = EEqnFvO.source(); scalarField& heSource = volHeatSource; heSource = eqnSource/V; // volHeatSource = EEqnFvO.source()/mesh.V(); // forAll(mesh.C(), celli) // { // volHeatSource[celli] = eqnSource[celli]/V[celli]; // } EEqn.relax(); fvOptions.constrain(EEqn); EEqn.solve(); fvOptions.correct(he); thermo.correct(); radiation->correct(); Cp1 = thermo.Cp(); // 2019-09-16 // rhok = rho*(1.0 - beta*(T - TRef)); // rhok.relax(); } However I found the velocity along the channel axis shows checkerboard feature, which is located in porous zones. However, there is no checkerboard out of the porous zones. And the temperature is good everywhere. In details, there are a lot of channel (in which the cross section is square) hang in the space. The fluid will flow into the channel from bottom and flow out the channel from top, as shown in the figure. There are spaces between two channels, and we call it as interchannel. The fluid is flow from top to bottom in interchannels, which is also shown in the figure. Besides, there are holes in the channel wall so that the fluid may flow through the holes aside. The channel wall is modelled as baffles. The major difference between channel and interchannel lies in: 1. the fuild in channels is heated by volumetric source and interchannel no vol heat sources; 2. there are solid structure in channels so the channel is modelled as porous media, whereas interchannel is pure fluid. Uz: p_rgh: T: I also investigated the reason of checkerboard feature in the forum. The problem should be solved by Rhie-Chow interpolation, make a "pseudo-staggered" version of the pressure equation in co-located grid, as discussed in, Differences in solution method for pisoFoam and buoyantBoussinesqPisoFoam Causes of checkerboarding I swithed the divSchemes from upwind to bounded Gauss limitedLinear but it doesnt help. I dont understand why velocity field has checkerboard feature only in porous zones. The source of porous resistance is add as Code:
pZones.addResistance(UEqn); Code:
forAll(cells, i) { const label celli = cells[i]; const label j = this->fieldIndex(i); const tensor D = dZones[j]; const tensor F = fZones[j]; AU[celli] += mu[celli]*D + (rho[celli]*mag(U[celli]))*F; }
__________________
Kai |
|
Tags |
baffles, checkerboarding, natural convection, porous domain |
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Problem with Velocity Poisson Equation and Vector Potential Poisson Equation | mykkujinu2201 | Main CFD Forum | 1 | August 12, 2017 14:15 |
Superimposing 2 Velocity fields For Model inputs in Comsol | f.black | COMSOL | 0 | August 2, 2017 10:35 |
atmBoundaryLayerInletVelocity - Velocity Profile not continuous through domain | sdfij6354 | OpenFOAM Running, Solving & CFD | 3 | July 26, 2017 17:16 |
Neumann pressure BC and velocity field | Antech | Main CFD Forum | 0 | April 25, 2006 03:15 |
what the result is negatif pressure at inlet | chong chee nan | FLUENT | 0 | December 29, 2001 06:13 |