Guidelines for Equation-Based Modeling in Axisymmetric Components

October 5, 2016

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

\rho C_p\mathbf{u}\cdot \nabla T + \nabla \cdot (-\kappa \nabla T)=Q,

If the translational velocity u is uniformly zero, the heat capacity Cp 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.

Solving the axisymmetric heat conduction problem with the Heat Transfer in Solids interface in COMSOL Multiphysics.
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

e_a\frac{\partial^2 u}{\partial t^2} + d_a\frac{\partial u}{\partial t} + \nabla \cdot \Gamma = f,

and we have to provide the mass coefficient ea, the damping coefficient da, 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?

Screenshot showing the use of the General Form PDE interface, without including the effect of a curvilinear coordinate system.
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

\frac{\partial \Gamma_r}{\partial r} + \frac{\partial \Gamma_z}{\partial z}.

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

\rho C_p\frac{\partial T}{\partial t}+\rho C_p\mathbf{u}\cdot \nabla T + \nabla \cdot (-\kappa \nabla T)=Q,

where T is temperature, which is the primary variable whose evolution we want to track; p, Cp, 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

\nabla T = \frac{\partial T}{\partial x}\mathbf{e}_x + \frac{\partial T}{\partial y}\mathbf{e}_y+\frac{\partial T}{\partial z}\mathbf{e}_z,

whereas in the cylindrical coordinate system, we have

\nabla T = \frac{\partial T}{\partial r}\mathbf{e}_r + \frac{1}{r}\frac{\partial T}{\partial \phi}\mathbf{e}_{\phi}+\frac{\partial T}{\partial z}\mathbf{e}_z.

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 ex, ey, and ez 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.

\nabla \cdot (\kappa \nabla T) = \frac{\partial }{\partial x}(\kappa \frac{\partial T}{\partial x}) + \frac{\partial }{\partial y}(\kappa \frac{\partial T}{\partial y}) + \frac{\partial }{\partial z}(\kappa \frac{\partial T}{\partial z}) = \frac{1}{r}\frac{\partial }{\partial r}(r \kappa \frac{\partial T}{\partial r}) + \frac{1}{r}\frac{\partial }{\partial \phi}( \frac{\kappa}{r}\frac{\partial T}{\partial \phi})+\frac{\partial }{\partial z}(\kappa \frac{\partial T}{\partial z}).

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

\nabla \cdot \Gamma = \frac{1}{r}\frac{\partial (r \Gamma_r)}{\partial r} + \frac{1}{r}\frac{\partial \Gamma_{\phi}}{\partial \phi}+\frac{\partial \Gamma_z}{\partial z}.

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

(1)

\nabla \cdot \Gamma = \frac{1}{r}\frac{\partial (r \Gamma_r)}{\partial r} +\frac{\partial \Gamma_z}{\partial z} = \frac{\partial \Gamma_r}{\partial r} + \frac{\partial \Gamma_z}{\partial z}+\frac{\Gamma_r}{r}.

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

\nabla \cdot \Gamma = f,

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

\frac{\partial \Gamma_r}{\partial r} + \frac{\partial \Gamma_z}{\partial z}+\frac{\Gamma_r}{r}=f.

This is equivalent to

\frac{\partial \Gamma_r}{\partial r} + \frac{\partial \Gamma_z}{\partial z}=-\frac{\Gamma_r}{r}+f.

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.

Screenshot depicting the correct General Form PDE 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.

Screenshot showing the use of nondimensionalization.
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

Comments (15)

Leave a Comment
Log In | Registration
Loading...
Ivar Kjelberg
Ivar Kjelberg
October 6, 2016

Thanks for a interesting and nicely illustrative BLOG 🙂

For me you give us here a few important lessons:
1) for the Engineers: buy the modules with the physics done by COMSOL, globally it’s cheaper for you than making the modules from scratch for your need, with the risk of getting it wrong, and you need to re-verify and validate your own modules for each new release. Put down the true effort number and show them to your boss, then he might/will understand 😉
2) for the Students: do the math/physics yourself once or twice, and compare well to the COMSOL modules, this is the best way for self learning and verifications.
3) If in doubt: use the Forum and COMSOL Support 🙂

If I could add a wish, for COMSOL Multiphysics: get it ready once to allow us to enter directly full tensor math as well as for the postprocessing, this would gain time compared to now when we need to rewrite everything for real/imaginary and only scalar components one by one.

Sincerely
Ivar

Eric Monsu Lee
Eric Monsu Lee
October 13, 2016

Good article Teme,

Without the coordinate compensation, the general/coefficient PDE interface will usually not converge and even if it converges, the solution is always incorrect when comparing to an analytical solution.

Krishna Kumar Soundararajan
Krishna Kumar Soundararajan
May 12, 2017

This article is truly useful. Then, how about the convection term? Should we change ux*d(T,x)+uy*d(T,y) into ur*d(T,r)+uz*d(T,z)+ur*T/r or ur*d(T,r)+uz*d(T,z)? Thanks a lot.

Temesgen Kindo
Temesgen Kindo
June 21, 2017

Dear Krishna,
Thank you for your comment and question.

If you are using the general form PDE, the convection term contributes u*T to \Gamma. As a result in \Gamma_r there will be a convection contribution of u_r*T. So you can add a source term -u_r*T/r.

When using the coefficient form you have two convection contributions. The conservative part should be handled just like what we discussed above in the general form PDE. The non-conservative one is \beta \cdot \grad T. In an axisymmetric problem the angular part does not contribute to this term and you do not have to do anything about it.

Cheers,
Temesgen

Xiang Sun
Xiang Sun
August 31, 2018

Very useful article. So, how to set the flux boundary for this problem. Just as -n_r*Gamma_r-n_z*Gamma_z? Thanks a lot.

Kostis Lazaridis
Kostis Lazaridis
September 18, 2018

This is a very useful article.
I have a question regarding the boundary conditions on the axis of symmetry in general.
What are the default boundary conditions on the axis of symmetry (I see in this case the axis of symmetry is not included on the domain). I’ve noticed that even if a boundary condition is not prescribed on the axis of symmetry, comsol doesn’t set zero flux on the axis. But the simulation can start. So, what is the default boundary condition on the axis of symmetry?

Thanks a lot
Kostis

Temesgen Kindo
Temesgen Kindo
September 19, 2018

Thanks Kostis for your comments.

After you enter your equation in the strong form, the software will convert it to the weak form before solving the latter. The weak form contains surface integrals of test(u)*(n dot Gamma). So leaving a surface alone is equivalent to applying the zero flux boundary condition. See the following blog post for more on weak forms.
https://www.comsol.com/blogs/brief-introduction-weak-form/
In the built-in physics interfaces we understand what axisymmetry means for each physics. With user defined equations, however, we can not predict what axisymmetry means and as such it can be inaccurate for us to force a certain boundary condition on the axis of symmetry. Thus we leave it alone but as explained above leaving it alone is the same as zero flux for advection-diffusion type equations.

Regards,
Temesgen

Sumith S
Sumith S
March 5, 2019

This is a very useful article.
I have a question regarding the slip boundary conditions on the axis of symmetry in Coefficient Form PDE
How Can I impose the Slip boundary Condition in Coefficient form PDE, either by apply as a flux souce or by using constraint, which one is more effective?

Brianne Christopher
Brianne Christopher
March 11, 2019

Hello Sumith,

Thank you for your comment.

For questions related to your modeling, please contact our Support team.

Online Support Center: https://www.comsol.com/support
Email: support@comsol.com

Eric Williams
Eric Williams
July 14, 2019

Is the statement true about cylindrical differences with version 5.4? I see a different equation in the coefficient form PDE that includes several additional terms which should account for the actual divergence instead of the cartesian only as you indicate. Can you please clarify if this is still accurate?

Wasiu Abidemi Hamzat
Wasiu Abidemi Hamzat
February 14, 2024

Your link not opening

Thomas König
Thomas König
August 3, 2022

The behavior of solving the PDEs on Cartesian coordinates when the user has specified cylindrical symmetry is highly dangerous – it generates the worst kind of error, silent bad results which look plausible at first glance.

I fell into this trap, and it was only through re-checking that I found that I had solved the wrong equation.

The least I would expect is a prominent warning that Comsol is doing the wrong thing here.

Kaleem Iqbal
Kaleem Iqbal
August 16, 2022

Hi, How do we give a sliding wall velocity in coefficient form PDE or General form PDE? Thank you

Sayali Jadhav
Sayali Jadhav
December 30, 2022

Thank you for the blog. I found the explanation regarding using the source term in axisymmetric case very useful. Regarding the same, can you suggest me how to write derivatives of interface normal in the axisymmetric models, while using a general form PDE on a boundary?
I have posted this query in detail here – https://www.comsol.com/forum/thread/318031/how-to-define-derivatives-of-interface-normal-in-moving-mesh-model.

Thanks a lot 🙂

EXPLORE COMSOL BLOG