Defining Curvilinear Coordinates for Anisotropic Materials

Nancy Bannach March 11, 2014
Share this on Facebook Share this on Twitter Share this on Google+ Share this on LinkedIn

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.
Anisotropic material fibers woven in an epoxy matrix
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.
Temperature gradients in a curvilinear coordinate system with anisotropic materials
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.
Velocity and pressure in a curvilinear coordinate system
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 showing the coil direction and magnetic flux density
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.

Different methods for modeling sharply curved geometries
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
Diffusion Method
Elasticity method
Elasticity method
Flow 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
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.

Using curvilinear coordinates for anisotropic materials
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

You can download the model presented here from our Model Gallery.



Comments

  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!

  3. Numi Rumi October 10, 2017   1:32 am

    Can we use curvilinear coordinates for time domain study?

  4. Nancy Bannach October 10, 2017   3:41 am

    Dear Numi,

    there should be no problem when you first solve the Curvilinear Coordinates (stationary) to get your coordinate system and then use this in a time dependent study for heat transfer, for example.
    A time changing coordinate system is not possible, as far as I know. If you are interested in something like this I recommend to discuss this with our support team in more detail via https://www.comsol.com/support

    Best Regards,
    Nancy

  5. Mst Nazmunnahar February 14, 2018   2:15 pm

    Hi,

    Thanks for this helpful blog. Can we use this anisotropic thermal conductivity under Heat transfer in solid(ht) sub-node solid 1 node or we need to add another node? Which node we should add for anisotropic thermal conductivity?

    Thanks

  6. Nancy Bannach February 15, 2018   11:11 am

    Dear Mst Nazmunnahar,

    you can also use it for solid, of course. Whenever you have the settings as shown in the image with the caption “Definition of anisotropy in the Heat Transfer node.”: An option to select a coordinate system and the option to define anisotropic material properties. So, you can use this approach in other physics as well.

    Best regards,
    Nancy

  7. Anonymous April 25, 2018   3:12 pm

    Dear Nancy,
    Hi! For me I applied creeping flow module with the porous media domains on, and I can only define the permeability of the porous media in “Fluid and Matrix Properties” interface and there is no option to choose the coordinate system like you shown to define the material property based on defined Curvilinear Coordinates as presented in the attachment. Could you help me to find a solution for this problem? Thanks a lot!

  8. Nancy Bannach April 26, 2018   2:59 am

    Hi,
    thanks for your feedback. You are right, there is no option to select another coordinate system for the Fluid and Matrix properties.
    The only option I see at the moment is, that you define the components manually. You can have a look at the model carbon_fibers.mph (https://www.comsol.com/model/anisotropic-heat-transfer-through-woven-carbon-fibers-16709).
    The definition of each component is written in the Equation View below the “Solid 2″ node. You can adapt this definition for your problem.
    If you need further assistance, please contact our Technichal Support team https://www.comsol.com/support.
    In the meantime, I will forward this to our developers.

    Best regards,
    Nancy

  9. Sergey Kuznetsov September 12, 2018   1:50 pm

    Hello!

    I was wondering where I can find the detailed comments in the creating Woven Carbon Fibers Geometry.

    The pdf document https://www.comsol.com/model/download/449671/models.heat.carbon_fibers_geom.pdf
    says that instructions are in the “Anisotropic Heat Transfer through Woven Carbon Fibers.”
    But the that document (https://www.comsol.com/model/download/469631/models.heat.carbon_fibers_infinite_elements.pdf) doesn’t have steps on creating geometry, instead it offers to open existin geometry “carbon_fibers_geom.mph”.

    Can i find anywhere actual instruction to create the geometry?

    Thank you!

  10. Nancy Bannach September 13, 2018   3:15 am

    Dear Sergey,

    you are right, there is no section on how to create the geometry. The sentence in the pdf file refers to the whole model set-up (not geometry set-up).
    If you are interested in how to create geometries in general I recommend to have a look at our video gallery
    https://www.comsol.com/videos?type%5B%5D=videotype-tutorial&workflow%5B%5D=workflowstep-geometry&s=

    Otherwise you can always have a closer look at the Geometry node and every single node that is used to build the geometry. E.g. when you right-click on “Work Plane 1″ in the geometry of the carbon fiber model you can choose “Build selected”. Then the geometry is build from the top node up to the Work Plane node.
    You can expand the Work Plane node to see, what is done there. Then go to “Sweep” and do the same and so on.

    I hope this helps. You can also contact our support team (https://www.comsol.com/support) in case you get stuck with your own geometry.

    Best regards,
    Nancy

  11. Sergey Kuznetsov September 15, 2018   9:54 pm

    Dear Nancy, thank you for explaining this!

Loading Comments...

Categorias


Tags