Good evening to all, I'm a phD student, and I'm trying to make some diffusion calculations of hydrogen in metals (applied to permeation tests, that is, basically, taking a metal plate, exposing one face to hydrogen, and recording the hydrogen getting out of the sample at the plate exit on the other face).
First of all I'm trying to simulate pure fickian diffusion through the plate; but though I've read and read again the documentation, I still haven't understood how the "real time during which the equation must be solved" has to be given to the solver, or put in the program. The program I've written (end of this e-mail), is similar to the first mesh1D example of the documentation, but when comparing the hydrogen profile to the erfc solution, great discrepancies arise. It's as though the solver time went faster than the analytical solution time, since the solver has already attained the steady state when the hydrogen penetration is just at its beginning according the analytical solution. So, my question is: how should the "real time" duration be given to the solver and/or the program? - is "duration" equal to the number of iterations calling the 'eq.solve(...)' times the 'dt' given to the solver, which is what I understood from the documentation (as in the example mesh1D: t=Steps*TimeStepDuration, at least for explicit method); and in that case, where is my mistake ? - or is it something else, and in that case, could someone please put me on the right track? Thanking you in advance for your kind attention, With my best regard, Fanny Jambon mailto:[email protected] ------------ >> Program ##-----------------------------------------------------Mesh X_scale= n.linspace(0,2.0e+4,1000)# grid 0-> 20 um ; units [nm] Y_scale= n.zeros((shape(X_scale))) # zero concentration at t=0 Nbx=len(X_scale) deltax= 1.0/float(Nbx) mesh=Grid1D(nx=Nbx, dx = deltax) ##-----------------------------------------------------Data D = 2.0e+4 # diffusion coefficient in [nm2/s] = D(cm2/s)*1.E+14 Cs = 1.0 # H concentration at entry surface deltat= 1 # time step [s] duration=15*60 # H charging test duration [s] ##-----------------------------------------------------Problem definition cH=CellVariable(name="H concentration", mesh=mesh, value= Y_scale) eq= TransientTerm() == DiffusionTerm(coeff = D) # Fick BCs = (FixedValue(faces = mesh.getFacesRight(), value=0), FixedValue(faces=mesh.getFacesLeft(),value=Cs)) # entry side: left ##-----------------------------------------------------...and solving elapsed =0 while elapsed <= duration: eq.solve(var = cH, boundaryConditions = BCs, dt = deltat) elapsed=elapsed+deltat ##--------------------------------------Comparison with the erfc solution Sol=1-erf(n.array(X_scale)/(2*sqrt(D*duration))) ##---------------------------------------------------Graph plotting (pylab) Y_scale=tuple(cH) plt.plot(X_scale,Y_scale, marker='x') plt.plot(X_scale, Sol, marker='o') plt.show()
