FlexPDE treats all spatial coordinates on an equal footing, and tries to create meshes that are balanced in size in all coordinates.
Sometimes, though, there are problems in which one dimension is expected to have much less variation that the others, and fully meshing the domain with equilateral cells creates an enormous and expensive mesh. In these cases, it would be advantageous to scale the long dimension to bring the coordinate sizes into balance. Similarly, in semiconductor problems, for example, the structure is extremely thin, and would benefit from an expansion of the Z thickness coordinate.
It is possible that FlexPDE will eventually implement automatic coordinate scaling, but in the meantime, users can implement it manually.
Consider as an example the heat equation
div(k*grad(T))+Q = C*dt(T)
with k the conductivity, Q a source and C the heat capacity.
Define a coordinate transformation,
z = s*w
where w is the physical coordinate, z is the FlexPDE coordinate, and s is a scaling factor.
The expanded physical equation is then
dx(k*dx(T)) + dy(k*dy(T)) + dw(k*dw(T)) + Q = C*dt(T)
We can transform the heat equation using this transformation and observing that
dw(f) = (¶f/¶w) = (¶f/¶z)*(¶z/¶w) = s*dz(f)
The result is
(1) dx(k*dx(T)) + dy(k*dy(T)) + s*dz(k*s*dz(T)) + Q = C*dt(T)
In forming the finite element model for this equation, FlexPDE assumes continuity of the surface integrals generated by integration-by-parts of the second-order terms (equivalent in this case to the Divergence Theorem). This is the Natural Boundary Condition for the equation, as discussed elsewhere in the FlexPDE documentation.
The z-directed flux terms in the transformed equation therefore assume that s^2*k*dz(T) is continuous across cell interfaces. This is equivalent to flux conservation in the physical system as long as s is constant throughout the domain.
In order to guarantee conservation of flux in the presence of differing scale factors in layers, we must have the following equality across an interface between materials 1 and 2:
k1*dw(T)1 = k2*dw(T)2
k1*s1*dz(T)1 = k2*s2*dz(T)2
This will be satisfied if we divide our transformed equation by s :
(2) dx(k*dx(T))/s + dy(k*dy(T))/s + dz(k*s*dz(T)) + Q/s = C*dt(T)/s
where s is defined as s1 in material 1 and s2 in material 2.
Fluxes appropriate to the unscaled system can be recovered by the same modifications as those made in the PDE:
|•||Fluxes in the scaled direction must be multiplied by the scale factor. Integrals of these fluxes need not be further modified, as they are integrated over surfaces in true coordinates.|
|•||Fluxes in the unscaled directions are correctly computed in true coordinates, but when integrated over surfaces, they must be divided by the scale factor to account for the scaled area.|
Flux integrals then appear in the same form as in the scaled PDE:
Total_Real_Flux = Surf_Integral(NORMAL(-k*dx(T)/s, -k*dy(T)/s, -k*dz(T)*s)
Natural Boundary Conditions
The natural boundary condition defines the argument of the outermost derivative operator (or the argument of the divergence). In the conservative equation (2):
|•||Components in the unscaled direction have been divided by s. Therefore the natural boundary conditions for these components must be divided by s. (e.g. NATURAL(T) = x_flux/s on x-normal surfaces.)|
|•||In the scaled direction, the value defined by the natural is k*s*dz(T) which is in fact k*dw(T), the flux in the physical coordinate system. The natural in the scaled direction is therefore unmodified by the scaling. (e.g. NATURAL(T) = z_flux on z-normal surfaces.)|
"Samples | Usage | Coordinate_Scaling | Scaled_Z.pde" shows the implementation of this technique.
"Samples | Usage | Coordinate_Scaling | UnScaled_Z.pde" provides an unscaled reference for comparison.