Can int0d deal with robin boundary condition?

General Discussion
1

Hi all,

I am trying to use FreeFem++ for calculation of carburizing operations (steel heattreament).

On a basic level these processes can be modelled as axisymmetric 1D.

The attached code works with Dirichlet boundary condition.
+ on(2, cC=cP);

But with int0d (code works with int1d in 2D axisymmetric version) I get the runtime error:

current line = 54
Assertion fail : (0)
line :4433, in file problem.cpp
Assertion fail : (0)
line :4433, in file problem.cpp
err code 6 , mpirank 0

Am I doing anything wrong when switching back to 1D or can int0d simply not yet deal with convection type BC?

I appreciate your support.

Kind regards
Gerald

load "msh3"

// Declare Variable

real t = 0;
real tend = 3*3600;
real dt = 1;

// Define the Domain
// Define the Mesh

real Ra = 4/2.0;

meshL Th = Sline(100, [Ra*x,0,0]);

// Define the Set of Finite Elements

fespace Vh(Th, P2);
Vh cC, v, cCOld;

// Define the Variational Problem + Initial / Boundary Condition / Material Properties

//IC

real cC0 = 0.18;
cC = cC0;

// BC

real beta = 1.25e-4;

real cP = 1.2;

//Material

real DC = 1.35e-5;

// Variational Problem

problem Carburize(cC, v) =
	  int1d(Th)(cC*v/dt*x + DC * dx(cC) * dx(v) * x)
	- int1d(Th)(cCOld*v/dt*x)
	+ int0d(Th, 2)(beta*cC*v*x)
	- int0d(Th, 2)(beta*cP*v*x);

// Solve the Problem

while(t <= tend){

	t += dt;

	cCOld = cC;
	Carburize;
}
2

Yes to day the code are missing, I will add in next version.

But you can do this:

problem Carburize(cC, v) =
  int1d(Th)(cC*v/dt*x + DC * dx(cC) * dx(v) * x)
- int1d(Th)(cCOld*v/dt*x)
//	+ int0d(Th, 2)(beta*cC*v*x)
//	- int0d(Th, 2)(beta*cP*v*x);
    + int1d(Th,qfe=qf1pElump) ( (abs(x-Ra)< 1e-5)* beta*cC*v*x) 
    - int1d(Th,qfe=qf1pElump) ( (abs(x-Ra)< 1e-5)* beta*cP*v*x) ;
// Solve the Problem

F. Hecht,

3

In the last commit, I have correct this problem…
in https://github.com/FreeFem/FreeFem-sources/commit/d4ea76cbda2f537d57781bd7c80f92595f128016

5

Thank you Prof. Hecht.

In the code posted previously I neglected to include the 2piRa factor (circumference) in the terms of the robin boundary condition (as the model is meant to represent 1D axisymmetry).
I think this will be necessary also when using in0d operator.

In the attached example the max. concentration reaches ~1% after 260 min of carburizing.

load "msh3"
load "gsl"

// Declare  Variables

real t = 0;
real tend = 260*60;
real dt = 1;

real[int,int] TDataPoints = [
	[0, 280*60, 310*60, 350*60],
	[930, 930, 850, 850]
];

gslspline fT(gslinterplinear, TDataPoints);
real T = fT(0);

// Define the Domain
// Define the Mesh

real ri = 0/2.0;
real Ra = 25/2.0;

meshL Th = Sline(250, [Ra*x,0,0]);

// Define the set of Finite Element

fespace Vh(Th, P1);
Vh cC, v, cCOld;

// Define the Variational Problem + Initial / Boundary Condition / Material Properties

//IC

real cC0 = 0.175;
cC=cC0;

// BC

real beta = 1.25e-4;

real[int,int] cPDataPoints = [
	[0, 40*60, 41*60, 260*60, 261*60, 350*60],
	[0.7, 0.7, 1.14, 1.14, 0.7, 0.7]
];

gslspline fcP(gslinterplinear, cPDataPoints);
real cP=fcP(0);

// Material

func real DcC(real c, real T){ return 47 * exp((-1.6*c)-(154811-27615*c)/(8.314*(T+273.15))); };
real Dc = DcC(cC0, T);

// Variational Problem

problem Carburize(cC, v, solver = CG) =
	  int1d(Th)(cC*v/dt*x + Dc*(dx(cC) * dx(v) * x + dy(cC) * dy(v) * x))
	- int1d(Th)(cCOld*v/dt*x)
	+ int1d(Th, qfe=qf1pElump)((abs(x-Ra) < 1e-5) * beta*cC*v*2*pi*Ra*x)
	- int1d(Th, qfe=qf1pElump)((abs(x-Ra) < 1e-5) * beta*cP*v*2*pi*Ra*x);
//	+ int0d(Th, 2)(beta*cC*v*2*pi*Ra*x)
//	- int0d(Th, 2)(beta*cP*v*2*pi*Ra*x);


// Solve the Problem

while( t <= tend ){

	cCOld = cC;

	T = fT(t);
//	Dc = DcC(cC,T);
	cP = fcP(t);

	Carburize;

	cout << "t " << t/60 << "T " << T << "cP " << cP << endl;

	t += dt;
6

In release 4.7 all is ok now.