Deal with Diffusivity as function of concentration and temperature in diffusion equation

General Discussion
1

Hi all,

in case hardening simulation the diffusion coefficient of Carbon in Austenite Dc should depend on local Carbon concentration and Temperature:

My first questionis: Is it really as easy as simply updating Dc in the time loop?
(Dc = DcC(cC,T):wink:
if you not could somebody give me suggestions on how to do it correctly.

Attached example code for 1D axisymmetric test case which gives a runtime error

call DcC at line 71
Exec error : Try to get unset x,y, …
– number :1
Exec error : Try to get unset x,y, …
– number :1

but works in 2D version.

Thank you for your support.

Gerald

load "msh3"
load "gsl"

// Declare  Variables

real t = 0;
real tend = 280*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*x)
	- int1d(Th, qfe=qf1pElump)((abs(x-Ra) < 1e-5) * beta*cP*v*x);
//	+ int0d(Th, 2)(beta*cC*v*x)
//	- int0d(Th, 2)(beta*cP*v*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;
//	plot(cC, value = true, fill = true);
}