Including Enthalpy Transport
This is looking a bit hard to comprehend without proper symbols but I've done my best to make it clear...
I'm trying to include an enthalpy transport term in my solver with binary diffusion, that term being:
divergence(sum_over_all_species(enthalpy_species_i * diffusion_mass_flux_species_i))
diffusion mass flux for species i is defined as :
rho * binary_diffusion_coefficient * gradient(species_concentration_i)
what I've done so far to try to implement this is store the enthalpy of each chemical species into fields then sum the laplacian that arrises from combining the top two things:
sum_over_all_species( laplacian ( rho*binary_diffusion_coefficient*enthalpy_specise_ i, species_concentration_i))
this compiles fine but when I try to solve I get an error that OF cannot sum over two different species concentration fields:
--> FOAM FATAL ERROR:
incompatible fields for operation
[N2] += [Kr]
From function checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)
in file /opt/openfoam211/src/finiteVolume/lnInclude/fvMatrix.C at line 1303.
Is there anyway around this? I tried using mag(Yi) in my laplacian but that throws a different error. Some way to remove the species tag when I take the laplacian?
Can you post mathematical formulation in a) differential form, b) integral form and c) numerical integral form?
The term I'm trying to implement is somewhat from the fluent user manual:
neglecting the thermal diffusion effects. The enthalpy transport term I'm trying to include in the energy equation is:
where Ji is the mass flux of species i given by:
hi is the enthalpy of species i
Di,m is the diffusion coefficient of species i in the mixture
rho is the density
Yi is the mass fraction of species i
I should note that the base energy equation I'm adding this to is essentially the rhoCentralFoam energy equation (re-written for sensible enthalpy of course) and reactingFoam's species equation with the relevant diffusion coefficient modeling implemented.
I can solve for each value of the summation fine, as mentioned above, but when I try to sum the resulting laplacians (if you substitute Ji into the top equation) I get the error mentioned.
I was able to get it to solve by commenting out the relevant check method... (not exactly the proper way to do it) but I'm getting some interesting behavior at the boundaries. At any boundary with a uniform species concentration equal to 1 (not sure if its specific to 1, probably not) the enthalpy field drops to zero. If it's in a small region, the simulation is able to recover and the enthalpy field returns to normal. At least I think it does, I don't have a simulation run to steady state yet but when I last looked at it the region in question was returning to normal.
A zeroGradient boundary condition does not exhibit this behavior, it seems to be unique to the uniform boundary condition. This seems odd because the gradient of the species concentration should be/is zero for a uniform boundary condition and so the enthalpy transport term should also drop to zero.
It would seem that the issue is that the enthalpy transport term for the boundary during the first (few?) time steps is roughly equal to the negative total enthalpy on the boundary for a uniform boundary condition. I think it might be due to the 2nd derivative of species concentration not being available for the first timestep? Does anyone else have any theories to why I'm seeing this behavior and/or a way around it in case I end up with a simulation where it can't recover.
My implementation is the following:
outside of my species equation loop I find the inertIndex (index of the inert species) then initialize a temporary fvscalarmatrix:
then inside my species equation loop, under the if(Y[i].label != inertSpecies) I sum the other species terms
hsJ() += fvc::laplacian(-rho*Dij*hs*Y[i],Y[i])
The hsJ term is then added to the enthalpy corrector equation in the main code.
I don't have access to my code at the time of writing so this is from memory, I'll post my Yeqn.h at my next opportunity.
You have error in dimensions, but i can't find it still, i need some time to switch to your case
Thanks for taking a look at the issue.
If it helps, the only dimensions I specified are for the diffusion coefficient when I initialized the field, actually its rho*Dij initialized as rhoDij, with the following dimensions:
dimensionedScalar("rhoDij", dimensionSet(1,-1,-1,0,0), 0.0)
which should be consistent with rho (kg/m^3) multiplied by D (m^2/s). The rest are defaults from the original code.
Although, given the error, and the section of the checkMethod function it seems like physical dimensions are not the issue, it's looking at the psi fields:
if (&fvm1.psi() != &fvm2.psi())
what does fvm.psi() call?
I also think I've figured out the enthalpy dropping to zero issue, I found the solution here:
looks like I wasn't giving the diffusion coefficient a face value for the boundary condition, leading to incorrect (large) enthalpy transport at the boundary.
I think I shouldn't have been using fvm::laplacian and an fvmatrix for my enthalpy transport term as hs is not the laplacian variable, therefore the enthalpy transport term boils down to a source/sink term for each cell. I have it as a volScalarField with fvc::laplacian solving the terms. Can anyone confirm if this is correct? This is getting into the realm of calculus where my memory is a bit fuzzy.
|All times are GMT -4. The time now is 11:32.|