## Defining Curvilinear Coordinates for Anisotropic Materials

##### Nancy Bannach March 11, 2014

A lot of materials have anisotropic properties and, in many cases, the anisotropy follows the shape of the material. COMSOL Multiphysics 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:

\mathbf{F}=\left(\begin{matrix} F_x\\F_y\\F_z\end{matrix}\right)=F_1\mathbf{e}_1+F_2\mathbf{e}_2+F_3\mathbf{e}_3=\left(\begin{matrix}
e_{11} & e_{12} & e_{13} \\
e_{21} & e_{22} & e_{23} \\
e_{31} & e_{32} & e_{33}
\end{matrix}\right)\cdot\left(\begin{matrix} F_1\\F_2\\F_3\end{matrix}\right)=\mathbf{M}\left(\begin{matrix} F_1\\F_2\\F_3\end{matrix}\right)

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}.

The following are available in COMSOL Multiphysics:

• Diffusion method
• Flow method
• Elasticity method

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

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 compared to the other two. 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).

#### 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:

-\nabla\cdot((\nabla\cdot\mathbf{e})\mathbf{I}+(\nabla\mathbf{e}+\nabla\mathbf{e}^{T}))=\lambda\mathbf{e}

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. The eigenvalue formulation was devised initially to model multi-turn coils (bundles of wires) in 3D magnetic applications with the AC/DC Module. For Multi-turn coils, the current density should be roughly constant in cross section, since it is assumed that each wire carries the same current and the wires are evenly spaced.

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 three methods lead to the same coordinate system. Most commonly, you want to choose a method that incurs low computational costs. In that case, the best choice is the diffusion method. That said, some shapes require special attention and the choice of the method can lead to significantly different results when the coordinate system is used by a physical application.

#### Strong Curvature

Let’s have a look at a geometry with a sharp bend, plotted below. If we look on the right-hand side of the plot, from top to bottom, we see that with the diffusion method, the streamlines accumulate at the inner radius because they take the shortest path. The elasticity method shows some wobbly curves, while the flow method provides a solution where the streamlines follow the expected path.

Results of different methods for sharply curved geometries.

#### 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 behavior is observed with the diffusion method. The flow method provides the best results here, but is most computationally expensive.

Diffusion 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 2nd-order tensor — has zero 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 associated Heat Transfer 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 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 non-diagonal 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 has described the different methods for defining the curvilinear coordinate systems that are available in COMSOL as well as when to choose one over another.

• Diffusion Method
• Pros: low computational cost
• Cons: computed vector field tends to take the shortest path in bends
• 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: most robust method, supports cross section changes and sharp bends
• Cons: computational cost often greater

#### Post Tags

Technical Content

1. Ivar Kjelberg February 17, 2015   12:20 pm

Hi Nancy
Thanks for å detailed discussion, I have one issue though:
I need to consider aniso-tropic thermal expansion, but I have no coordinate frame selection list to choose from in the thermal expansion sub-node, in contrary to your example of the thermal conduction k.
So what are the correct equations & variable names to use in this case? V5.0

Sincerely
Ivar

2. Yangguang Ou December 16, 2015   3:40 pm

Thank you for the helpful blog!