| Author | Message | ||
     Britta Hagmeyer (britta) Member Username: britta Post Number: 6 Registered: 11-2006 |
Hello everyone! I've got some questions about the results being produced by FlexPDE. I've tried to validate the model I've created, and so far by checking the results of my model, a few questions have turned up to which I have no satisfying answers. I hope you can help me out! First of all, here is the code I've constructed. I've built a 3D Model of a diffusion system with the particle concentration as a variable and a particle source (constructed according to Fick's 1st law) at the bottom of the 3D cylinder. The source represents the flux of the particles after diffusing from a reservoir through a membrane. In all 2.5e-17 mMol of particles can diffuse into the cylinder. The particle concentration inside the cylinder changes according to Fick's 2nd law. The unit of the dimension of the system is µm. --------------------------------------------------- title 'Diffusion of dye in water at 25°C' coordinates cartesian3 variables u(threshold=0.1) definitions D = 5.34e2 { diffusivity µm**2/s } { dimensions of the system } r_pore = 21.5 !radius of the particle source meshheight = 500 r_mesh = 250 !thickness of membrane d_mem = 6 !flux of the source: V = (D*((2.5e-14-integral(u))/2.5e4-u)/d_mem)*1.178/100 { unit: [mmol/µm^2*s] } !timestep of the simulation timesteps = 1 initial values u = 0 equations div(D*grad(u)) = dt(u) { diffusion equation (Fick's 2nd law) } extrusion { 3D-structure } surface 'fluid_bottom' z = 0 layer 'fluid' surface 'cylinder_top' z = meshheight { Construction of the system } boundaries region "cylinder" start(0,r_mesh) arc(center=0,0) angle 360 limited region "bottom_pore" { restricting 'bottom_pore' to 'fluid_bottom'. } surface 'fluid_bottom' natural(u) = V { Boundary condition: Flux at 'fluid_bottom' in the region 'bottom_pore' = D*grad(u) } start(0,r_pore) arc(center=0,0) angle 360 time 0 to 20 by timesteps plots for t = 0 by timesteps to endtime contour(u) on x = 0 zoom (-100,0,200,200) range (0.01e-20,4e-20) surface(u) on z = 0 zoom (-100,-100,200,200) range (0.01e-20,4e-20) elevation(u*1e18) from (0,0,0) to (200,0,0) fixed range (0,35e-3) contour(u) on z=0 zoom (-100,-100,200,200) range (0.01e-20,4e-20) Summary export file='Summary.txt' merge noheader report(t) as 'time' report(V) as 'flux from source' report(u) as 'concentration' report(integral(u)) as 'particles in the cylinder' end --------------------------------------------------- There are some things about the results that I do not quite understand: Do the dimensions of the system have any influence on the calculation of the concentration? I've made some tests with systems of varying sizes and unmodified parameters and equations, and I observed that the concentration at the same coordinates in the respective systems remained constant, no matter the system dimensions. Does this mean that the particle concentration is being calculated with respect to a reference unit and doesn't take the dimensions of the system into account? Where does the reported concentration u (with REPORT(u) in the summary) come from? Is it the mean concentration of the system? How is it being calculated? What is its unit? How can I check how many mMol of particles are being released into the cylinder by the particle source? Thanks a lot! Britta | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 779 Registered: 06-2003 |
1) The meaning of the spatial coordinate numbers is totally determined by the coefficients you enter in your equations. In your case, the values of D and V determine what the coordinate values mean and how the material will diffuse. As long as the system boundaries are far enough away that the boundary conditions do not interact with the diffusion wave, the penetration pattern will be unaffected by system size. Integral(U) is simply the amount of U present in the system, and if it is concentrated near the input disk, the position of the boundaries will not matter. (This is merely physics; it is not special to FlexPDE.) 2) Since U is a field variable and Report(u) does not specify an evaluation position, the report is meaningless. The location of evaluation is unknown. REPORT(VAL(U,<x>,<y>,<z>)) will pick a specific place. REPORT(Integral(U)/Integral(1)) will give you the mean, but this is not very meaningful, as it will depend on the size of the volume as well as the amount diffused. 3) The units of U are whatever is implied by your choice of numeric values. FlexPDE does not know or care. 4) You can calculate the total input flux as Tintegral(Surf_Integral(-D*dz(U),"Fluid_Bottom","Bottom_Pore")). This should agree with Integral(U). [See "surf_integral" in the Help Index for integral forms.] | ||
| Britta Hagmeyer (britta) Member Username: britta Post Number: 7 Registered: 11-2006 |
> You can calculate the total input flux as Tintegral(Surf_Integral(-D*dz(U),"Fluid_Bottom","Bottom_Pore")). > This should agree with Integral(U). Dear Mr. Nelson, as you have posted before, I understand that "the area (or volume) integral of a divergence is equal to the surface integral of the normal component of flux across the bounding surface". I've set the threshold value of u from 0.1 to 0.003 and calculated several values. Whereas the values of integral(u) and Tintegral(Surf_Integral(V)) are almost identical, the values of Tintegral(Surf_Integral(-D*dz(U),"Fluid_Bottom","Bottom_Pore")) and integral(u) do not match. The value of the flux across the whole surface 'Fluid_Bottom', calculated with Tintegral(Surf_Integral(-D*dz(U),"Fluid_Bottom")) and integral(u) also don't match. What is the reason? What do I have to change in the system? Thanks again! Regards, Britta | ||
| Britta Hagmeyer (britta) Member Username: britta Post Number: 8 Registered: 11-2006 |
I've performed these calculations on the above code as well as on a system with a constant source, e.g. V = 6e-6. The differences between the integrals were not as big, but still noticeable. Britta | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 781 Registered: 06-2003 |
The mesh on the bottom pore is quite crude, so there is a smearing of the input flux (since the finite element interpolator cannot model the discontinuities). This means that an integral over only the pore leaves out some of the flux that is actually entering the system. The more appropriate integral is over the entire bottom surface (or enough of it to include the smearing of the flux). You can plot an elevation of the bottom flux to see this happening. Particles in cylinder=1.7346 Flux integral on pore=1.4967 Flux integral on bottom surface=1.7098 What you really care about, though, is the flux you asked for, which is TINTEGRAL(pi*r_pore^2*V). If you compare against this, you will see: Particles in cylinder=1.7346 Programmed flux integral=1.7470 This shows that although the approximated flux distribution is smeared, the total entering the system is held very close to the number you requested (this is because the finite element equations are based on conserving flux integrals). In any numerical approximation, there is a tradeoff between accuracy and expense. In FlexPDE we have tried to choose the default behavior to give a "reasonable" balance. If you want more accuracy, you can simply reduce ERRLIM. The accuracy and the cost will both go up. You can also manually force higher mesh density in the places you know will need it. MESH_SPACING=r_pore/10 on the boundary of the pore will reduce the smearing of the interpolated flux. The results are: Flux integral on pore=1.6669 Flux integral on bottom surface=1.7285 This has improved the spatial resolution, but has done nothing to increase the accuracy of the time integration. ERRLIM=1e-4 together with the dense mesh on the pore boundary increases the running time substantially. The results are: Flux integral on pore=1.7066 Flux integral on bottom surface=1.7395 Notice also that in this case a smaller error per step is balanced by a larger number of steps, so that the total error in the time integral is not improved: Particles in cylinder=1.7312 Programmed flux integral=1.7809 I should also point out that a cleanly discontinuous flux input is never possible in the real world. There will always be some sort of rounding or smearing at the edge of the hole. So pursuing an exact answer to a physically unrealizable setup doesn't make much sense. Try putting a ramp at the edge of the input flux distribution. |