[gionni]

Hi,
I'm developing a set of BC for a plasma simulation (2 fluids model with drift-diffusion). Implementing the set using Swak I get the expected solution (validated with COMSOL).
Implementing the same set using a codedMixed BC I get different values for density of particles and temperature (~30% difference). An exmple in the picture attached.

The BCs displaying this anomaly are electron charge density (rho_e) and energy density (n_eps): the other BCs are behaving in the right way. This is not completely surprising, since the electron dynamic is much faster than ions and neutral particles. This means though that the coded BC is correct, in general. As an additional information, in this kind of simulation the gradient of rho_e at the boundary is very high, I think this might be the source of the problem.
I tried different level of mesh refinement and timestep, the swak BCs always give the right results, always different from the codedMixed BC. The mesh is 2D, axisymmetric; the BC is applied to the external patches.

The swak version of rho_e BC:

    patch1
    {
        type            groovyBC;
        gradientExpression "0";
        valueExpression "0";
        fractionExpression "fraction";
        variables	(
			"M_ion=6.683e-26;"
			"q=1.60218e-19;"
			"alpha=De&normal()&normal();"
			"beta=(mue&E&normal())+sqrt(q*Te/M_ion);"
			"fraction=beta/(beta+alpha/mag(delta()));"
			);
        value           uniform 1.00011e-01;
    }

The codedMixed version of rho_e BC:

rhoeRR	{
code 
#{
// VARS INITIALIZATION
			const dictionary& physicalProperties = db().lookupObject<IOdictionary>("physicalProperties");
			scalar& q = dimensionedScalar(physicalProperties.lookup("q")).value();
			scalar& M_ion = dimensionedScalar(physicalProperties.lookup("M_ion")).value();
						
			// FIELDS LOOKUP
			const volTensorField& De 			= db().lookupObject<volTensorField>("De");
			const volScalarField& Te			= db().lookupObject<volScalarField>("Te");
			const volTensorField& mue 			= db().lookupObject<volTensorField>("mue");
			const volVectorField& E			= db().lookupObject<volVectorField>("E");
			const vectorField norman = patch().nf();
			int c; 
			scalar beta, alpha;
			scalarField fracField = scalarField(patch().size());
			const labelUList& celle = patch().faceCells();

			forAll(celle, f)
				{
					c = celle[f];
					
					alpha = De[c]&norman[f]&norman[f];
					beta = ((mue[c]&E[c])&norman[f])+sqrt(q*Te[c]/M_ion);
					fracField[f] = beta/(beta+alpha*patch().deltaCoeffs()[f]);
				}
				
						
			refGrad() = Zero;
			refValue() = Zero;
			valueFraction() = (fracField);
#};	 

}

    wall
    {
		type            codedMixed;
		name			rhoeRR	;
		refValue		uniform 1.00011e-01;
		refGradient	uniform 0;
		valueFraction 	uniform 0;
		
    }

I tried to define the drift (mue&E&normal()) and the flux (sqrt(q*Te[c]/M_ion)) in the solver to give both BCs the same input field, same result. I even rewrote normal and delta to get the same versors and face-center distance for both; same results.

I found an old thread where it was shown a situation in which swak approximates the BC imposing zeroGradient and changing the boundary cell; i tried to do this in my coded BC to see if that could be the case, and it's not, the results were horrible.

If I implement the BC using a Neumann formulation, taking the value of the field explicitly from the previous time step, swak gives the same result, while my codedBC seems to differ, which is another unexpected behaviour.

I understand that it's a very specific case, so any suggestion about ways in which swak behaves differently than a coded BC would be extremely helpful.

[JJJ94]

Since the dynamic of the electrons is so much faster, did you try to see some algebraic models able to increase the stability?

[gionni]

Thank you for the reply!
I'm actually working on that idea, but it remains the fact that the same BC implemented in the 2 ways should lead to the same result.
Let's say that I want to understand this difference from swak first, without changing the BC.
Maybe the way I'm implementing the coded is wrong, I want to be aware of this.