| Author | Message | ||
     Kyle (kyle) Member Username: kyle Post Number: 15 Registered: 08-2006 |
Hi Robert: The diffusion equation in FlexPDE YCYLINDER coordinates can be written as: dt(c) = (1/r)*dr(r*D*dr(c)) + dz(D*dz(c)). I believe Natural(c) for planes normal to the z-axis should be defined as: Natural(c) = - (Flux in z-direction) = - (-D*dz(c)) = +D*dz(c). 1) Is the above correct? 2) What is Natural(c) for the radial direction? 3) What if the equation is re-written to an alternate, but equivalent form, dt(c) = (1/r)*(D*dr(c) + r*dr(D*dr(c))) + dz(D*dz(c)), where the radial component now has both first order and second order derivative terms. Does this change the form of Natural(c)? Thank you very much, Kyle | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 722 Registered: 06-2003 |
All equations are integrated over volume to form the finite element equations. In cylindrical geometry, this means multiplying by the volume differential r*dr*dz and integrating over r and z. The Divergence theorem says that the volume integral of a divergence is equal to the surface integral of the normal component of the argument of the divergence. This is the same result you get if you expand everything out an integrate by parts, which is what FlexPDE does. The Natural BC supplies the values of the surface integrands to allow FlexPDE to compute these terms. Applying the divergence theorem to Vol_Integral[Div(D*grad(c))] gives the surface term Surf_Integral(Normal(D*grad(c)). The Natural BC specifies the value of Normal(D*grad(c)). Therefore the value is the negative of the flux D*dn(c), regardless of the surface orientation. Actually, FlexPDE expands the divergence to (1/r)*dr(r*D*dr(c))+dz(D*dz(c)), so the form you have used is the same as using div(D*grad(c)). Multiplying by the r weighting of the volume element gives Integral[dr(r*D*dr(c))+dz(r*D*dz(c))]dr*dz This is converted by integration by parts to line_integral[(r*D*dr(c))*nr*dl + (r*D*dz(c))*nz*dl], the integral of the normal component of the (negative)flux. If you rewrite the equation as in your item 3), the r from the volume element cancels the leading (1/r) giving D*dr(c) + r*dr(D*dr(c))+r*dz(D*dz(c)) FlexPDE converts the second term to a pure derivative using the rule dr(f*g)=f*dr(g)+g*dr(f). This converts the equation back into the same form as the div(D*grad(c)) equation. The bottom line is that all three forms generate the same internal code. You can modify your script by adding SELECT DEBUG(FORMULAS) and FlexPDE will print out the final internal forms. In the printout, DX and DY are generic identifiers for derivatives wrt the first and second coordinates, not the Cartesian derivatives. |