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

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)
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.

Gerald

``````load "msh3"

// 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);
}``````