| Author | Message | ||
     Britta Hagmeyer (britta) New member Username: britta Post Number: 1 Registered: 11-2006 |
Hello, I have a problem with setting up a particle flux into a 3D cylinder through a limited region by a natural boundary condition... I'm trying to create a 3D-simulation of the time-dependent diffusion of particles from one compartment (comp1) through a tunnel into another compartment [comp2). The compartments and the tunnel are being modelled as cylinders. The particle concentration in comp1 is u_tr, the concentration in comp2 and the tunnel is u = 0 at t = 0. The system variable is u. The diffusion of the particles from comp1 through the tunnel should follow Fick's 1st law (Flux = -D*(grad(u)) ) and the diffusion of the particles in the second compartment comp2 should follow Fick's 2nd law (dt(u) = div(D*grad(u)), with D being the diffusion constant. I've already tried several system setups: 1. One cylinder with three compartments, being comp1, the tunnel (realised by a limited region in a third compartment of the cylinder) and comp2. Here, the mesh around the tunnel is being constructed with a very high resolution... how do I define my system so that this doesn't happen? 2. One cylinder with just one compartment, being comp2, and the tunnel and comp1 being replaced by a limited region in the bottom of the cylinder. Therefore, the particle flux into comp2 through the limited region, which I set using a natural boundary condition, should be following Fick's 1st law. Which system setup is better for my problem? I've been working on the second system setup. I've set the boundary conditions to natural=0 at all sides of the cylinder and now I need a boundary condition for the flux through the limited region. But how do I set the flux of particles according to Fick's 1st law into the system through the limited region? Which boundary condition should I use? What does it look like? I've already tried several conditions: natural(u) = -D*dz(u_tr) natural(u) = -D*normal(grad(u_tr)) natural(u) = -D*(u_tr - u) There doesn't seem to be any substance flowing into the cylinder through the pore... where is the problem? What have I been missing or misunderstanding? Should I use a value boundary condition?? I've already read the help manual, but I cannot see the problem... I hope you can help me!! Here is the program code: title 'Diffusion of fluids in water at 25°C' coordinates cartesian3 variables u(threshold=0.1) definitions D = 5.34e2 { diffusivity micron^2/sec Eosin in water} { variables for the size of the tunnel(pore), the radius and the height of the mesh } r_pore = 0.025 { Pore radius in micron } meshheight = 1 { meshheight } r_mesh = 0.56 { Radius of the mesh } V_tr_t0 = 0.25 { volume of comp1 } u_tr_t0 = 200 { Concentration of particles at the beginning of the simulation in comp1 } n_tr_t0 = u_tr_t0 * V_tr_t0 { Number of particles at the beginning of the simulation in comp1 } n_cyl = u*pi*r_mesh*r_mesh*meshheight { number of particles in comp2 = concentration in comp2 * comp2 volume } n_tr = n_tr_t0 - n_cyl { number of particles in comp1 = particle concentration in the beginning of the simulation in comp1 - number of particles in comp1} u_tr = n_tr/V_tr_t0 timesteps = 2e-1 { time steps of the simulation } bc = -D*(u_tr - u) !-D*dz(u_tr) -D*normal(grad(u_tr)) { for the natural boundary condition } initial values u = 0 equations div(D*grad(u)) = dt(u) { Diffusion according to Fick's 2nd law } extrusion { 3D-Structure } surface 'fluid_bottom' z = 0 layer 'fluid' surface 'cylinder_top' z = meshheight { Creating the cylinder and the limited region } boundaries surface 'cylinder_top' natural(u) = 0 { Isolation for 'cylinder_top' } region "cylinder" start(0,r_mesh) natural(u) = 0 { Isolation for the sides of the cylinder } arc(center=0,0) angle 360 limited region "bottom_pore" surface 'fluid_bottom' natural(u) = bc !10*(1-u) { Refining 'bottom_pore' to 'fluid_bottom'. } { boundary condition.: Flux at 'fluid_bottom' in region 'bottom_pore' } start(0,r_pore) arc(center=0,0) angle 360 time 0 to 1 by timesteps plots for t = 0 by 0.4 to endtime contour(u) on x = 0 surface(u) on z = 0 elevation(u) from (0,0,0) to (0,0,0.5) Thanks a lot! Britta | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 708 Registered: 06-2003 |
1. Fick's Law is automatically satisfied both internally and at internal material interfaces. The Natural BC is computed by integrating the second-order terms of the PDE by parts. The surface terms that arise define the meaning of the Natural BC. In your case, this is the outward normal component of D*Grad(u). At an internal material interface, this value will be continuous across the interface. In other words, flux is conserved. 2. This definition of the Natural BC means that (for your model 2) at the inlet spot you want to define Natural(u) as the net flux entering the cylinder (ie, the outward normal component of D*grad(u)). If the external cylinder has concentration U1, then your (commented out)BC 10*(1-u) is a reasonable one, and the "10" represents D0/deltaz, the diffusivity in the outer chamber divided by an estimated fall-off distance. 3. To model it more correctly, you can use two cylinders with the transfer opening only on the interface surface. You do not need to specify any boundary condition on the internal face; Fick's law is the default continuity condition. 4. If you construct the two-chamber model with a discontinuous initial condition (eg, 1 in the outer and 0 in the inner chamber), you will get a massive mesh refinement on the transfer port. FlexPDE is trying to represent the physically impossible discontinuity. It will do the same with the timestep, cutting it to tiny values while calculating the motion of the infinite slope. Don't worry, these things will go away after the front has diffused a little. Or, you can specify an initial condition that is not discontinuous; it's up to you. 5. I don't know what you are doing with your u_tr. You have effectively defined the source concentration as a function of the interior concentration. I don't know what this means. | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 709 Registered: 06-2003 |
PS. You can do this problem more cheaply in 2D cylindrical geometry. | ||
| Britta Hagmeyer (britta) New member Username: britta Post Number: 2 Registered: 11-2006 |
Thank you for the advice! I'm afraid I have another thing to ask: I've now built a system with two cylinders and a transfer opening in the interface surface between them. But how do I get this interface surface to be impermeable to the diffusion substance, except for the transfer opening in the middle? I've tried VALUE(u) = 0 on the interface surface, but this doesn't seem to work... Thanks! Britta | ||
| Britta Hagmeyer (britta) Junior Member Username: britta Post Number: 3 Registered: 11-2006 |
Thank you for the advice! I'm afraid I have another thing to ask: I've now built a system with two cylinders and a transfer opening in the interface surface between them. But how do I get this interface surface to be impermeable to the diffusion substance, except for the transfer opening in the middle? I've tried VALUE(u) = 0 on the interface surface, but this doesn't seem to work... Thanks! Britta |