Guidelines for Equation-Based Modeling in Axisymmetric Components

Cylindrical coordinates are useful for efficiently solving and postprocessing rotationally symmetric problems. The COMSOL Multiphysics® software has built-in support for cylindrical coordinates in the axisymmetry physics interfaces. When defining custom partial differential equations (PDEs) using the mathematical interfaces, paying close attention to their meaning is important. The PDE interfaces assume partial differentiation in a Cartesian system, requiring manual coordinate transformations to change to a cylindrical system. See how to account for such coordinate transformations when using your own PDEs.

An Axisymmetric Heat Conduction Problem

Today’s example involves solving an axisymmetric heat conduction problem on a cylinder. This model was taken from the NAFEMS benchmark collection and solved in COMSOL Multiphysics using the Heat Transfer in Solids interface. We will take dimensions and material properties from that model and reproduce the result using the General Form PDE interface. Essentially we are now benchmarking our manual PDE interface approach to the problem by comparing it with the Heat Transfer in Solids interface, where the axial symmetry is automatically taken care of.

The equation for a stationary temperature distribution T on a rigid solid is

If the translational velocity u is uniformly zero, the heat capacity C_p has no effect and the only material property we have to specify is the thermal conductivity κ. Using the parameter kappa for thermal conductivity, we obtain the settings shown below with the Heat Transfer in Solids interface. The solution for a certain set of temperature and flux boundary conditions is also depicted.

An axisymmetric heat conduction problem solved with the Heat Transfer in Solids interface. This benchmark model is included in our Application Gallery.

Let’s now solve the same problem on the same finite element mesh using the General Form PDE interface. This interface is one of the several mathematics interfaces in COMSOL Multiphysics that facilitate solving custom equations. When the equation you want to solve is not built into one of the physics interfaces, you can use the mathematics interfaces to solve algebraic equations, ordinary differential equations, and partial differential equations. The new equations you add can be linear or nonlinear, and they can be solved either on their own or coupled with the prebuilt physics interfaces. The equations for heat transfer are a good starting point since those are already available in the built-in interfaces and we can easily compare the results with those where we have used our own PDE.

The General Form PDE interface has the template

and we have to provide the mass coefficient e_a, the damping coefficient d_a, the flux Γ, and the source f to specify our mathematical model.

If we compare this template with the equation from the Heat Transfer in Solids interface, with the dependent variable u standing for temperature, our flux specification should be that \Gamma = -\kappa \nabla u. To achieve this, we’ll specify the r– and z-components of the flux in the PDE interface as -kappa*ur and -kappa*uz, respectively. The mass and damping coefficients are not included in a stationary analysis.

As the image below indicates, the solution using the above settings differs from the solution we obtained using the Heat Transfer in Solids interface. What is the reason for this?

An axisymmetric heat conduction problem solved with the General Form PDE interface. The effect of a curvilinear coordinate system is not accounted for here.

The short explanation to such discrepancy is that the Heat Transfer in Solids interface understands \nabla \cdot (-\kappa \nabla T) as the divergence of the thermal flux, whereas the PDE interfaces understand \nabla \cdot \Gamma as

This expression is not the divergence of Γ in the cylindrical coordinate system. While we’ll demonstrate how to fix this in a later section, we’ll first provide a more detailed explanation of what causes the discrepancy.

Conservation Laws

The equations that the various physics interfaces in COMSOL Multiphysics solve are mathematical abstractions of the laws of physics. Often stated as conservation laws or accounting principles, these laws describe how a certain quantity changes on account of activities on a domain and across the domain’s boundary. While conservation laws hold for all materials, the extent of a given material’s response to domain forces and boundary fluxes differs from one material to another. The material responses are specified by so-called constitutive equations or equations of state. We have to make sure that constitutive equations do not violate the second law of thermodynamics. Together, conservation laws and valid constitutive equations provide enough information to derive a well-posed mathematical model.

For instance, thermal conduction in a rigid solid object is governed by the law of conservation of thermal energy. The rate of change of thermal energy equals the rate at which heat is supplied by sources in the domain plus heat flux through the boundary. From empirical observations, the heat flux in solids is proportional to the temperature gradient and is directed from hotter areas to colder areas. This is the constitutive equation. From here, we can use multivariable calculus to write the heat transfer equation for an isotropic material, for example, as

where T is temperature, which is the primary variable whose evolution we want to track; p, C_p, and K are material parameters describing density, heat capacity, and thermal conductivity, respectively; u is the translational velocity; and Q is the heat source per unit volume. Boundary conditions are also used to describe what happens on the boundary.

Coordinate Systems

After deriving mathematical models like the above equation, the next step is to solve them for the primary dependent variable and other quantities of interest. In heat conduction, for example, we want to obtain the temperature in time and space. To decide how we will identify a point in space, we select a coordinate system. Making a wise choice for a coordinate system can facilitate our analysis, whereas a bad choice can make our work unnecessarily complicated. No matter which coordinate system we choose, it is important to make sure that the physical meaning of the equation stays the same.

The heat conduction equation, for instance, contains the gradient \nabla T of temperature. In the Cartesian coordinate system, we have

whereas in the cylindrical coordinate system, we have

In the first case, partial derivatives of a scalar, with respect to independent variables, provide components of the gradient. This is not the case in a curvilinear coordinate system like the cylindrical coordinate system. As referenced earlier, the physical meaning of the gradient of temperature, as the vector that points in the direction of the greatest increase of temperature with a magnitude equal to the rate of increase, should stay the same. To ensure such invariance, we have to use covariant derivatives instead of regular partial derivatives. In the Cartesian coordinate system, covariant derivatives and partial derivatives overlap; in curvilinear systems, they do not.

The other differential operator we have is the divergence operator \nabla \cdot. In the heat transfer equation, this operator acts on the flux -\kappa \nabla T. When using the Cartesian coordinate system, taking the divergence of the flux involves differentiating only the components -\kappa \frac{\partial T}{\partial x}, -\kappa \frac{\partial T}{\partial y}, and -\kappa \frac{\partial T}{\partial z} of the flux. The basis vectors e_x, e_y, and e_z remain the same from one point to another. For the cylindrical and other curvilinear coordinate systems, the basis changes and taking the divergence thus involves taking derivatives of the basis vectors as well. Covariant derivatives take this into account. Simple sums of partial derivatives, when used in curvilinear coordinate systems, lose the physical meaning of the divergence reflecting how much a vector spreads out from a given point.

Similarly, the dot product \mathbf{u} \cdot \nabla T should always coincide with the product of the magnitudes of the two vectors and the cosine of the angle between them. In the Cartesian coordinate system, this is simply the sum of the products of the corresponding components of the two vectors. In curvilinear coordinate systems, the metric tensor of the coordinate system is added into the mix.

Because of the mathematical complexity that arises in curvilinear coordinate systems, you might wonder why we would want to use anything other than the Cartesian coordinate system. But, as we’ll highlight here, there are some applications where curvilinear coordinate systems are particularly useful.

When to Use Curvilinear Coordinate Systems

Consider a three-dimensional problem where the geometry, material properties, boundary conditions, and heat sources are symmetric about an axis. The solution is rotationally symmetric about that axis. If we use a coordinate system with the z direction along the symmetry axis, all partial derivatives with respect to Φ vanish. What that leaves us with is a two-dimensional problem in the rz-plane. This results in an equation that is easier to solve than the one in the Cartesian coordinate system, where all three spatial partial derivatives remain in the equation. In the finite element modeling of such problems, using an axisymmetric formulation facilitates the use of 2D meshes rather than 3D meshes, which leads to significant savings for both memory and time.

Curvilinear coordinate systems are also useful when material properties, such as the thermal conductivity κ, are not isotropic. If the anisotropy occurs in certain directions, we can define coordinate systems that align with those preferential directions and simplify the material property input. Note that COMSOL Multiphysics features a Curvilinear Coordinates interface that can be used to generate nonstandard coordinate systems. You can learn more about how to use the Curvilinear Coordinates interface and how to solve anisotropic problems on curvilinear coordinate systems on the COMSOL Blog.

For our particular example, we’ll focus on using the standard cylindrical coordinate system.

Physics Interfaces Vs. PDE Interfaces in Axisymmetric Components

In COMSOL Multiphysics, there are several physics-based interfaces that solve equations arising from one or more conservation laws. For instance, in the heat transfer in rigid solids example referenced above, we follow the conservation of thermal energy. In isothermal stress analysis problems, we follow the physical laws for conservation of mass, linear momentum, and angular momentum. When you use one of these physics interfaces, the software maintains the physical meanings of differential operators by using the corresponding expression for the Cartesian or curvilinear coordinates. That is, covariant derivatives are used instead of partial derivatives and, as a result, the coordinate system invariance is maintained.

But when you use the Coefficient Form PDE or General Form PDE interfaces, the software uses partial derivatives.

In other words, the physics interfaces understand \nabla u, \nabla \cdot \Gamma, and \nabla \times \Gamma to be the gradient, divergence, and curl, respectively, of a physical scalar u or higher-order tensor Γ. On the other hand, in the PDE interfaces, these tensorial meanings no longer apply and the operators are replaced by partial derivatives with respect to the independent variables. For example, the independent variables may not represent physical coordinates.

Consider \nabla \cdot \Gamma. In the physics interfaces, this means divergence. But in the PDE interfaces, this means [\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}]\cdot[\Gamma_x,\Gamma_y,\Gamma_z] in a 3D component and [\frac{\partial}{\partial r}, \frac{\partial}{\partial z}]\cdot[\Gamma_r, \Gamma_z] in a 2D axisymmetric component. The first results in \frac{\partial \Gamma_x}{\partial x} + \frac{\partial \Gamma_y}{\partial y} + \frac{\partial \Gamma_z}{\partial z}, which overlaps with the divergence, and the latter results in \frac{\partial \Gamma_r}{\partial r} + \frac{\partial \Gamma_z}{\partial z}, which is not the divergence of Γ.

In a cylindrical coordinate system, the divergence is given by

In an axisymmetric problem, the second item in this sum vanishes, leaving

(1)

The \nabla \cdot \Gamma expression in the PDE interfaces does not contain the last term, as the \nabla \cdot operator is interpreted simply as the sum of partial derivatives. For the divergence of a physical quantity in the cylindrical coordinate, we have to compensate for the last term in (1). We will use an example to highlight one approach for doing so using the source term.

Solution to Our Axisymmetric Heat Conduction Problem

Let’s now go back to our initial problem. What we want to solve is the stationary problem

where \nabla \cdot is the divergence operator. For an axisymmetric problem, we have

This is equivalent to

The left-hand side of this equation is what the General Form PDE interface in an axisymmetric component understands to be \nabla \cdot \Gamma. The right-hand side of the equation is added as an extra source term, as shown in the screenshot below. The solution now matches the solution that we obtained using the Heat Transfer in Solids interface.

An axisymmetric heat conduction problem solved with the General Form PDE interface. The effect of a curvilinear coordinate system is accounted for here.

Of course the item added in the source term is not a physical heat source. It simply compensates for the missing term between the covariant differentiation that we need for divergence and the partial differentiation the PDE interface does. A good practice is to add this term in the PDE settings window and add physical heat sources to the Model Builder by right-clicking the PDE interface and selecting “Source”.

In this example, only the divergence operator was used. For other differential operators, such as the curl, you should make similar compensations when solving a physics-based problem.

These adjustments are for curvilinear coordinate systems embedded in a Euclidean space. If the underlying space you are working with is not Euclidean, there are no Cartesian coordinates at all. In such cases, you need to be vigilant, even when using 2D or 3D components without axisymmetry. COMSOL Multiphysics features built-in physics interfaces for these types of problems. They include the Shell and Membrane interfaces in structural mechanics, the Thin-Film Flow interfaces in fluid mechanics, the Electric Currents, Shell interface in electromagnetics, and more. For a more extensive list, take a look at our product specification chart.

Using a Physics Interface as a PDE Template

Another approach to equation-based modeling is to use a physics interface as a PDE template for a problem with a similar mathematical structure. For example, to solve a physical problem that has a convection-diffusion-reaction nature, we can use either the Heat Transfer or the Transport of Diluted Species interfaces. These interfaces keep the tensorial meanings of differential operators in an axisymmetric component. All you need to do is to use the boundary and domain conditions that mathematically match the items you want to add.

The only downside to this strategy is that the units for the variables may not match the units you have. In such cases, you can use nondimensionalization. Once you obtain the dimensionless equations, you can make your COMSOL model dimensionless by going to the root of the Model Builder and setting Unit System to None.

With nondimensionalization, you can use a physics interface to solve a different type of problem with a similar mathematical structure but different dimensions.

Concluding Thoughts on Equation-Based Modeling in Axisymmetric Components

With the Coefficient Form PDE and General Form PDE interfaces in COMSOL Multiphysics, you can implement partial differential equations to solve novel problems not yet built into the software. These partial differential equations may or may not be derived from a physical problem. Therefore, the differential operators in the PDE interfaces are by design kept simple and not converted to tensorial operators automatically. When solving a physical problem using curvilinear coordinate systems — say when you want to exploit axisymmetry — make sure that the items you enter in the software’s templates accurately represent your physical problem. In the physics-based interfaces, COMSOL Multiphysics takes care of this. But in the PDE interfaces, the software can not attach any meaning to your equations, leaving the responsibility to you.

One way to address this is to add extra source terms to balance the difference between partial derivatives and covariant derivatives. Another way is to use an existing physics interface that has the same mathematical structure as your equation. If you have any questions about these strategies or other questions pertaining to this topic, please do not hesitate to contact us.

For more details on differential operators in curvilinear coordinate systems and partial differential equations on surfaces, you can turn to various books on tensor calculus or differential geometry. Here are some of my personal favorites, in no particular order:

• Rutherford Aris, Vectors, Tensors, and the Basic Equations of Fluid Mechanics, Dover Publications, Inc., 1989
• Pavel Grinfeld, Introduction to Tensor Analysis and the Calculus of Moving Surfaces, Springer Science+Business Media, 2013
• Mikhail Itskov, Tensor Algebra and Tensor Analysis for Engineers: with Applications to Continuum Mechanics, Springer-Verlag, 2013