Defining Curvilinear Coordinates for Anisotropic Materials

A lot of materials have anisotropic properties and, in many cases, the anisotropy follows the shape of the material. The COMSOL Multiphysics® software offers different methods for defining curvilinear coordinate systems. Here, we discuss the concepts of each and when to use which method.

Anisotropic Properties

Anisotropic properties are found in a wide variety of areas, such as rock formations with anisotropic seismic properties; liquid crystals used in LCD displays; materials for the aerospace industry, which must be lightweight and still withstand high loads; or soft tissues that biomedical replacements should mimic for optimal performance.

The Basics of Curvilinear Systems

In a previous blog post, we saw how to use the Curvilinear Coordinates interface and how to apply it to account for anisotropic thermal conductivity.

Let’s revisit this application and consider a carbon-fiber-reinforced polymer. The woven fibers embedded in an epoxy matrix have high thermal conductivity along the fiber axis and low conductivity in the cross section. It is almost impossible to express the anisotropy referring to the well-known Cartesian coordinate system. If we had a coordinate system that follows the fibers, it would be straightforward to specify the anisotropic properties.

Woven fibers in an epoxy matrix.

How can such a coordinate system be determined? Physically, there are numerous effects that result in a vector field following the shape of the geometry; for instance, flow through the fibers or heat conduction from one end to the other or even a bundle of current-carrying wires that produce a magnetic field. These are precisely the methods that are used in the COMSOL® software to compute the curvilinear system. All methods compute a vector field, \mathbf{v}, which forms the first basis vector. Since most applications require a normalized vector field, COMSOL Multiphysics automatically normalizes by dividing with |\mathbf{v}|. A second vector field is specified manually, and one of the Cartesian coordinates is often a good choice. Starting from this, the second basis vector \mathbf{e}_2 is reconstructed to ensure that it is perpendicular to \mathbf{e}_1 and normalized. The cross product of these two gives the third base vector, \mathbf{e}_3.

Internally, the software uses the Cartesian coordinate system (\mathbf{e}_x,\ \mathbf{e}_y,\ \mathbf{e}_z) for computation and converts all quantities referring to a different coordinate system (\mathbf{e}_1,\ \mathbf{e}_2,\ \mathbf{e}_3). A direction given by a vector, \mathbf{F}=(F_1,\ F_2,\ F_3), in an arbitrary coordinate system can always be transformed into Cartesian coordinates, as follows:

Here, \mathbf{M} is the transformation matrix. For the inverse transformation, simply use the inverse, \mathbf{M}^{-1}, and if (\mathbf{e}_1,\ \mathbf{e}_2,\ \mathbf{e}_3) is orthonormal then \mathbf{M}^{-1}=\mathbf{M}^{T}.

Next, I will illustrate the different methods available in COMSOL Multiphysics that can be used to calculate a curvilinear coordinate system, which include the:

• Diffusion method
• Adaptive method
• Flow method
• Elasticity method

Let’s pick a single fiber and have a closer look.

Diffusion Method

The diffusion method solves Laplace’s equation: -\Delta U=0. The solution, U, is a scalar potential, and its gradient forms the first base vector. Because you solve for a single scalar potential only, this method is computationally inexpensive. The direction of the vector field is specified with the inlet and outlet boundary conditions. If the geometry is a closed loop, you can set the Jump boundary condition on an interior boundary to specify the direction.

The diffusion method is equivalent to solving the stationary heat conduction equation with constant temperatures at the inlet and outlet boundaries. The temperature gradient then forms the first base vector, as illustrated below.

Curvilinear coordinate system (arrows), temperature gradient (streamlines), and temperature (surface).

Adaptive Method

The adaptive method, which is similar to the diffusion method, is based on Laplace’s equation as well. In addition, the resulting vector field is adapted to the shape of the geometry such that the streamline density is constant across the cross section of the geometry. This formulation is used to model multiturn coils (bundles of wires) in 3D magnetic applications with the AC/DC Module, an add-on to COMSOL Multiphysics. For multiturn coils, the current density should be roughly constant in the cross section, since it is assumed that each wire carries the same current and the wires are evenly spaced.

Flow Method

Here, you solve the incompressible Stokes equation for a vector field and a scalar. Thus, this method is the most computationally expensive. The boundary conditions are the same as for the diffusion method. A physical analogy would be incompressible creeping flow with a constant normal velocity at the inlet and a fixed pressure at the outlet. The resulting velocity field gives you the first base vector.

Curvilinear coordinate system (arrows), velocity field (streamlines), and pressure (surface).

Elasticity Method

The elasticity method solves the following eigenvalue equation:

where \mathbf{e} is the vector field, \mathbf{I} the identity matrix, and \lambdathe eigenvalue.

This method is slightly cheaper in terms of computational costs compared to the flow method because you solve for a vector field only. This difference in performance is more apparent in 2D models. The inlet and outlet boundary conditions are identical, \mathbf{e}\times \mathbf{n}=0. Prior to the implementation of the adaptive method, this method was used for modeling multiturn coils because it provided the best results in terms of constant streamline density in cross section.

Curvilinear coordinate system (arrows), coil direction (gray streamlines) and magnetic flux density (red streamlines).

Apart from these predefined methods, the COMSOL® software also provides a user-defined input, as usual. You may encounter other scenarios where you want to implement curvilinear coordinates manually, such as anisotropic hyperelastic material for modeling collagenous soft tissue in arterial walls.

Which Method Should I Choose?

At first glance, all of the methods lead to the same coordinate system. However, some shapes require special attention, and the choice of the method can lead to significantly different results when the coordinate system is used for a physical application. Pay attention to geometries that have at least one of the following characteristics.

Curvature

Let’s take a closer look at the streamlines for the various methods. Remember, the streamlines follow the vector field, which defines the first base vector. They start from equidistant points but follow different paths, as can be seen below.

Diffusion method: Streamlines follow the “shortest” path.

Adaptive method: Evenly distributed streamlines.

Elasticity method: Streamlines tend to accumulate at the convex bend.

Flow method: Streamlines tend to accumulate at the convex bend.

If bends are very sharp, the characteristics of each method are more distinctive, and also the adaptive method can develop nonuniform streamline density.

Variable Cross Section

In this case, the elasticity method can fail, and you obtain an eigenvalue with eigenvectors that do not produce the required coordinate system. Alternatively, you may have to manually search for the correct eigenvalue. We can see that the streamlines also do not follow the shape perfectly at the upper part of the geometry. Similar, but not as distinctive, behavior is observed with the diffusion and adaptive methods. The flow method provides the best results here but is also the most computationally expensive.

Diffusion method.

Adaptive method.

Elasticity method.

Flow method.

Streamlines along the center plane of the geometry.

Application to Heat Transfer

Let’s go back to our model with the fibers that have an anisotropic thermal conductivity of 60\ W/mK in the fiber direction and 4\ W/mK perpendicular to this. If these directions coincide with the axes of the coordinate system, the thermal conductivity — a second-order tensor — has 0 off-diagonal elements.

k=\left(\begin{matrix}
k_{xx} & k_{xy} & k_{xz} \\
k_{yx} & k_{yy} & k_{yz} \\
k_{zx} & k_{zy} & k_{zz}
\end{matrix}\right)=\left(\begin{matrix}
60 & 0 & 0 \\
0 & 4 & 0 \\
0 & 0 & 4
\end{matrix}\right)

To be able to use this diagonal form, the curvilinear coordinate system for the fibers must be calculated prior to solving for the heat transfer. Because the geometry doesn’t have sharp bends or a changing cross section, the diffusion method gives a fast solution for the curvilinear coordinate system.

After that, you can refer to this coordinate system in the Heat Transfer in Fibers node. The anisotropy of the thermal conductivity can be defined in the Materials node, with the syntax k=\{60, 4, 4\}. Alternatively, you can select the user-defined input in the associated heat transfer node.

Definition of anisotropy in the associated heat transfer node.

In the model, a boundary heat source in the form of a Gaussian pulse is applied to the center of the geometry and the temperature spreads along the fibers.

Streamlines indicate the vector field obtained with the Curvilinear Coordinates interface.

If you want to visualize, for example, the xx-component of the thermal conductivity (k_{xx}), keep in mind that you plot the xx-component in Cartesian coordinates (\mathbf{e}_x,\ \mathbf{e}_y,\ \mathbf{e}_z). Then, the thermal conductivity tensor, k, for the fibers is of nondiagonal form, according to the transformation described above. The local base vector system (\mathbf{e}_1,\ \mathbf{e}_2,\ \mathbf{e}_3), which was used to define k, is now varying through space and so does k_{xx}. In this model, you can plot the components of the thermal conductivity vector in a Slice plot, for example. You can select them from the Expressions menu in the corresponding Settings window or simply type ht.kxx (where ht is the tag for the Heat Transfer in Solids interface that was used for this model).

Concluding Remarks

This blog post has described the different methods for defining the curvilinear coordinate systems that are available in COMSOL Multiphysics as well as when to choose one over another.

In summary, the adaptive method provides the best solution for most cases at comparatively low computational costs. The diffusion method has even lower computational costs and is suitable for simple geometries without bends or varying cross sections. The other methods can be advantageous in some cases and are interesting for certain applications.

• Diffusion method
  • Pros: low computational cost
  • Cons: computed vector field tends to take the shortest path in bends
• Adaptive method
  • Pros: low computational cost, provides the best solution for most cases
  • Cons: variable cross sections not always handled perfectly
• Elasticity method
  • Pros: computational cost lower than the flow method, better representation of moderate bends than diffusion method
  • Cons: often requires manual selection of the eigenvalue and is not robust in all cases
• Flow method
  • Pros: robust method, supports cross-section changes and sharp bends
  • Cons: computational cost is often greater

To try the carbon fiber model presented here, click the button below.

Get the Tutorial Model

 

Editor’s note: This blog post was updated on December 5, 2018, to include information on the adaptive method.