mhamadmahdialloush / ufvm Goto Github PK

grad_f(:,1,iComponent) = (1-g_f).*grad(neighbours_f,1,iComponent) + g_f.*grad(owners_f,1,iComponent);
grad_f(:,2,iComponent) = (1-g_f).*grad(neighbours_f,2,iComponent) + g_f.*grad(owners_f,2,iComponent);
grad_f(:,3,iComponent) = (1-g_f).*grad(neighbours_f,3,iComponent) + g_f.*grad(owners_f,3,iComponent);

I think the code above is wrong， it should be changed as
grad_f(:,1,iComponent) = (1-g_f).*grad(owners_f,1,iComponent) + g_f.*grad(neighbours_f,1,iComponent);
grad_f(:,2,iComponent) = (1-g_f).*grad(owners_f,2,iComponent) + g_f.*grad(neighbours_f,2,iComponent);
grad_f(:,3,iComponent) = (1-g_f).*grad(owners_f,3,iComponent) + g_f.*grad(neighbours_f,3,iComponent);

Am I right?

1. In cfdCorrectNSSystemFields, we need ti run function 'setupPressureCorrection' before we correct pressure and velocity? What's the purpose of this function? Why don't add dphi to pressure directly?
   function setupPressureCorrection
% Get correction
dphi = cfdGetDPhi;
% Check if needs reference pressure
pRefValue = cfdReferencePressure;
if cfdNeedPressureLevel
pRefCell = cfdGetRefCell;
pRefValue = pRefValue + dphi(pRefCell);
end
% Update correction
dphi = dphi - pRefValue;
cfdSetDPhi(dphi);
end

2. When we correct mdot_f in cfdCorrectMdot.m, what's the purpose of coefficient 0.75? Is it a kind of under relaxation factor?
   mdot_f(iFaces) = mdot_f(iFaces) + 0.75*(FluxCf(iFaces).*pp(owners_f) + FluxFf(iFaces).*pp(neighbours_f));

In the book 'the FVM in CFD', the discretized form of source term of k eqution in formula(17.57) is:
FluxCc = rhoepsilon/kvolume,
FluxVc = Pkvolume.
And the discretized form of source term of epsilon eqution in formula(17.59) is:
FluxCc = C2rhoepsilon/kvolume,
FluxVc = C1epsilon/kPk*volume.
I don't understand how to get this expression by linearize the source term in formula(17.25) and formula(17.26). Can you show me the derivation? Thank you very much!

When I try the 'bump' case which is in the folder 'tutorials/compressible', I found that the function 'cfdGetKeyValueFromBlock.m' may return a value which is not the value of the key in the block.
For example, in the function 'cfdReadTimeDirectory.m', we will get the type of boundary by the code in line 197:
type = cfdGetKeyValueFromBlock('type', ['boundaryField/', theBCInfo.name], fieldFilePath);
and then we will get the value of boundary by the code in line 203:
value = cfdGetKeyValueFromBlock('value', ['boundaryField/', theBCInfo.name], fieldFilePath);
When a boundary's type is 'zeroGradient', such as the boudary 'bump' in file'T', there is no 'value'. But the code still return a value which is actually the value of 'inlet'.
The content in the file ‘P’ is as follows:

    bump
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    inlet
    {
        type            fixedValue;
        value           uniform 300;
    }

After I run the code in line 203, the results is as follows:

which doesn't match the file 'T'.

Hello, uFVM is a very good learning material. I have benefited a lot from its contents.

but I have a small problem in cfdAssembleStressTerm. m.

local_FluxVf = local_FluxVf + mu_f.*dot(tUGrad_f(:,:,iComponent)',Sf')'; % Add transpose term
What puzzles me in this code is why "+" is here？ I think it should be "-" here. Here are my reasons：

\int_V \nabla\cdot\left(\mu(\nabla\vec{u}+\nabla \vec{u}^\top)\right) dV
= \sum_f \mu_f \left[(\nabla\vec{u})_f\cdot\vec{E}_f +\nabla\vec{u}_f\cdot\vec{T}_f+ (\nabla\vec{u})_f^\top\cdot\vec{S}_f\right]

Because “\mu_f (\nabla\vec{u})_f\cdot\vec{E}_f” this term contributes to the matrix, and "\nabla\vec{u}_f\cdot\vec{T}_f+ (\nabla\vec{u})_f^\top\cdot\vec{S}_f " this term contributes to the source term.

and " \nabla\vec{u}_f\cdot\vec{T}_f+ (\nabla\vec{u})_f^\top\cdot\vec{S}_f" two formulas of this term are calculated by "+"

so I think the code should be
local_FluxVf = local_FluxVf - mu_f.*dot(tUGrad_f(:,:,iComponent)',Sf')'; % Add transpose term

This is the point that makes me confused. Did I miss anything?

I am looking forward to your reply and best wishes : )

Hello, when I look at the code in the cfdAssembleMassDivergenceTerm.m file, I found that the average and redistributed body forces(I guess the reason for ignoring this term is to ignore the body force) and transient flux(shown below) are not present in the code.

$$\frac{\bar{a^oD_f^v}}{V_f}(v_f^o-\bar{v_f^o})$$
compared to the formula (shown below).

$$\begin{aligned} v_f = & \bar{v_f} - \bar{D_f^v}(\nabla p_f - \bar{\nabla p_f}) + \bar{D_f^v}(\bar{B_f^v} - \bar{\bar{B_f^v}})\\ &+ \frac{\bar{a^oD_f^v}}{V_f}(v_f^o-\bar{v_f^o}) + (1-\lambda^v)(v_f^{(n)}-\bar{v_f^{(n)}}) \end{aligned}$$

Is this because uFVM is mainly used to solve steady-state problems or the omission of transient terms is also simplified ?

I am looking forward to your replay, best wishes : )

In function cfdRunFalseTransientCase, why cfdUpdateGradients is called only for one-time outside the main loop?

This problem results in that the gradient used throughout the whole computation process is the same! (Notice that in cfdAssembleAndCorrect_xxx_Equation functions, gradient is used to calculate non-orthogonal correction terms.) Thus, I think an additional call to cfdUpdateGradients should be placed at the end of the main loop body.

Or some changes needed?

Hello,

When going through the book and the code, the assembly of the stress term seems to have changed a lot from the equations in the textbook. Would it be possible to provide some references to understand the equations further?

The code snippets that I am confused is on
local_FluxCb = mu_b.magSb./wallDist_b.(1-nx2);
local_FluxVb = -mu_b.magSb./wallDist_b.(u_b.*(1-nx2)+(v_C-v_b).*ny.*nx+(w_C-w_b).*nz.*nx);

Initially it looks very similar to geoDiff however I am unsure about where the (1-nx2) term comes from

Thanks,
Marcus

Hi, I don't understand why the FluxVf in cfdAssembleMassDivergenceAdvectionTerm.m is -mdot_f.
local_FluxCf = drhodp_f./rho_f.*max( mdot_f,0);
local_FluxFf = -drhodp_f./rho_f.*max(-mdot_f,0);
local_FluxVf = -mdot_f;
We have already considered the mdot_f in cfdAssembleMassDivergenceTerm.m. I think FluxVf here should be zero. Is that right?

Look at the code for calculating geoDiff_f
In the 78 row,''n = cfdUnit(Sf);''
In the 79 row,''e = cfdUnit(CF);''
In the 81 row, ''CF_limited = max(dot(n',CF'), 0.05*cfdMag(CF)')';''
In the 84 row,''magCF = cfdMag(CF_limited);''
In the 105 row,''geoDiff_f = cfdMag(Ef)./magCF;''

I think the CF_limited is wrong, it should be ''CF_limited = max(dot(e',CF'), 0.05*cfdMag(CF)')';''

In the line

and

Should it change to
ac(nei) = ac(nei) - theFluxes.FluxCf(iFace);
anb{nei}(iNeighbourOwnerCoef) = anb{nei}(iNeighbourOwnerCoef) - theFluxes.FluxFf(iFace);

For

Cp_old = cfdGetPrevIterSubArrayForInterior('Cp');

Should we also use

Cp_old = cfdGetPrevTimeStepSubArrayForInterior('Cp'); like the density rho_old and temperature T_old?

Best.

Look at the code in the function cfdAssembleMassDivergenceTermInterior,
''
DUTf = DUSf - DUEf;

local_FluxCf = local_FluxCf + rho_f.*geoDiff;
local_FluxFf = local_FluxFf - rho_f.*geoDiff;
local_FluxVf = local_FluxVf - rho_f.*dot(p_grad_f(iFaces,:)',DUTf')'; %here maybe wrong

local_FluxVf = local_FluxVf + rho_f.*dot(p_grad_bar_f(iFaces,:)',DUSf')';
''
''local_FluxVf = local_FluxVf - rho_f.*dot(p_grad_f(iFaces,:)',DUTf')';'' this code maybe wrong, I think it is part of the formula (15.100) in book <>, so it should be changed to
''local_FluxVf = local_FluxVf - rho_f.*dot(p_grad_f(iFaces,:)',DUSf')';''

In the book, the FVM in CFD, it says there is a set of basic test cases that can be used to learn how to setup and run problems using the code.

testAdvection.m
testDiffusion.m
testDiffusion01.m
testDiffusion02.m
testDiffusion03.m
testDiffusion04.m
testFlow01.m
testFlow02.m
testGradients.m
testShearRate.m
testSmithHutton.m
testStepProfile.m
testStepProfile2.m
testTemperature.m
testTransient.m
testTurbulence.m

But I didn't find them in tutorals. I wonder where I can find these test cases, especially 'testTurbulence.m'.

When uFVM assembles the convection term of momentum equation, there is a function, cfdAssembleMomentumDivergenceCorrectionTerm.m, to claculate the divergence correction term. I can't understand the purpose of this function. Why we need to calculate the divergence correction term and add them to coefficient? Are cfdAssembleMomentumConvectionTerm.m and cfdAssembleMomentumDCSchemeTerm.m not enough for convection term?

    if strcmp(theTermName, 'Convection')
        cfdZeroFaceFLUXCoefficients;
        cfdAssembleMomentumConvectionTerm(iComponent);
        cfdAssembleMomentumDCSchemeTerm(iComponent);
        cfdAssembleIntoGlobalMatrixFaceFluxes;
        cfdZeroElementFLUXCoefficients;
        cfdAssembleMomentumDivergenceCorrectionTerm(iComponent);  % ？
        cfdAssembleIntoGlobalMatrixElementFluxes; 

I think this code is used for update velocity 'U' in boundaries.
However, for noslip boundary, the code is the same as zeroGradient boundary. I think it's wrong.

if strcmp(thePhysicalPatchType,'wall')
if strcmp(theBCType,'fixedValue')
updateFixedValue(iBPatch,theFieldName);
elseif strcmp(theBCType,'zeroGradient') || strcmp(theBCType,'noSlip') || strcmp(theBCType,'slip')
updateZeroGradient(iBPatch,theFieldName);
else
error([theBCType ' bc not defined']);
end

As Matlab has been able to run with Parallel Computing and compile to as C++ code ,I think it possible in code add LES.
Could you add an LES component?

Hi, I have a question about assembing stress term in uFVM.
When we assemble stress term at interior face,we need to add non-orthogonal part and velocity gradient transpose part to FluxVf, which is as followed:
% Non-orthogonal part
local_FluxVf = - mu_f.*dot(UGrad_f(:,:,iComponent)',Tf')';
% Add transpose term
local_FluxVf = local_FluxVf + mu_f.*dot(tUGrad_f(:,:,iComponent)',Sf')';

But we didn't do the same thing at boundary:
% non-linear fluxes
local_FluxVb = zeros(size(local_FluxFb));
Why?

Hello, I have a question in file cfdAssembleMassDivergenceAdvectionTerm.m.

The code(show blow) in the cfdAssemblyMassDivergenceTerm. m file, I think, is for solving this part $\sum_f\dot{m}^*_f$.

% Assemble Coefficients

% assemble term I
%     rho_f [v]_f.Sf

U_bar_f = dot(U_bar_f',Sf')';
local_FluxVf = local_FluxVf + rho_f.*U_bar_f;

% Assemble term II and linearize it
%      - rho_f ([DPVOL]_f.P_grad_f).Sf

DUSf = [DU1_f.*Sf(:,1),DU2_f.*Sf(:,2),DU3_f.*Sf(:,3)]; magDUSf = cfdMag(DUSf);
eDUSf = cfdUnit(DUSf);

DUEf = [magDUSf.*e(:,1)./dot(e',eDUSf')',magDUSf.*e(:,2)./dot(e',eDUSf')',magDUSf.*e(:,3)./dot(e',eDUSf')'];
geoDiff = cfdMag(DUEf)./cfdMagCF;

DUTf = DUSf - DUEf;

local_FluxCf = local_FluxCf + rho_f.*geoDiff;
local_FluxFf = local_FluxFf - rho_f.*geoDiff;
local_FluxVf = local_FluxVf - rho_f.*dot(p_grad_f(iFaces,:)',DUTf')';

%  assemble term III
%    rho_f ([P_grad]_f.([DPVOL]_f.Sf))

local_FluxVf = local_FluxVf + rho_f.*dot(p_grad_bar_f(iFaces,:)',DUSf')';

However, if the fluid being solved is a compressible fluid, the cfdAssemblyMassDivergenceAdvanceTerm function will be run after running the cfdAssemblyMassDivergenceAdvanceTerm function. However, I found the following code(The less important part in the middle has been omitted) in cfdAssemblyMassDivergenceAdvanceTerm. m which is also a part of $\sum_f \dot{m}^*_f$.

local_FluxVf = -mdot_f;                        
......

theFluxes.FluxVf(iFaces) = theFluxes.FluxVf(iFaces) + local_FluxVf;

This means that this term $\sum_f \dot{m}^*_f$ has been calculated(plus or minus) twice.

The above is my question, and I look forward to your answer and best wishes! : )

In the Region.mesh, there is a variable named wallDistLimited. What's the use of this variable? Is wallDist not enough for calculation?
wallDist(iBFace) = max(dot(faceCf(iBFace,:),n), cfdSMALL);
wallDistLimited(iBFace) = max(wallDist(iBFace), 0.05*cfdMag(faceCf(iBFace,:)));

I am a graduate student studying uFVM with the book FVM in CFD, and I find something puzzle.

In function cfdAssembleIntoGlobalMatrixFaceFluxes (in the file src/assemble/cfdAssembleIntoGlobalMatrixFaceFluxes.m),
lines include that

...
bc(own)                       = bc(own)                       - theFluxes.FluxTf(iFace);
...
bc(nei)                       = bc(nei)                       + theFluxes.FluxTf(iFace);
....
bc(own) = bc(own) - theFluxes.FluxTf(iBFace);
....

I suppose that bc(...) here represents the right hand side (RHS) of the algebraic equations (just like vector b of system Ax=b). But as the book writes in equation (8.80) in Chapter 8,

bc = Qc*Vc + non-orthogonal terms

where non-orthogonal terms should be FluxVf calculated in function cfdAssembleDiffusionTerm.

So I am wondering why this is not adding or subtracting the term theFluxes.FluxVf(...) but rather FluxTf(...) ? I notice that FluxTf(...) calculated in function cfdAssembleDiffusionTerm represents the total flux through the face, and could not understand why it should be included in bc(...).

for ZeroGradient boundary condition in cfdAssembleStressTermInletZeroGradientBC and cfdAssembleStressTermOutletZeroGradientBC,
the code is

'local_FluxCb = mu_b.*geoDiff_b;'
'theFluxes.FluxCf(iBFaces) = local_FluxCb;'

but according to the book FVM in CFD, the theFluxes.FluxCf(iBFaces) should be 0. I feel puzzle here

The imresize() function is part of the Image Processing Toolbox. Replacing this functionality with a custom function would increase accessibility of this code. I was able to comment out lines which use imresize() – for me this resulted in blank icons, but the GUI seemed to work fine otherwise.

Thanks for this contribution. Use of the Mathworks rms() function makes this code dependent on the Signal Processing Toolbox. Accessibility of this code could be increased by removing this dependency.

One option is to add a custom 'rms()' function to the source files. This two-line function worked for me:

function out = rms(in)
out = sqrt( mean(in(:).^2) );

Hello,

I wanted to clarify how the boundary face contributes to the gradient in the function cfdComputeGradientGaussLinear0.
In the for loop, phi_b(k) simply iterates across the first column of the matrix (x velocity), whereas phiGrad iterates in all three directions.

phiGrad(owners_b(k),:,iComponent) = phiGrad(owners_b(k),:,iComponent) + phi_b(k)*Sb(k,:);

I was wondering why this is the case or if I perhaps got the wrong idea.

Thanks,
Marcus

While assembing term III
The code is
"local_FluxVb = local_FluxVb + rho_b.*dot(p_grad_b',DUSb')'; "

But according formula(15.159) in Book FVM ib CFD, the term III should be
"local_FluxVb = local_FluxVb + rho_b.*dot(p_grad(owners_b)',DUSb')'; "

line 17 of cfdGetPrevTimeStepSubArrayForInterior checks for 'scfdUrfaceScalarField', I believe this should be 'surfaceScalarField'. Is this correct?

In "cfdComputeFormulaAtLocale" there is a call to "cfdGetFormulaForComponent" but I cannot find the last function, please check.

Hello, I have some questions about the discretization of the energy equation in uFVM, which have troubled me for some time.

In the book(The book written by your research group is excellent), the energy equation for an ideal Newtonian fluid is

$$ \begin{aligned} \frac{\partial}{\partial t}(\rho c_p T) + \nabla\cdot[\rho c_p\vec{v} T] &= \nabla\cdot[k\nabla T] + \rho T\frac{\mathrm{D}c_p}{\mathrm{D}t} + \frac{\mathrm{D}p}{\mathrm{D}t}+ \lambda\Psi + \mu \Phi+\dot{q}_V \\ p &= RT\rho \\ \Psi &= \nabla\vec{v}\bullet(\nabla\cdot\vec{v})\pmb{\mathrm{I}}=(\nabla\cdot\vec{v})^2, \quad \lambda = -\frac{2}{3}\mu \\ \Phi &= (\nabla\vec{v}+\nabla\vec{v}^\top)\bullet\nabla\vec{v} \end{aligned} $$

My question is about discrete implementations of the transient terms, convective terms, and diffusion terms of the energy equation that can be found in uFVM, but what about the other terms? Has the equation been simplified？

I look forward to your answer and best wishes : )