What is a Superellipse?

An ordinary ellipse is commonly described by the equation

where a and b are called the semiaxes of the ellipse. When a = b, the equation of a circle with radius a is recovered.

The superellipse is a generalization, where an exponent other than 2 is used, so that

The fountain in the center of Stockholm uses n = 2.5 and the shape factor a/b = 6/5.

The fountain at Sergels Torg (left) and an outline of its shape (right).

Examples of other superellipses are shown below.

You can find more details about the mathematics of superellipses here.

Stress Concentration Factor

The stress concentration factor, K_t, is an important concept in solid mechanics. It is used to describe the increase in stress around a geometric discontinuity in a structure. K_t relates the maximum stress to a suitably defined nominal stress (stress without the geometrical discontinuity) through

Traditionally, stress concentration factors for common cases have been tabulated in handbooks. With today’s easy-to-use finite element (FE) programs, computing a stress concentration factor is often faster and more accurate than looking it up in a graph or table.

The most well-known result is likely that K_t = 3 for a circular hole in a uniaxially loaded large plate.

A model in COMSOL Multiphysics^® of von Mises equivalent stress around a circular hole in a plate subjected to a unit axial stress in the horizontal direction.

Another analytical result is the stress concentration factor for an elliptic hole in a large plate:

Here, the semiaxis with length b is the one perpendicular to the stress field so that the stress concentration increases with the ratio b/a.

A model of von Mises equivalent stress around an elliptical hole with b/a = 3 in a plate subjected to a unit axial stress in the horizontal direction.

Stress Analysis of the Superelliptic Hole

A superelliptic curve can easily be constructed in the COMSOL^® software using the Parametric Curve feature.

A superelliptic curve in the Model Builder made using the Parametric Curve feature settings.

Here, a representation of the superellipse in polar coordinates is used for the parameterization. It would also be possible to use a simpler parameterization based on the original equation expressed in x and y, but it becomes less accurate, since the relation between yand x is highly nonlinear.

The parameterization used (left), and an alternative, simpler but less efficient, parameterization (right).

First, let’s take a look at the results for a hole with a geometry similar to the aforementioned fountain. The peak stress depends on the orientation of the larger semiaxis, but in both cases the stress concentration factor is lower than the value of 3 for a circular hole.

Models of von Mises equivalent stress in a plate with a superelliptic hole having n = 2.5 and a/b = 6/5.

This looks promising. From here on, it is easy to set up a parametric sweep and study the stress pattern for many different shapes.

For any value of the superellipse exponent n less than 2, there will be sharp corners on the hole edge. That would cause stress singularities. Such cases are not of interest in this context. Thus, n is kept in the range 2–8. At the highest values, the hole is almost rectangular with corner fillets. For the axis ratio b/a, values ranging from 0.2 to 5 are tested. The computed stress concentration factors are shown in the diagram below.

Stress concentration factors for superelliptical holes for a range of exponents n and axis ratios q = b/a. The markers on the curves show the minimum values.

For n = 2, it can be seen that the values are as expected for an ellipse:

If the curves are normalized using this factor, we can see how a superellipse differs from an ordinary ellipse having the same axis ratio:

Normalized stress concentration factors for superelliptical holes.

As can be seen, there is always a superellipse with the same axis ratio as a certain ellipse that will give a smaller stress concentration factor. In particular, if we replace a circle with the best possible symmetric superellipse (q = 1), then it is possible to reduce the stress concentration factor by 14% using n = 3. This is quite a significant improvement. Such a decrease in stress could improve the fatigue life by a factor of 2.

The optimal replacement for a circular hole.

Strengthening Structures While Increasing Hole Area

The area enclosed by a superellipse can be expressed as

where Γ is the gamma function.

For comparison, the area of an ordinary ellipse is

This means that the relation between the area of a superellipse and an ordinary ellipse having the same semiaxes is independent of the axis ratio and can be expressed as

In the plot below, the function is shown. 

Relative area increase as a function of n when compared to an ellipse. The asymptotic value 4/π is indicated by the dashed line.

It can be seen that the area of any superellipse with n > 2 is always larger than that of the ellipse with the same semiaxes. Somewhat surprisingly, this means that it is always possible to lower the stress concentration factor by changing a circular or elliptic hole to a superelliptic one that has a larger area, that is, by removing material. In my previous blog post, Making Structures Stronger by Removing Material, some other cases where material removal is beneficial are presented.

It should, however, be noted that we have only investigated uniaxial stress states that are aligned with one of the semiaxes of the hole. For some other stress states, a stress reduction effect cannot be obtained. The reduced radii of the corners of the superellipse will instead raise the stresses.

Real-World Applications of Superellipses

This is, of course, mainly a fun theoretical discussion. In practice, it is often much easier to drill circular holes. But if, for example, additive manufacturing or casting is used, then it can be beneficial to choose other shapes.

The approach of using a superelliptic shape is not only applicable to holes. The same idea can be used to reduce stress concentrations at fillets.

You can download the model used in the examples above by clicking the button below.

Precursors to Glasses

Reading Stones

One of the earliest examples of people using a lens to assist vision came not in the form of glasses but in the form of a small quartz stone. Reading stones, as they were dubbed, are small, hemispherical stones carved out of quartz, beryl, or glass and shaped in such a way to function as a convex lens. When placed upon paper, this lens would magnify the text written on it for those whose vision was failing. The invention of these stones is often attributed to Ibn Sahl, a 9^th-century Persian mathematician, but they have been found all across the world, including in Viking graves in Sweden.

A glass reading stone in Archeon historical theme park in the Netherlands. Image licensed under the Creative Commons Attribution-Share Alike 4.0 International license via Wikimedia Commons.

Sunglasses

Another early precursor to the corrective glasses we know today? Sunglasses, which were developed by multiple cultures around the world, some of the earliest of which were found in 12^th-century China and featured flat panes of smoky quartz. Another example of early noncorrective eyewear is traditional Inuit snow goggles. These goggles, typically carved from antler, driftwood, or ivory, reduced exposure to sunlight and helped prevent glare from the reflective snow.

Two pairs of Inuit snow goggles from Alaska. Image licensed under the Creative Commons Attribution-Share Alike 4.0 International license via Wikimedia Commons.

Early Corrective Lenses

Most sources attribute the invention of the first corrective lenses to late 13^th-century Italy, although their inventor is unknown. Early evidence lines up with this, as paintings and sermons involving lenses began appearing around this time throughout Europe. There is further evidence of rules and regulations having been set in place for early lens makers by the 14^th century, during the Venetian Renaissance. These “glasses” looked a good deal different than those that exist today; they were often held by hand, made of two magnifying glasses attached to a pair of riveted handles to form an almost scissor-like shape.

A pair of French “scissors glasses”. Image licensed under the public domain via Wikimedia Commons.

How Do Glasses Work?

How Do We See?

Before looking into how glasses work, it helps to first understand vision itself. Seeing is a complex process that takes place in milliseconds. At a high level, the eye can be loosely thought of as two focusing lenses that project incoming light into an image on the retina. The cornea, which is the outer layer of the eye, refracts the incoming light toward the crystalline lens, which in turn focuses the light onto the retina. The retina is where the light is absorbed into our eye’s photoreceptors, the rods and cones. These photoreceptors then send electrical signals to our brain, which can then interpret what we see. To focus on objects at different distances, tiny muscles deform the crystalline lens to change its focus. The last crucial elements are the iris and pupil, which are between the cornea and the crystalline lens. The iris is muscular tissue in the shape of an annular disc, and the pupil is the hole in its center. The iris expands or shrinks the size of pupil to adjust the amount of light that reaches the retina.

Issues with vision can happen during any part of this process but most commonly stem from natural imperfections in eye shape. These include things like the eyeball being misshapen (e.g., excessively oblong), the cornea itself being uneven or not convex enough, or the crystalline lens being damaged. As we age, the crystalline lens also becomes stiffer, causing presbyopia.

How Do Glasses Help?

Glasses use lenses to bend beams of light before they reach the eye, refocusing them based on the eye condition. In the case of myopia, or nearsightedness, the light naturally comes to focus before hitting the retina, causing objects that are far away to appear blurry. Lenses with a concave shape help send the focal point of the light farther away from the retina. In cases of hyperopia, or farsightedness, the opposite happens, so a convex lens is used to move the focal point of the light closer to the retina. When the cornea is uneven, also known as having an astigmatism, corrective lenses are shaped in more complex ways.

Top: By moving the focal point of light beams closer to the retina, a convex lens corrects farsightedness. Bottom: By moving the focal point farther from the retina, a concave lens corrects myopia. Image licensed under the Creative Commons Attribution-Share Alike 4.0 International license via Wikimedia Commons.

How Have Glasses Evolved?

Bifocals

The invention of bifocals is often attributed to Benjamin Franklin, who used them in court to read the lips of French speakers while also taking and reading notes. While Franklin is rumored to have fashioned his bifocals by sawing in half the lenses of two pairs of his glasses and then attaching them together in a new frame, modern bifocals are made by molding a high magnification segment onto a regular corrective lens.

A pair of bifocal glasses. Image licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license via Wikimedia Commons.

Progressive lenses take bifocals a step further, providing three or more magnification levels in a single pair of glasses and allowing wearers to see at multiple distances.

Modern Lens Technology

By the late 19^th century, glass had long since replaced costly quartz and beryl as the primary material for lenses, and in Germany, physicist Ernst Abbe and chemist Otto Schott experimented with adding other elements into the melted glass to increase factors like durability, clarity, and weight. In the 21^st century, many lenses are made of specialized plastics that further reduce cost and weight while increasing durability.

Additional modern developments include blue-light, scratch-resistant, anti-glare, and photochromic technologies that involve coating and/or impregnating lenses with various materials. Anti-glare lenses were originally developed for military optics during World War II but were made available for public use in the late 1940s. These coated lenses help glasses to not appear distorted in photographs and improve visual acuity in situations with inconsistent or bright light, such as driving at night or viewing a computer screen. Photochromic lens material (which may be applied as a coating or within the lens) is UV reactive and darkens in sunlight, removing the need for a person to own a pair of prescription sunglasses and glasses.

An uncoated pair of glasses (top) compared to one with an anti-glare coating (bottom). Notice the reflection of the photographer in the top frame and the blue-tinted lights reflected in the bottom pair. Image licensed under the Creative Commons Attribution-Share Alike 3.0 Unported via Wikimedia Commons.

Contact Lenses

Compared to glasses, contact lenses offer better peripheral vision and are less impacted by weather conditions. The precursor to contact lenses was proposed in the 16^th century by Leonardo da Vinci, who suggested wearing water-filled glass hemispheres over the eyes. The first contact lenses, developed in the 19^th century, were rigid and made of glass and wax. It wasn’t until the early 1960s that modern soft lenses were invented.

Looking Forward

While we’ve mentioned a few recent developments in corrective eyewear here, other innovations are on the horizon. For example, researchers have used simulation to develop liquid lenses that dynamically change optical power and reduce eye fatigue when using AR/VR devices. Others are using virtual modeling to investigate lighter and more durable frames and specialized coatings for glasses. As advancements to eye care continue to be made and wearable tech devices like smart glasses come closer to everyday availability, we only have to wait and see where this technology will take us in the future.

To drive innovation in eye care technology and other optical systems, many engineers turn to ray optics modeling. The Ray Optics Module, an add-on to COMSOL Multiphysics^®, enables the modeling of electromagnetic wave propagation through ray tracing. With this approach, propagating waves are represented as rays that can reflect, refract, or be absorbed. Learn more about the module via the button below.

Reference

1. O. Wilson, “What Percentage of People Wear Glasses? The Global Vision Correction Reality,” GLASSON, 18 Dec. 2025; https://www.glasson.app/blog/what-percentage-of-people-wear-glasses-the-global-vision-correction-reality/.

In dielectric materials with ultrahigh thermal conductivity, heat can violate diffusion and behave fluid-like, forming vortices, displaying temperature waves, and even locally backflowing against the temperature gradient. In this blog post, we introduce a phonon hydrodynamics custom physics interface that we developed for the COMSOL Multiphysics^® software and show how continuum modeling can capture viscous heat transport phenomena beyond Fourier’s law.

Why Model Heat Beyond Fourier’s Law?

Fourier’s law has long been the standard framework for describing heat conduction as it successfully captures thermal transport in common materials and devices at room temperature and at millimeter or larger scales. However, recent experimental and theoretical studies have shown that Fourier’s picture can break down in microscale devices made of ultrapure, high–thermal-conductivity materials such as diamond, graphite, and hexagonal boron nitride.

In these systems, phonons — the primary heat carriers — undergo frequent momentum-conserving collisions within certain temperature ranges, which can extend from cryogenic (~70 K) up to near room temperature. This regime is called phonon hydrodynamics. As a result, heat transport becomes collective and fluid-like, leading to phenomena such as temperature waves (“second sound”), Poiseuille-like heat flow (i.e., faster in the center of a channel and slower at the boundaries), thermal backflow, and heat vortices.

Capturing these effects requires a model that goes beyond standard Fourier diffusion while remaining practical for realistic geometries. This is precisely the motivation behind the viscous heat equations (VHE) (Refs. 1 and 2) and their implementation in COMSOL Multiphysics^®.

Origin of the Viscous Heat Equations

At a deeper level, the VHE originate from a systematic coarse-graining of the linearized phonon Boltzmann transport equation (LBTE), where the complex microscopic phonon dynamics are projected onto a small set of local-equilibrium fields: the temperature T(r,t) and a drift velocity u(r,t). This drift velocity emerges because, in the hydrodynamic regime, momentum-conserving (normal) phonon collisions dominate over momentum-relaxing (Umklapp) processes, allowing phonons to collectively carry crystal momentum. As a result, heat transport is no longer governed solely by temperature gradients but also by the evolution of this momentum field. The VHE therefore consist of two coupled equations:

• An energy balance equation for T(r,t)
• A momentum balance equation for u(r,t), analogous to the linearized Navier–Stokes equations

Crucially, two transport coefficients arise from the underlying LBTE symmetries: the thermal conductivity (associated with the odd-parity part of the phonon distribution) and the thermal viscosity (associated with the even-parity part). While the former describes diffusive heat flow, the latter accounts for momentum diffusion and viscous stresses in the phonon fluid, enabling phenomena such as heat vortices, backflow, and nonlocal transport. This framework naturally interpolates between regimes — it reduces to Fourier’s law when momentum is strongly dissipated and recovers previously proposed hydrodynamic models such as the Guyer–Krumhansl and dual-phase-lag equations as limiting cases, while remaining fully applicable to realistic materials with complex phonon dispersions.

From Microscopic Theory to a Continuum Model

At the microscopic level, phonon transport is described by the Boltzmann transport equation (BTE). While predictive, the BTE is extremely expensive to solve in complex geometries and is therefore ill-suited for device-scale modeling.

The viscous heat equations provide a mesoscopic alternative. They extend Fourier’s law by introducing an additional field: the phonon drift velocity. In this sense, they represent the phonon analogue of the linear Navier–Stokes equations for classical laminar fluids but applied to quantum phonon fluids. In this framework, temperature evolution is coupled to momentum balance equations for phonons, allowing the model to naturally interpolate between diffusive and hydrodynamic regimes.

A key feature of this approach is that all material parameters — such as thermal conductivity, viscosity, and relaxation rates — are determined from a single solution of the linearized BTE using first-principles calculations. These parameters are then incorporated into the continuum model, enabling efficient simulations in arbitrary geometries without compromising physical fidelity while retaining full quantum mechanical and ab initio accuracy.

Implementing the Viscous Heat Equations of Phonon Hydrodynamics

To make this framework accessible, we developed a custom physics interface using the Physics Builder in the COMSOL Multiphysics^® software. Once installed, the interface appears as Viscous Heat Equations (VHE).

Features Available After Implementing

• A coupled temperature–velocity formulation of heat transport
• Steady-state and time-dependent simulations
• Boundary conditions tailored to hydrodynamic heat flow
• Compatibility with complex geometries and mesoscopic devices

Installing the Interface

The installation follows the standard Physics Builder workflow:

• Enabling the Physics Builder in COMSOL preferences
• Importing the provided builder file
• Adding the Viscous Heat Equations (VHE) interface from My Physics Interfaces

Get step-by-step information on how to install this interface by expanding the section below.

Expand or collapse section…

Example 1: Fourier Diffusion vs. Hydrodynamic Heat Flow

To illustrate the qualitative differences between diffusive and viscous heat transport, consider a simple two-rectangle geometry. A vertical temperature gradient is applied across the main domain, while a smaller side region is laterally connected.

Under Fourier’s law, heat flows directly from hot to cold, producing smooth isotherms. When the viscous heat equations are used instead, a very different picture emerges: Heat recirculates in the smaller domain, forming vortex-like patterns and producing a small but finite thermal backflow. Below we show the comparison of steady-state temperature profiles obtained with Fourier’s law (left) and the viscous heat equations (right), showing hydrodynamic recirculation and thermal backflow.

Example 2: Thermal Vortices in the Incompressible Limit

The hydrodynamic nature of heat transport becomes even clearer in the incompressible flow limit of the VHE. In this regime, the divergence of the phonon drift velocity vanishes, and heat transport is dominated by momentum flow rather than diffusion.

Simulations show the spontaneous formation of steady-state thermal vortices (Ref. 3) — a behavior entirely absent in purely diffusion models.

Below we show the steady-state temperature profiles in the diffusive regime (top) and incompressible hydrodynamic limit (bottom), highlighting vortex formation and heat backflow (Ref. 3).

Time-Dependent Simulations: Hydrodynamic Heat in Motion

The interface also supports fully time-dependent simulations. By applying transient boundary conditions — such as oscillating drift velocities or alternating temperature gradients — it is possible to observe the formation, motion, and reversal of heat vortices over time.

These simulations reveal heat behaving much like a driven viscous fluid, responding dynamically to external forcing.

Why This Matters!

The interface provides a practical bridge between microscopic transport theory and continuum modeling. It enables the exploration of nondiffusive heat transport in realistic geometries that are inaccessible to fully microscopic methods.

By making hydrodynamic heat transport accessible within COMSOL Multiphysics^®, the interface opens new possibilities for studying thermal phenomena in next-generation materials and devices.

To learn more about the interface and download it, visit its Application Exchange entry: Phonon Hydrodynamics Interface. 

About the Guest Authors

Enrico Di Lucente is a postdoctoral research scientist in the Department of Applied Physics and Applied Mathematics at Columbia University in the City of New York, with a joint affiliation in materials science and engineering. He joined the research group of Prof. Michele Simoncelli in December 2025, immediately after completing his PhD in the Theory and Simulation of Materials (THEOS) group at EPFL, under the supervision of Prof. Nicola Marzari. He defended his doctoral thesis titled Theoretical and Computational Advances in Quantum and Hydrodynamic Thermal Transport.

His research focuses on thermal transport beyond Fourier’s law, phonon hydrodynamics, and first-principles modeling of heat transport and other fundamental and coupled excitations in condensed matter, ranging from magnons to light–matter interactions. His work bridges microscopic transport theory and continuum modeling, with the goal of enabling predictive simulations of nondiffusive heat transport in realistic materials and device geometries. It also aims to guide experimental efforts toward the design of innovative devices capable of detecting nonstandard quantum transport phenomena, with potential impact across emerging technologies including electronics, energy storage, fusion shielding, spintronics, and hypersonics. His broader research interests include condensed matter physics, computational physics, quantum transport, and magnetism.

Michele Simoncelli has been an assistant professor in the Department of Applied Physics and Applied Mathematics at Columbia University since January 2025. His group develops the theoretical and computational framework to understand, quantitatively describe, and control quantum transport phenomena in materials involving, e.g., charge, heat, light and spin, their possible synergies or conflicts, and related macroscopic signatures. Prior to joining Columbia, he held the Crone Research Fellowship in the Physics Department at the University of Cambridge (2021-2024). There, he worked on fundamental quantum theory and computational methods to describe the emergence of hybrid crystal–glass properties in materials with controlled degrees of atomistic disorder in, for example, chemical composition, bond network topology, or geometry. He received his PhD from EPFL (Switzerland) in 2021 under the supervision of Nicola Marzari, presenting in his thesis novel microscopic and mesoscopic theories of thermal transport in solids: the Wigner transport equation, generalizing the semiclassical Peierls–Boltzmann equation, and the viscous heat equations, generalizing Fourier’s law.

Acknowledgments

The development of the custom physics interface and the underlying theoretical work benefited from discussions and collaborations with Zhiyi Wang (Université Grenoble Alpes; THEOS, EPFL) and Prof. Nicola Marzari (THEOS, EPFL; Theory of Condensed Matter, Cavendish Laboratory, University of Cambridge).

References

1. M. Simoncelli, N. Marzari, and A. Cepellotti, Generalization of Fourier’s law into viscous heat equations, Physical Review X 10, 011019 (2020); https://journals.aps.org/prx/abstract/10.1103/PhysRevX.10.011019
2. Dragašević, B. Rajkov, and M. Simoncelli, Viscous heat backflow and temperature resonances in extreme thermal conductors. Physical Review Letters 136, 186302 (2026); https://doi.org/10.1103/nbbn-56hr
3. Di Lucente, F. Libbi, and N. Marzari, Vortices and backflow in hydrodynamic heat transport, Physical Review Letters 136(5), 056307 (2026); https://journals.aps.org/prl/abstract/10.1103/g9dx-hjyn

COMSOL Origins and Software Evolution

COMSOL AB was founded in Stockholm by Svante Littmarck and Farhad Saeidi. The company was originally the distributor of MATLAB^® software for all Nordic countries, but it grew into an independent entity in 1998 with its own software release.

“Svante and Farhad were grad students who believed mathematical modeling software was going to be something that engineers and scientists would need in the future,” said Bernt Nilsson, president of COMSOL, Inc. “They were right, and they built a sales organization in the Nordic countries and then released the company’s first multiphysics simulation software package.”

A CFD model from the first version of COMSOL Multiphysics^®.

COMSOL Multiphysics^® is a physics-based modeling and simulation platform used to predict how devices, processes, and systems behave under real-world conditions. It brings together multiple interacting physical phenomena in a single model and solves the underlying mathematical equations within a unified modeling environment. Advances in numerical methods, such as the finite element method, made it possible to apply these techniques to complex engineering problems involving realistic geometries, materials, and physical effects.

“These solutions did not happen by software alone,” said Svante Littmarck, CEO of COMSOL Group. “Moore’s law gave us rapidly increasing computer power, and at the same time, numerical algorithms improved a lot. Together, this meant that many previously out-of-reach problems could be solved. Once results such as stresses, temperatures, flows, electromagnetic fields, and chemical concentrations could be seen, simulation became not only a calculation tool but also a design tool.”

2000 to 2009: Independence and Growth

With independence came growth and change. COMSOL established its first additional offices in the United States and the United Kingdom in the late 90s and early 2000s, respectively. Additional locations were added across Europe through the 2000s, with offices in Germany, France, the Netherlands, and Switzerland.

This period introduced versions 1.2 through 3.5a. Our developers also created some of our users’ favorite modules, including the AC/DC Module, Acoustics Module, and Heat Transfer Module, as well as the Material Library.

This era also saw the release of the first edition of COMSOL News, our annual user-focused magazine that includes stories of how engineers and researchers use COMSOL Multiphysics^® in their daily work. The publication highlights how the software is being used in the development of batteries, loudspeakers, aerospace systems, medical devices, and many other areas.

Additionally, in 2005, COMSOL held a conference in Cambridge, Massachusetts, and in Las Vegas, Nevada. This conference offered a space for COMSOL staff and users to meet in person and discuss modeling and industry trends.

2010 to 2019: Software Improvements and Simulation Apps

The next era of COMSOL saw a focus on enhancing the software and extending the reach of simulation.

We released our LiveLink™ set of products, including LiveLink™ for MATLAB^® and LiveLink™ for AutoCAD^®, as well as numerous add-on modules, including the Battery Design Module, Nonlinear Structural Materials Module, and Semiconductor Module.

The software underwent significant change during this decade, with major additions that remain to this day. The Model Builder was introduced in version 4.0 and is the central workspace for building simulation models. In this workspace, users can define the geometry, materials, physics, mesh, studies, and results using the structured, easy-to-follow workflow. Like the rest of the COMSOL Multiphysics^® user interface, it has the same appearance and organization regardless of the engineering application or physical phenomena being modeled.

Introduced in version 4.2 was the Physics Builder, a graphical environment in which users can create physics interfaces through an interactive user interface without writing code. It is used extensively in the development of the COMSOL add-on modules.

Version 5.0 brought users the Application Builder, which provides a workspace for creating and maintaining custom simulation apps based on COMSOL models.

In 2014, version 5.0.1 introduced COMSOL Server™, which enables organizations to deploy, manage, and run simulation applications built by their in-house simulation experts. Users can then manage and distribute the apps, control who has access to them, and run the apps from a web browser on any computer.

The Model Builder when it was first released in 4.0 (left) and an example simulation app from 5.0 (right).

In 2018, COMSOL Compiler™ was also added to the COMSOL product suite. This product enables users to compile their simulation apps into standalone executable files that can be run with or without internet access and with or without COMSOL licenses. Users can distribute these files for free or charge a fee.

During this same time period, we continued to grow by adding offices in China and India, and in 2013, we held the COMSOL Conference in Boston, Massachusetts; Bengaluru, India; and Rotterdam, the Netherlands.

2020 to 2026: Simulation Democratization and Acceleration

The 2020s have seen a greater democratization of simulation and an emphasis on reducing solving times for complex models. With version 6.0, we introduced the Model Manager, which enables colleagues to collaborate and centrally organize their custom models, apps, reports, presentations, and other project assets within simulation projects. Version 6.2 introduced a new surrogate model framework that enables users to create compact models that reconstruct simulation results in a fraction of a second. These surrogate models enable fast simulation apps, accelerate optimization and uncertainty quantification, and support the development of effective digital twins. This release also improved the speed of solving CFD models by up to 40%.

In version 6.3, GPU acceleration was added, offering up to 25x faster acoustics simulations and surrogate model training. This version also introduced the Java Shell window, an interactive environment for developing and modifying models using the COMSOL API. An optional AI-powered chatbot window can assist users with Java programming, model development, and many other aspects of using COMSOL Multiphysics^®.

These updates bring us to the latest and most powerful version of our software: 6.4. This version features solver performance that is greatly enhanced through the NVIDIA CUDA^® direct sparse solver (cuDSS) for NVIDIA GPUs, providing several-fold speedups for both single-physics and multiphysics simulations.

A 40× speedup in a wave-based room acoustics simulation enabled by NVIDIA cuDSS.

With this version, we also released the Granular Flow Module, which enables the simulation of granular materials by modeling the motion and interaction of individual particles in bulk solids processes.

A ribbon mixer model created with the Granular Flow Module, released with version 6.4.

This year, we are continuing to focus on highlighting our users’ innovations. We recently published the 2026 edition of COMSOL News, which features, among other organizations, Volvo Trucks and Resolvent, UK Fusion Energy, and the Technical University of Denmark. Additionally, the COMSOL Conference continues to be a staple for us and the simulation community, with this year’s tour including stops in Cambridge, UK; Boston; Tokyo; Shenzhen; Bengaluru; Málaga; and Taipei. The conference highlights the work of our users through posters, papers, and live presentations, and attendees are able to hear from industry leaders and learn more about modeling. (If you want to learn more about the 2026 COMSOL Conference tour, click here.)

Moreover, in 2026, we hosted our inaugural Simulation Summit in Santa Clara, California. This one-day event brought together engineers and researchers to attend technical sessions, keynote sessions on how leading organizations are using multiphysics simulation, and a panel discussion where experts discussed the future of simulation tools.

The People Who Drive COMSOL

COMSOL has grown to 16 offices worldwide. “Our employees built everything,” said Littmarck. “We have more than 650 people, and we need every single one. Everybody is contributing.”

Throughout our many years and the growth of our company, our vision has remained the same: to empower users with the tools they need to improve their designs and processes and bring their ideas to life.

In the age of AI and advanced computing, we are excited about what the future has in store, but we also know that whatever comes next will be shaped by the people who make it all possible: our users and employees. “It’s been an amazing 40 years,” said Nilsson, “but when I think about it, the best is yet to come.”

The Bell that Inspired the Simulation

While in Philadelphia recently for an Acoustical Society of America meeting, I saw the Liberty Bell display and was reminded of some of its history. At the conference, I was also pleasantly surprised to learn about past (Ref. 1) and ongoing (Ref. 2–3) research regarding structural acoustics of the bell.

The Liberty Bell. Image by William Zhan — Own work. Licensed under CC BY 2.0, via Flickr Creative Commons.

The Liberty Bell was ordered for the Pennsylvania State House and cast by the Whitechapel Bell Foundry in London (Ref. 4). The original imported bell cracked during its first test ring and was recast locally by John Pass and John Stow. The famous crack visible today is a later feature of the recast bell; according to the National Park Service (Ref. 4), its exact origin is not recorded, but a narrow split likely developed in the early 1840s and was widened in 1846 by metal workers to prevent it from spreading further. Today, the Liberty Bell is treasured as a national icon. It is currently on display at the Liberty Bell Center at Independence National Historical Park in Philadelphia, Pennsylvania.

Modeling the Liberty Bell

Our first goal as simulation engineers interested in structural acoustics is to compute the structural modes of the bell. To do this, we will need some baseline material properties, geometry, and physics with boundary conditions. The material properties for so-called bell bronze can vary dramatically, but for now we will assume the elastic modulus, density, and Poisson ratio are constant () and employ a linear elastic material model. Since the bells are designed to ring, the structural damping is quite low, so it will be neglected for computing the modes.

The geometry of the bell is very important. In fact, the acoustics of church bells have been analyzed extensively, and bell foundries will typically tune the first five modes with particular ratios in order to achieve optimal tuning. The names of the first five modes as well as the frequency ratios for a well-tuned bell (Ref. 5) are: hum ( ); fundamental, also known as prime, (); tierce (); quint (); and nominal (). It is not known if or how precisely the Liberty Bell was tuned, which is one of the questions the first simulation will investigate.

In order to obtain a realistic geometry, we will utilize a 3D scan of a Liberty Bell replica (see caption). The scan is a surface mesh that contains all the inscriptions, mounting structure, clapper, and some defects from the scan. However, using the geometry tools in COMSOL^®, we extract a representative cross section and revolve this to create the solid geometry. We will also ignore some of the features, since it would add complexity and should not significantly contribute to the ringing modes of the bell.

At left: Geometry of the Liberty Bell replica from a 3D scan by the VAR Lab at Penn State Behrend, licensed under CC BY 4.0 International. At right: The geometry was modified and defeatured for simulation.

For the physics, we use the Solid Mechanics interface and leave the default Free boundary condition on all boundaries. An eigenfrequency study is used to compute the first five modes of the structure only. The model predicts a fundamental of . Compared to a perfectly tuned bell, the replica with baseline properties is slightly off tune (see the harmonic ratios below). We should also note that each of these modes has an orthogonal pair (degenerate at the same frequency), which will be important to remember for later analysis.

First five eigenfrequencies and their corresponding mode shapes.

Uncertainty Quantification

So far, we have computed the modes assuming some baseline properties. The elastic and mass properties of bell bronze, however, are not well documented. Furthermore, its composition could vary depending on the specific alloy mix of the casting. A natural question to ask is: What effect does the uncertainty in the material properties have on the eigenfrequencies? This question can be answered by using the Uncertainty Quantification Module, an add-on to COMSOL Multiphysics^®.

Instead of defining discrete material properties, we will define distribution ranges for them. In particular, we set the properties to be uniformly distributed throughout the following ranges: . We will use the Uncertainty Propagation study type along with an eigenfrequency reference study to determine a 95% prediction interval for the modes of interest. The method utilizes a surrogate model in order to efficiently compute the uncertainties.

The result from the analysis shows us error bars with the expected range of eigenfrequencies. From these results, we can note that the variation increases with increasing mode number (going from +/- 17 Hz for the hum to +/- 60 Hz for the nominal). Furthermore, experimental measurement from Ref. 1 has modes that do fall within the range of the error bars, giving us confidence that we have reasonably bracketed the properties of the bell.

The mean eigenfrequencies (blue) and 95% confidence interval (black error bars) for the eigenfrequencies of the bell, given uniform distributions of the mass and elastic material properties.

What About the Crack?

A cracked bell can still vibrate and radiate sound, but the crack can strongly alter its tonal character. To investigate, we add a 1 mm crack to the geometry, qualitatively similar to the location of the large crack. The crack locally reduces stiffness and breaks the axisymmetry of the structure, leading to mode splitting. The formerly degenerate mode pairs split, producing 10 distinct eigenfrequencies over the frequency range considered (shown below at right). This, of course, results in modes that are vastly different from the design guidelines and lack harmony.

Cracked modes of the structure: cracked geometry and nominal (left) and all computed modes of the cracked bell and baseline bell (right).

Tuning the Bell

Can simulation be used to create a nearly perfect bell? We can set up a shape optimization problem to help find out. The boundaries of the bell are allowed to freely deform in an attempt to minimize the objective function — the squared difference between the current eigenfrequencies and the target ones. The shape optimization runs for 20 iterations, and the result is a new design that achieves the target frequencies within less than 0.02% error for all frequencies! Essentially, it is able to almost perfectly match the guidelines in Ref. 5. Interestingly enough, the new shape closely resembles the baseline shape, which truly provides appreciation for the high precision involved in casting and machining bells for high acoustic quality.

Result of the shape optimization: one half of an axisymmetric view of the bell cross section showing baseline geometry (red) and shape-optimized geometry (green).

Acoustic–Structure Interaction

Lastly, we model the sound radiation response from the defeatured Liberty Bell replica by using the Acoustic–Solid Interaction, Frequency Domain interface. We also assume a loss factor of . A boundary load condition is applied on a representative area of where the clapper would strike. The model computes the full field displacements, stresses, and strains, as well as acoustic pressure and sound pressure level (SPL) resulting from the force.

SPL vs. frequency for observer location adjacent to the bell at (0,0,-2m).

The results can be animated over an acoustic cycle to get the full picture of the vibroacoustics.

Further Resources

In this blog post, we have presented several structural acoustic simulations of the Liberty Bell. However, the possibilities do not have to end there. Simulation can be used to investigate many other aspects of the bell, including: sound propagation from the Liberty Bell and how far away the bell could be heard, dynamic crack propagation caused by repeated impacts, and molten metal solidification during the bell casting process.

Additional information on the technical topics mentioned in this blog post can be found below:

• Performing Optimization in COMSOL Multiphysics
• Introduction to Uncertainty Quantification
• Getting Started with Modeling Structural Mechanics
• Introduction to Modeling Acoustic-Structure Interactions

References

1. G.H. Koopmann et al., “Tuning a Replica of the Liberty Bell via Material Tailoring: An Application of a Method for Optimal Acoustic Design.” ASME International Mechanical Engineering Congress and Exposition, vol. 35517, American Society of Mechanical Engineers, 2001.
2. J. Young, S. Collier, A. Stearns, and E. Brown, “Modes of America: Computational Acoustics and the Sound of the Liberty Bell,” Presented at the 190th Meeting of the Acoustical Society of America (ASA), Philadelphia, PA, 2026.
3. S. Collier, J. Young, A. Stearns, and E. Brown, “Modal Analysis of a Liberty Bell Replica,” Presented at the 190th Meeting of the Acoustical Society of America (ASA), Philadelphia, PA, 2026.
4. “The Liberty Bell,” June 2025; https://www.nps.gov/inde/learn/historyculture/stories-libertybell.htm.
5. Thomas D. Rossing and Robert Perrin, “Vibrations of Bells,” Applied Acoustics, vol. 20.1, pp. 41–70, 1987.

Key Themes from OFC 2026

OFC continues to be the premier event where the optical ecosystem converges, from component manufacturers and system designers to network operators and researchers. With hundreds of companies presenting the latest solutions, the event provides a unique opportunity to see how research and development is being translated into real-world applications.

Over 17,000 industry members from 91 countries attended this year’s conference, and participants attended panels and tutorials led by representatives from companies like Nokia, NVIDIA, and Google.

The venue for OFC 2026.

Across the exhibition floor, the main focus was on AI and enabling technologies that increase bandwidth, decrease latency, and improve energy efficiency. This was reflected in the wide range of technologies on display, from hollow-core and multicore fibers to optical transceivers and modulators and even thermal management solutions. These innovations highlight how the industry is evolving to support the AI boom as it both scales up and scales out.

Insights from the Show Floor

Beyond showcasing the technologies themselves, the exhibition floor offered a valuable opportunity for attendees to connect directly with engineers, researchers, and industry professionals. Conversations throughout the event highlighted how rapidly the optics space is evolving, particularly as new demands are placed on network infrastructure. For COMSOL, exhibiting at OFC for the first time provided us with the chance to engage with attendees working across a wide range of applications.

Andrew Strikwerda, Tyler Tippens, and Morgan Castadoro at COMSOL’s exhibition booth.

Many discussions centered around the challenges of integration and packaging for these complex devices, especially as systems become more compact and performance expectations continue to rise. The variety of companies present also underscored how interconnected the industry has become. From component-level innovations like hollow-core fibers to enabling technologies like co-packaged optics, it was clear that advancements across the entire ecosystem are contributing to the next generation of communication technologies.

The Role of Simulation in Optical Design

As the technologies showcased at OFC continue to advance, the need for accurate and efficient design tools becomes increasingly important.

Engineers are often required to consider multiple physical effects simultaneously, particularly when working with high-performance optical and electronic systems. Simulation addresses these challenges by allowing teams to model and evaluate designs before moving to physical prototypes. The multiphysics modeling capabilities of the COMSOL^® software are especially valuable when working with components like photonic chips, connectors, and high-speed interconnects, where performance can be influenced by a combination of electromagnetic, thermal, and structural factors.

Andrew Strikwerda, senior application manager at COMSOL and attendee of OFC, further emphasized this point, saying, “With the growth of AI and the ever-increasing need for high bandwidth and low-latency data transfer, the field of fiber optics is growing rapidly. Simulation is more important than ever, not only for designing next-generation fibers like the hollow-core antiresonant fiber but also for complex optical packaging problems where thermal, mechanical, and optical effects interact.”

By enabling a deeper understanding of the interplay of these effects, multiphysics modeling software helps reduce development time while improving overall design accuracy. As systems continue to scale and evolve, this approach becomes even more critical.

Simulation in the Future of the Optics Industry

OFC 2026 highlighted the growing importance of optical technologies in supporting the demands of an increasingly connected and data-driven world. From advancements in hardware to broader system-level innovation, the event showcased how the industry is adapting to meet the needs of modern communication networks. For COMSOL, exhibiting at OFC marked an important step in connecting with optical communications professionals to best support the development of next-generation technologies.

Modeling Wave Optics

If you’re interested in learning more about how to use the COMSOL^® software to model optical technologies, take a look at the Wave Optics Module add-on product below.

Ahead of Mexico’s home-turf match against Czechia tomorrow at this legendary stadium, let’s look at its acoustics in detail.

Sacred Football Ground

Together with Rio de Janeiro’s Maracanã Stadium, Estadio Azteca (dubbed “Mexico City Stadium” by FIFA^® for the 2026 tournament) is one of only two stadiums to have hosted two World Cup finals. But the Maracanã of the 2014 World Cup final is essentially a different stadium from the one where Uruguay defeated Brazil in front of nearly 200,000 fans in the Maracanazo of 1950. The Azteca that hosted Pelé in 1970, on the other hand, is still recognizably the same stadium that’s hosting the World Cup today.

No stadium has witnessed more historic World Cup moments. Pelé won his third World Cup here. In 1986, Maradona scored both the “Hand of God” and the “Goal of the Century” here, arguably the most controversial goal and the greatest goal in football history, separated by only four minutes.

When it hosted the 2026 opening match between Mexico and South Africa, Estadio Azteca became the first stadium in history to host matches in three FIFA World Cups. Reverently referred to by fans as el gigante (“the giant”, as it was immortalized in song by Andrés Calamaro) and El Coloso de Santa Úrsula for its massive capacity, the Azteca is sacred ground for football.

Figure 1. A COMSOL Multiphysics^® representation of Estadio Azteca. Note that this is not an exact representation of the stadium, but it’s good enough for our investigation, which is mainly for fun.

The Acoustics of a Football Stadium

The sound of a football stadium is about much more than acoustics. Every supporter knows that their home stadium sounds better than any other stadium in the world when their team scores. This is one of the great experiences of football fandom.

As engineers, however, we need to quantify the sound of a stadium. FIFA specifies requirements for quantities such as the speech transmission index, reverberation time, and the uniformity of the sound field in the stands (Ref. 1). These quantities can be measured objectively and are used when designing the stadium sound system to ensure that spectators can clearly hear announcements and other information provided through the PA system.

FIFA specifies the following design targets for these quantities:

• Reverberation time should be no more than 4 s in the frequency range 125–4000 Hz.
• Speech transmission index should exceed 0.55 for a full stadium. (The recommendation is 0.75.)
• Nonuniformity of the sound field should be no more than ±3 dB.

These quantities can also be simulated, which is what our team did to estimate the quality of the sound from Estadio Azteca’s brand new PA system. We started by focusing on the sound from a single loudspeaker cluster, simulating how it propagates through the stadium.

Can You Hear the Giant?

The new sound system of Estadio Azteca was installed in 2026 and appears to comprise roughly 340 loudspeakers. Based on photographs from the renovation and information published on social media, the system appears to be supplied by d&b audiotechnik^®. The loudspeaker clusters hanging from the roof appear to consist of four loudspeaker cabinets and two subwoofers. Based on their appearance, I believe the loudspeakers may belong to the d&b audiotechnik Vi or Yi series.

For our simulations, we took one of these loudspeaker clusters and placed it just below the roof, as shown in Figure 1. Figure 2 shows a close-up of the loudspeaker cluster and the resulting total acoustic pressure in the stands below.

Figure 2. The loudspeaker cluster hanging from the roof of Estadio Azteca and the resulting total acoustic pressure in the stands below.

We created the geometry using the built-in geometry tools in COMSOL Multiphysics^® and used the Pressure Acoustics, Time Explicit interface to model the acoustics in the low frequency range and the Ray Acoustics interface for the high frequency range.

The Pressure Acoustics, Time Explicit interface automatically defined the numerical model using fourth-order discontinuous-Galerkin-based functions, which for a representative frequency of 100 Hz required the mesh shown in Figure 3. The mesh consisted of grid elements (a feature that will be available in the upcoming release) in the bulk, with pyramids and tetrahedrons close to the stands, the pitch, and the roof. The resulting system of equations contained 99 million degrees of freedom and was solved in 1 hr 55 min on two NVIDIA RTX^™ 6000 Ada Generation GPUs.

Figure 3. A cross section of the volume mesh and the surface mesh used in the numerical model of Estadio Azteca.

Figure 4 shows an animation of the total acoustic pressure (visualized on surfaces) during 0.4 s in the stands that are reached by the output of the central loudspeaker cluster. Based on available photographs of Estadio Azteca, there appear to be roughly 30–40 such clusters hanging from the roof and additional clusters located in the stands and along the sides of the pitch. If we placed all of these loudspeaker clusters around the stadium in the model, we could estimate diffraction effects and the combined contribution from all sound sources at low frequencies throughout the stands.

Figure 4. Animation of the total acoustic pressure generated by the central loudspeaker cluster hanging from the roof.

In the low frequency range, effects such as diffraction are important and are best captured using a wave-based method, as illustrated above. As the analyzed frequency increases, the computational cost also increases and it is common practice to switch to a high-frequency method such as ray tracing.

In the Ray Acoustics interface, we can easily define a source with a given spatial directivity. The source data can be imported from a file or created in a model, as in this case, where we used the model of the loudspeakers and subwoofers. A simplified model representation of the loudspeaker cluster at the Azteca is shown in Figure 5. The sound radiation pattern is computed using the Pressure Acoustics, Boundary Elements interface (based on the boundary element method, or BEM).

Figure 5. Radiation plot (bubble plot) of the simplified loudspeaker cluster used at Estadio Azteca at 1000 Hz, computed using the BEM. Note that the level shown is relative.

The propagation of rays from the speaker cluster, with the source characteristics shown in Figure 5, is illustrated in the animation in Figure 6. The animation includes just 10,000 rays, but a practical simulation could easily be performed using many more. The model equations solved in minutes.

Figure 6: Ray propagation from the speaker array.

The resulting sound pressure level map in a plane just above the seating area is shown in Figures 7 and 8, for one and two speakers, respectively. You can see in both plots that the nonuniformity is much larger than FIFA’s target of ±3 dB. However, adding just one additional speaker improves the coverage. In a real stadium design, additional loudspeaker clusters would likely be used to improve the coverage of the central part of the stands. From fan footage of Estadio Azteca, you can tell that at least 4–5 of the loudspeaker clusters hanging from the roof would likely cover the plot area in Figures 7 and 8.

Figure 7. Sound pressure level map just above parts of the seating area below the location of one speaker.

Figure 8. Sound pressure level map just above parts of the seating area below the locations of two speakers.

The impulse response, reverberation time, and speech transmission index could also be computed using this model. And to estimate the full acoustic coverage in Estadio Azteca, we could place all of the loudspeaker clusters throughout the model.

The Sound of the Beautiful Game

Several different types of models are needed to design the sound system of a stadium. In our case, we used a pressure acoustics model in the time domain for diffraction effects and the low frequencies, a BEM model to characterize the sound radiation from the loudspeaker cluster at higher frequencies, and a ray acoustics model to estimate the resulting sound pressure levels at higher frequencies. Together, these models can provide an accurate picture of the quality of the sound from the PA system.

One aspect of stadium acoustics that we didn’t consider? The roar of the crowd. That is the real sound of the giant.

More than 80,000 people experienced this thunderous wall of sound during the opening match between Mexico and South Africa. Some of them may have felt the same rush that Andrés Calamaro sings about in “Estadio Azteca”, where he describes the stadium as a giant that has “crushed” him. This feeling is not fear but awe. Anyone who grew up loving football recognizes that feeling when entering a huge stadium for the first time. In my case, it was walking into Estadio Centenario in Montevideo holding my father’s hand.

The sound of the crowd may very well be the subject of my next blog post.

For the Love of the Game (Only!)

Although the simulations presented here are based on established acoustics modeling techniques, they were created primarily for fun. A professional acoustics study of Estadio Azteca would require significantly more detailed information about the stadium geometry, loudspeaker system, materials, crowd distribution, and operating conditions than is publicly available.

The loudspeaker system used in the simulations was reconstructed from publicly available photographs and information from media reports and social media. We therefore make no claim that the model accurately represents the actual sound system installed in Estadio Azteca for the 2026 FIFA World Cup.

These investigations were performed independently of FIFA, Estadio Azteca, and d&b audiotechnik, and we do not claim any cooperation with these organizations.

Reference

1. A. Peretokin et al., “Acoustics Features of Sports Facilities on the Example of FIFA 2018 Football Stadiums in Russia,” Proc. 23rd Int’l Cong. Acoust., Integ. 4th EAA Euroregio (ICA 2019), pp. 811–818, 2019.

Adidas and Trionda are registered trademarks of adidas AG.

COMSOL AB and its subsidiaries and products are not affiliated with, endorsed, by, sponsored by, or supported by any of the foregoing trademark owners.

Defining IUQ

In our previous blog post “How Reliable Is Your Resistor?”, we discussed how forward UQ predicts how variations in material properties, geometry, or manufacturing processes affect a resistor’s performance. In many cases, however, measurement data is already available (e.g., resistance values from experiments), and it becomes more interesting to calibrate the input parameters, especially the material properties. In the COMSOL^® software, IUQ can be used for this task, as it works backward from observed data to estimate those unknown inputs, combining the finite element method with surrogate models to efficiently calibrate the input parameters.

While parameter estimation focuses on identifying the best fit parameter values, IUQ provides a probabilistic description of the calibration parameters, including their likely ranges, confidence levels, and correlations. By refining parameter distributions based on measurements, IUQ improves model accuracy and predictive capability.

Moreover, it’s often desirable to perform an IUQ study when experimental data is available and input distribution is unknown. For example, after identifying key parameters using screening or a sensitivity analysis, an IUQ study can be used to obtain posterior distributions for these parameters. These calibrated distributions can then be used as input parameters for a forward UQ study like uncertainty propagation and reliability analysis. This procedure creates a more realistic workflow since the uncertainties used for prediction are informed by data and the prior assumption.

Here, we will discuss how to perform an IUQ study and how to seamlessly use the resulting posterior distributions in a forward uncertainty propagation study in the COMSOL^® software.

Understanding IUQ with a Resistor Model

IUQ estimates calibration parameters by combining experimental data with prior knowledge. It can be viewed as parameter estimation in a Bayesian framework, where measurements guide the updating of parameter distributions.

The experimental data provides the reference quantities that the model must reproduce, such as the measured resistance at different applied voltages in a resistor. During an IUQ study, COMSOL^® evaluates how likely it is for each parameter set to occur by comparing simulated outputs with these measurements.

In the software, IUQ compares the experimental data with predictions from a surrogate model trained using finite element simulations. The finite element model defines the physics-based relationship between uncertain inputs and measurable outputs, whereas the surrogate model approximates this relationship and enables efficient sampling of the parameter space.

The Bayesian updating process combines prior distributions with a likelihood function that measures the consistency between surrogate predictions and experimental data. The resulting posterior distributions represent the most likely values of the calibration parameters together with their uncertainty. Through this process, IUQ provides calibrated material properties in a data-informed and computationally efficient way. Thus, experimental data directly influences the inferred parameter distributions.

Workflow: Performing IUQ in COMSOL^®

Define the Problem and Build the FEM Model and IUQ Study

The first step of setting up an IUQ study is to create a physics-based model representing the system at hand. In the resistor example highlighted in our previous blog post, the Electric Currents interface is used to compute the resistance based on specified material properties and geometry. Two conductivities, Sigma1 and Sigma2, defined on different regions of the resistor, are selected as the calibration parameters. A stationary study establishes the forward relationship between the input parameters and the output quantity, which is the resistance. For our IUQ example, we will keep these settings.

The next step is to add an Uncertainty Quantification study, using Study 1 as the reference, and selecting Inverse uncertainty quantification as the study type. Both Gaussian process and polynomial chaos expansion approaches can be used as surrogate models for IUQ. In this example, an adaptive sparse polynomial chaos expansion is used. The variable comp1.Res, which evaluates the resistance of the resistor, is defined as the quantity of interest (QoI) and refers to the simulation output.

Figure 1. The Uncertainty Quantification study settings, with Inverse uncertainty quantification selected.

Define Uncertainty Input Parameters and Their Prior Distributions

Under Input Parameters, include both the calibration parameters (Sigma1 and Sigma2) and the experimental parameter, which is the applied voltage, V0. Each value of V0 corresponds to a different experiment. Therefore, V0 is treated as an experimental parameter rather than a calibration parameter.

Next, provide a prior distribution for each parameter. In this demonstration, normal distributions are assumed for Sigma1 and Sigma2. For the experimental parameter V0, a uniform distribution is used because the measurements are performed over a range of voltages.

Figure 2. Input parameters and prior distributions.

The experimental parameter V0 must be included because the surrogate model represents the system response as a function of both the calibration parameters and the experimental condition. The bounds of the V0 distribution (e.g., from 4.9 to 35.1 V) should cover the voltage range used in the experiments, for example, from 5 to 35 V.

Prepare and Import Experimental Data

In the Experimental Data Settings section, experimental values can be provided in table format. External measurement data can be imported from files such as .txt or .csv. For this demonstration, pseudoexperimental resistance data is generated by running the Stationary study with an auxiliary sweep over the voltage V0. A measured resistance variable is defined as

where Res is the resistance evaluated from the finite element method (FEM) model and rn1() is a random function. The random function rn1() adds measurement noise to the simulated resistance. The evaluated values of Res_measured as a function of V0 are then selected in the Experimental data table in the IUQ settings.

Figure 3. The experimental data settings with Calibrated selected as the measurement uncertainty type.

Here, V0 (the voltage parameter) is selected as the experiment parameter, and the pseudoexperimental data, Res_measured evaluated under Global Evaluation, is defined as the QoI.

If Measurement uncertainty type is set to Calibrated, COMSOL^® estimates the measurement uncertainty as part of the IUQ process. If Experimental data is selected instead, the measurement uncertainty must be provided explicitly.

Samplings and Surrogate Model Settings

Under the sampling settings, the maximum number of input points can be specified, which directly controls the size of the training dataset generated from FEM simulations. These input points are used to train the surrogate model based on FEM simulations.

Note: In this demo model, before running the IUQ study, the Auxiliary sweep in the Stationary step should be disabled.

IUQ Results

The IUQ study produces the joint probability distribution and the calibrated confidence intervals shown in Figures 4 and 5, respectively. Figure 6 compares the prior and posterior distributions for Sigma1 and Sigma2, demonstrating how the experimental data refines the parameter estimates and reduces uncertainty, especially for Sigma2. Note that Sigma1 is the conductivity of the resistive material, and the resistance is insensitive to Sigma1, as demonstrated in the previous screening analysis. Thus, the prior and posterior distributions for Sigma1 are quite similar.

Figure 4. Joint probability distribution for Markov chain Monte Carlo (MCMC) samples.

Figure 5. The calibrated confidence interval.

Figure 6. A comparison between the prior and posterior distributions for Sigma1 and Sigma2.

Using Posterior Distributions for a Forward UQ Study

The posterior distributions obtained from an IUQ study can be directly reused in a forward UQ study. To implement the distributions, the first step is to add a new Uncertainty Quantification study of type Uncertainty Propagation from the IUQ study.

Figure 7. Adding a new UQ study for uncertainty propagation.

A new Quantities of Interest table is created automatically and is identical to the one used for the IUQ study. Since the same QoI, comp1.Res, is evaluated, either Analyze only or Improve and analyze can be selected. The former reuses the existing FEM data to train the new surrogate model, e.g., a new Gaussian Process model, while the latter allows additional simulations or input points to improve the surrogate model. For simplicity, we will use Analyze only for the forward UQ studies.

By default, a new Gaussian Process surrogate model is generated, which is used as the surrogate model for this new forward UQ study. From there, enable the new Gaussian Process and click Train Model. Through these steps, we can add the posterior distributions from the MCMC samples for the forward UQ study.

Figure 8. The settings for the Gaussian Process surrogate model.

Note that the distributions defined under the Input parameters section can be ignored for this forward UQ study. Instead, the posterior distributions obtained from the IUQ study fully define the sampling space. Therefore, no additional prior assumptions are required for the forward UQ sampling.

In the Surrogate-Based Monte Carlo Analysis section, select Manual for the Monte Carlo parameters source. Then, for each calibrated parameter, Sigma1 and Sigma2, select Data as the source type, Result table as the data source, and the corresponding column from the MCMC samples table for each parameter.

Note that we can select only one working condition for V0, which is the nominal value of 1 V.

Figure 9. The settings in the Surrogate-Based Monte Carlo Analysis section for using the posterior distributions.

After this forward UQ study is finished, the kernel density estimation (KDE) plot based on the posterior distributions can be generated. This KDE reflects the uncertainty propagated from the calibrated parameter distributions. Figure 10 shows that the highest probability density of the resistance is located around 50 . According to the QoI confidence interval table (not shown), the mean predicted resistance is 49.975 , with a standard deviation of 0.0957 .

Similarly, a reliability analysis can be performed with the posterior distributions specified in the Surrogate-Based Monte Carlo Analysis section (Figure 9). For example, the probability of the resistance that is larger than 50.25 and 52.5 is approximately 2.55 and 0%, respectively. This probability quantifies the risk of exceeding specified resistance thresholds. Another way to say this is that we know that this resistor is more reliable for any applications requiring a precise 50 load.

Figure 10. The KDE plot using posterior distributions.

A Data-Informed Workflow for Calibration and Uncertainty Propagation

Inverse uncertainty quantification provides a systematic way to calibrate uncertain input parameters using experimental data. By combining finite element modeling, surrogate modeling, and Bayesian updating, IUQ refines prior assumptions into posterior distributions that better represent the physical system.

When these posterior distributions are subsequently used in a forward UQ study, the predictions become more realistic and data informed. Instead of relying on assumed parameter variations, the uncertainty propagation is based on calibrated parameter distributions that reflect both measurements and model physics.

This combined workflow of IUQ followed by forward UQ enables more reliable prediction, improved confidence in simulation results, and a clearer understanding of how parameter uncertainty influences system performance. This workflow integrates model calibration and uncertainty propagation and can be completed within the COMSOL^® software’s dedicated user interface, without needing to switch tools.

Next Steps

Want to learn more and try out the model discussed above? Download the related MPH file below:

The Sensors on the Trionda

On May 3 this year, Manchester United secured a dramatic 3–2 win against Liverpool at Old Trafford. In the 14^th minute, Benjamin Šeško scored the 2–0 goal for Manchester United. The goal became controversial because Šeško may have touched the ball with his hand just before it crossed the goal line.

However, because the Premier League does not use balls with embedded inertial measurement unit (IMU) chips, the VAR officials spent several minutes reviewing broadcast camera footage before ruling the images inconclusive and allowing the goal to stand.

This type of situation is less likely to occur at the 2026 World Cup. The Trionda contains a 500-Hz IMU chip consisting of MEMS accelerometers and gyroscopes. Even relatively small changes in acceleration and angular velocity caused by a touch of the ball can therefore be detected and analyzed in real time.

The sensors in the Trionda have already played a decisive role in the World Cup this year. In Sunday’s game between Sweden and Tunisia, Mattias Svanberg’s 4–1 goal for Sweden was initially disallowed because he was called offside. It appeared that Svanberg had received the ball directly from the free kick, in which case he would have indeed been offside.

However, the sensors inside the ball detected a tip-of-the-toe touch by Svanberg’s teammate Alexander Isak between the free kick and the shot on goal. Whether Isak touched the ball was difficult, if not impossible, to confirm from the camera footage alone. Because the sensors could determine the exact moment of contact, VAR was able to establish that the pass leading to Svanberg’s shot came from Isak rather than directly from the free kick. At the moment Isak touched the ball, Svanberg was no longer in an offside position, and the goal was therefore allowed to stand.

To add to the excitement, Svanberg had come off the bench mere seconds before the play and the goal was his very first touch of the match. Talk about timing, by both the coach and the player!

Schematic illustration created in COMSOL Multiphysics^® of the Adidas Trionda, with the sensor chip placed in one of the ball’s four surface panels. Counterweights in the other three panels make the ball mechanically symmetric.

The chip inside the Trionda that’s making such call reversals possible was developed by Adidas in cooperation with KINEXON Sports. They have not publicly disclosed which IMU product was used to develop the chip. However, the published specifications of commercially available devices such as the TDK^® InvenSense^® ICM-20649 and ICM-45686 are strikingly similar to what one would expect from an IMU designed for football applications. These devices support angular velocities up to ±4000 degrees per second (dps) and accelerations up to ±32 g. Interestingly, the datasheet of the ICM-20649 explicitly mentions “soccer ball kicks” as a target application. TDK’s MEMS MotionTracking^® devices integrate three MEMS accelerometers and three MEMS gyroscopes, together with signal conditioning electronics, analog-to-digital converters, temperature sensors, and communication interfaces, into a compact 2.5 × 3 × 0.81-mm hermetically sealed package.

A MEMS gyroscope can be modeled as a vibratory rate gyroscope, where Coriolis forces couple a driven vibration mode to a sensing mode. Although the geometry of the IMU used in the Trionda is proprietary and thus its gyroscopes can’t be recreated exactly for simulation, COMSOL Multiphysics^® models of comb-drive tuning fork gyroscopes illustrate this principle.

In addition to the IMU, the electronics package inside the Trionda also contains a local positioning transmitter that sends timing and positioning data to stadium anchors, an ultra-wideband (UWB) RF antenna for communication with the anchor infrastructure, and a battery. Like the IMU, the RF electronics and batteries can be analyzed using COMSOL Multiphysics^®.

The IMU and positioning system can be used by the VAR team for more than just determining whether the ball was touched in a potential handball situation, like the Šeško goal I mentioned earlier. It can also be used to determine the exact moment when the ball is touched during a pass, making it possible to determine whether the receiving player was offside at the exact moment of the pass. The high sampling frequency allows the system to distinguish between successive touches during rebounds and deflections, which can be important in crowded situations inside the penalty area. Additionally, the positioning system can determine where the ball is relative to the goal line (presumably in combination with image processing).

The gyroscopes can also be used to measure the angular velocity and spin axis of the ball in real time during flight, making it possible to analyze curl, Magnus effect trajectories, and low-spin “knuckleball” shots in far greater detail than before.

It would be even nicer if we, in front of our TV screens, could also see statistics for the hardest shots and the shots with the most curl during a game!

The Sensors on the Vests

The vests that the players wear under their jerseys contain even more sensors and electronics than the ball. While the IMU inside the ball is designed to measure and withstand hard shots and rapidly spinning balls, the IMUs in the players’ vests are optimized for long-duration motion tracking and sensing during running and changes in direction.

If you look closely at the players during a game, you may notice a small protrusion under their shirts between their shoulder blades. That little hump is the pod of the vest, where most of the sensors and electronics are located. The ECG electrodes that measure heart rate and heart rate variability are the only sensing devices placed outside the pod. They are positioned in the front lower part of the vest.

You can often see the vest after games, when the players exchange shirts.

Schematic image of the sensors in the vest and in the pod. In the real vest, the pod is hidden inside a pocket, and the ECG electrodes are embedded inside the vest and thus cannot be seen (unless the vest is turned inside out). This image was created with ChatGPT^™.

In addition to the IMU with accelerometers and gyroscopes measuring acceleration and angular velocity in all three directions, the pod also contains a global navigation satellite system (GNSS) receiver together with a local positioning system, similar to the one used in the ball. It also contains a triaxial magnetometer that measures the magnetic field direction for orientation correction, as well as antennas, RF transmitters, microcontrollers, batteries, and other electronic components.

You can imagine the importance of the thermal management design of the pod, with all of these electronics packed inside. The vest itself probably impedes cooling of the player, and if the pod is not well designed, it may even contribute as a nasty heat source.

A COMSOL Multiphysics^® simulation of a capacitively actuated surface micromachined accelerometer, which is typically used in consumer products.

Nonetheless, the data extracted from the vest is valuable for improving the players’ performance over time. The schematic below shows the output from a fictional player’s vest, which coaching staff can view live during the game and analyze afterward.

If a player is, for example, taking fewer sprints with lower acceleration and showing other signs of fatigue, the coach may decide to substitute the player. After the game, the player and coaches may analyze how the player moved across the pitch using the heat map and may also make tactical adjustments.

The data can also reinforce tactics that work well. For example, if the heat map shows good possession statistics and a high number of offensively successful crosses from one side of the pitch (“field”) for a central midfielder, the coaches may decide that the player should prioritize attacks on this side while remaining more central when defending. This might be the left side for a left-footed player, for example.

Since the vests were introduced, their statistics have been kept internal to the team and not shown on TV — much to my frustration. Sometimes, when a player is substituted, we get to see the total distance they covered, but that is usually all. The ball and vest data are not normally processed together either. Statistics such as player contact with the ball, expected goals, and shots on goal are usually obtained through image processing independently of the vest and ball data.

Imagine being able to see the full statistics from the ball and vests live during the game! Being able to visualize the most successful and least successful parts of a player’s game in real time, and not just as a table shown after the match or during halftime, would be fascinating.

Schematic view of the “UI” presenting the data from the vest of player “John Doe”. The heat map reveals that the player is an attacking central midfielder, since there are few defensive high-speed sprints. The heat map almost always shows the attacking direction from left to right. This image was created with ChatGPT.

The Crowd

So far, I’ve talked about the sensors inside the ball, how they make the VAR team’s job easier, and how they can also make the game more interesting for us viewers.

While this latest tech is changing the game for viewers, one aspect of football that always has and always will make it so exciting to watch is the crowd. And the best way to experience football is, of course, live at the stadium.

Some of the 2026 World Cup matches will be played at the iconic Estadio Azteca (officially “Estadio Banorte”) in Mexico City, the stadium where Carlos Alberto scored in the 1970 World Cup final with a fantastic thunderous instep shot from a pass by Pelé and where Diego Maradona scored the Goal of the Century (“el barrilete cósmico”) during the 1986 World Cup.

Estadio Azteca has been extensively renovated for the World Cup (and temporarily renamed “Mexico City Stadium” by FIFA). Among other things, a completely new sound system has been installed to improve the live experience. For fun and to help imagine how it sounds in the stadium, the COMSOL team simulated the acoustic wave propagation of one of the system’s speaker arrays (below); stay tuned for the third blog post in this World Cup series, where I’ll share more simulation results as I sound off on the acoustics of the beautiful game.

Simulation of the acoustic wave propagation from one of the speaker arrays hanging from the roof in the newly renovated Estadio Banorte (Estadio Azteca). The full sound pressure distribution can be obtained by superposition of the sound fields from roughly 70 distributed speaker arrays suspended from the roof structure and placed throughout the stadium.

Until then, let’s continue rooting for our favorite teams in what will hopefully be the best World Cup ever!

For the Love of the Game (Only!)

Although the models and simulations presented here are state-of-the-art, they were created just for fun. A serious scientific study would investigate the involved parameters in much greater detail. For example, the geometry of Estadio Banorte would need to be modeled in far greater detail and the simulation results validated against measurements.

These investigations were performed independently of Adidas, Kinexon, and TDK, and we do not claim any cooperation with any of these organizations.

Adidas and Trionda are registered trademarks of adidas AG.

COMSOL AB and its subsidiaries and products are not affiliated with, endorsed, by, sponsored by, or supported by any of the foregoing trademark owners.

Underground hydrogen storage is becoming essential for the global energy transition. Lined rock caverns (LRCs) provide a flexible and geographically adaptable solution, but their safety depends on complex interactions between hydrogen gas, structural linings, and fractured rock masses. Understanding these coupled effects requires advanced numerical modeling. In this blog post, we demonstrate how the COMSOL Multiphysics^® software can be used to model LRC behavior during hydrogen pressurization.

Why Use COMSOL Multiphysics^® for LRC Hydrogen Storage Modeling?

Modeling lined rock caverns for hydrogen storage presents significant computational challenges. The system involves multiple interacting materials (hydrogen gas, steel, reinforced concrete, and fractured rock), each governed by distinct physical behavior. The surrounding rock mass contains numerous preexisting fractures, introducing pronounced geometric discontinuities and constitutive nonlinearities that strongly influence the system’s response. Deformation of the fractured rock mass interacts closely with the cavern structure, governing stress redistribution and overall system stability. Under high-pressure hydrogen exposure, the steel lining may be susceptible to embrittlement, while the concrete lining may develop cracking that redistributes stresses and influences the deformation of the steel lining. These tightly coupled processes span multiple spatial and temporal scales, making realistic analysis a demanding multiphysics problem.

COMSOL Multiphysics^® is well suited for solving such multiphysics problems due to its exceptional capabilities for:

• Simultaneously solving fully coupled multiphysics equations, enabling direct interactions between mechanical, hydraulic, and transport processes
• Defining model parameters as functions of other field variables, enabling indirect couplings such as stress-dependent material properties or damage-driven stiffness evolution
• Explicitly representing discrete fractures within rock masses and resolving nonlinear hydromechanical processes within complex fracture networks
• Handling multiple interacting materials within a unified framework, supported by an extensive library of built-in constitutive models as well as flexible user-defined material formulations
• Implementing custom governing equations and application-specific constitutive models, enabling advanced formulations such as hydrogen embrittlement models to be directly incorporated into the analysis

Below, we discuss the process for building numerical models in COMSOL Multiphysics^® for LRC analysis and give a simulation example.

Overview of Modeling Steps

When using COMSOL Multiphysics^®, there are three key simulation steps:

1. Generating the geometry and mesh
2. Coupling parameters and implementing material properties, boundary conditions, etc.
3. Calculating the solution

Let’s go over these steps in more detail.

Geometry and Mesh

First, the model geometry is constructed following a multiscale strategy. On the large scale (Figure 1a), a two-dimensional domain is defined to represent a horizontal cross section of the lined rock cavern embedded within a fractured rock mass. The cavern geometry, including the steel lining, reinforced concrete layer, and surrounding shotcrete, can be created directly in COMSOL Multiphysics^® or imported from CAD software. Discrete fracture networks are geometrically represented as line segments within the rock domain and can be generated in COMSOL^® or using external tools such as MATLAB^® or CAD software and then exported as DXF™ files for direct import into COMSOL^®. On the small scale, a dedicated model of the steel lining is constructed to resolve hydrogen diffusion and embrittlement processes. The two models are coupled at the steel–concrete interface by enforcing displacement compatibility, enabling consistent interaction between structural deformation and material degradation.

Once the geometry is defined or imported, the large-scale domain is discretized using an unstructured mesh of triangular finite elements generated via Delaunay tessellation (Figure 2). Mesh refinement is applied near the cavern boundary and around fracture intersections to accurately resolve stress concentrations and damage evolution. To represent natural fractures, joint elements between neighboring finite elements are implemented, enabling explicit analysis of nonlinear fracture slip and opening within the rock mass. Thin structural components such as the steel lining and shotcrete layer are modeled using interface elements in the large-scale model in order to maintain computational efficiency. In the small-scale model, however, the steel lining thickness is explicitly represented (Figure 1b) and discretized using solid elements to resolve hydrogen diffusion and embrittlement processes across the thickness. This multiscale discretization strategy makes it possible for the global structural response to be captured efficiently while locally resolving material degradation mechanisms within the steel lining.

Figure 1. A schematic showing the model design and boundary condition of the multiscale model, including (a) a large-scale model representing an LRC situated in a fractured rock mass and (b) a small-scale model capturing the response of the steel lining subject to cyclic internal pressurization and boundary displacement constraints, as well as hydrogen diffusion. The dimensions and scales shown follow the model configuration adopted in our previous study (Ref. 1) and are presented here to illustrate the model setup strategy. The actual model configuration may vary depending on specific applications and site conditions.

Figure 2. Mesh discretization for the LRC model (Ref. 2). The reinforced concrete layer is discretized using structured quadrilateral elements, while the rock matrix is discretized using unstructured triangular elements. Steel lining and shotcrete are discretized using line elements, whereas fractures are discretized using joint elements.

Interfaces and Couplings

We use the Solid Mechanics interface in COMSOL Multiphysics^® to simulate deformation of the fractured rock mass, concrete lining, and steel components under internal hydrogen pressurization. The rock matrix is modeled using a scalar damage formulation to capture crack initiation and propagation, while discrete fractures are implemented as joint elements with nonlinear normal and shear constitutive relations to represent fracture slip and opening.

The reinforced concrete lining is modeled using the Mazars’ damage model, available in COMSOL^®, which captures tensile cracking and progressive stiffness degradation. The influence of reinforcement is incorporated through an enhanced elastic modulus and residual strength parameters, consistent with the adopted concrete class.

In the large-scale model, the steel lining is represented using interface elements due to its thin geometry. Its mechanical behavior follows an elastoplastic constitutive law with exponential hardening. To capture hydrogen embrittlement effects, a separate small-scale model explicitly represents the steel lining thickness using solid elements. In this model, a customized interface is developed with the Physics Builder in COMSOL^® to simulate the hydrogen diffusion across the lining while accounting for the coupling with solid mechanics. The mechanical behavior of steel is coupled to hydrogen concentration by defining the yield stress as a concentration-dependent variable, enabling analysis of hydrogen-induced strength degradation.

When pore pressure effects in the fractured rock mass are considered, the Darcy’s Law interface in the Subsurface Flow Module is incorporated to simulate fluid flow within the rock matrix and fractures. The Poroelasticity coupling can be activated to achieve direct coupling between mechanical deformation and pore pressure evolution. In this case, hydraulic properties such as permeability or fracture aperture may be defined as stress-dependent variables, enabling indirect hydromechanical coupling consistent with our previous modeling framework.

Material properties and constitutive equations are defined separately for the rock matrix, fractures, concrete, and steel. Direct multiphysics coupling ensures consistent interaction between deformation and hydrogen transport (and pore pressure when activated), while additional indirect couplings are introduced by defining model parameters as functions of evolving field variables such as stress, damage, or hydrogen concentration. Mechanical boundary conditions include in situ stresses applied at the outer rock boundary and internal hydrogen pressure applied at the cavern wall.

Calculating the Solution

The analysis is performed in two stages: In the first stage, the large-scale model is brought to equilibrium under the prescribed in situ stresses using a ramped loading procedure. In the second stage, hydrogen pressurization is applied at the cavern boundary, either monotonically or cyclically, to simulate storage operation. The resulting displacement field from the large-scale model is imposed on the small-scale steel model, where hydrogen diffusion and concentration-dependent mechanical degradation are solved in a time-dependent manner. Nonlinear solution schemes are used to resolve fracture reactivation, damage evolution, and plasticity.

LRC Simulation Example

We apply the multiscale model to simulate cyclic hydrogen pressurization of an LRC embedded in a fractured rock mass (Ref. 1). The large-scale modeling results show cyclic radial displacement of the concrete lining and progressive damage development in both concrete and surrounding rock (Figure 3). Damage localizes primarily in tensile regions and near fracture intersections (Refs. 1–2), highlighting the strong control of fracture distribution in rock on the LRC’s structural response (Figure 4).

The small-scale model captures hydrogen diffusion and embrittlement within the steel lining (Figure 5) (Ref. 1). Hydrogen concentration increases from the inner surface and evolves over loading cycles, leading to local strength degradation that correlates with stress concentration zones. These results demonstrate the coupled interaction between cyclic pressurization, fracture reactivation, stress redistribution, and hydrogen-induced degradation across scales.

Figure 3. Simulation results (Ref. 1) showing the distribution and evolution of (a) radial displacement in the concrete lining and (b) damage development in the concrete and surrounding rock over multiple hydrogen pressurization cycles (cycle period T0 = 24 h).

Figure 4. Simulation results (Ref. 1) showing the distribution of damage and local maximum principal stress in the vicinity of the LRC during cyclic hydrogen gas pressurization (cycle period T0 = 24 h).

Figure 5. Simulation results (Ref. 1) showing the spatial distribution of (a) hydrogen concentration, (b) maximum principal stress variation, and (c) strength degradation in the steel lining at different loading stages.

In addition, the framework has been extended to include hydromechanical coupling in the fractured rock mass (Ref. 3) and time-dependent rock creep (Ref. 4). These extensions enable us to evaluate how fluid pressure diffusion and viscoelastic deformation in the surrounding rock mass influence the LRC’s long-term performance.

References

1. C. Zhao et al., “Modelling lined rock caverns subject to hydrogen embrittlement and cyclic pressurisation in fractured rock masses,” International Journal of Hydrogen Energy, 2025; 152: 150027.
2. C. Zhao, Q. Lei, Z. Zhang, “Impact of fracture networks on the structural deformation of lined rock caverns under high internal gas pressure,” Underground Space, 2025; 21: 252-269.
3. C. Zhao, Z. Zhang, Q. Lei, “Coupled hydro-mechanical simulation of the interaction between adjacent lined rock caverns subject to internal gas pressurisation,” Geomechanics for Energy and the Environment, vol. 43, 2025: 100701.
4. C. Zhao et al., “Influence of rock creep on the performance of lined caverns under cyclic pressurization and hydrogen embrittlement,” International Journal of Rock Mechanics and Mining Sciences, vol. 199, 2026; 106401.

About the Author

Qinghua Lei is an associate professor at Uppsala University, Sweden. He earned his BSc (2009) and MSc (2012) in civil engineering from Tongji University, China, and his PhD (2016) in rock mechanics from Imperial College London, UK. He worked as a postdoctoral researcher in fluid mechanics at Imperial College London (2016–2018) and later as a senior researcher and lecturer in engineering geology at ETH Zurich, Switzerland (2018–2023). He is the recipient of the 2025 ERC Consolidator Grant, 2024 Chin-Fu Tsang Coupled Processes Award, 2019 Rocha Medal, 2016 NGW Cook PhD Dissertation Award, and 2015 Rock Mechanics Research Award. He is a fellow of the Young Academy of Europe. His research interests include rock mechanics, hydrogeology, geophysics, natural hazards, and geotechnical engineering.
 
MATLAB is a registered trademark of The MathWorks, Inc. Autodesk, the Autodesk logo, AutoCAD, and DXF are registered trademarks or trademarks of Autodesk, Inc., and/or its subsidiaries and/or affiliates in the USA and/or other countries.

Understanding the Pain Points of Data Transfer

We are frequently asked how data that is in different formats can be brought into COMSOL^®. This data often lives in text files, and the file formats can be unique, with no translation tools available. Our customers may understand those formats quite well, while we have a strong understanding of how such data should be used within COMSOL^®. In this situation, AI can significantly reduce the effort needed to implement a data translation workflow. Let’s look at this in the context of a specific scenario.

A Sample Scenario

Suppose we have been given a text file containing a description of a lumped thermal model. A lumped model is a bit different from the finite element models you may be familiar with in COMSOL Multiphysics^®. Lumped models use nodes that represent volumes of material. Each volume has a known density and specific heat, so each node has an associated thermal mass. The temperature of a node can increase due to an applied heat load, and heat can flow between nodes by either conduction or radiation, computed based on the temperature difference between two nodes. We can also consider nodes that are at a fixed temperature, called temperature sinks. To make this more concrete, we will put together a very simple lumped model of a satellite in deep space.

The satellite model that we will work with is shown in the figure below. It is a box structure with two solar panels protruding from either side. The six sides of the box and the two solar panels are each represented by a single thermal node. There is conductive heat flux between adjacent faces of the box structure and radiative heat flux between the side faces and the solar panels. There is also radiative heat flux from all faces to deep space. The objective of this model is to compute temperature over time, starting from an initial temperature. We will assume that solar and planetary loads are negligible, such as when a geostationary satellite goes into eclipse, and that the only load is due to a heater on one node.

Schematic of an eight-node lumped thermal model of a satellite. There is conduction between the nodes of the box, as well as radiation from all nodes to deep space.

These types of lumped models are quite common, and a number of different file formats are used to represent them. To keep the example format-agnostic, we will generate a file that is similar in spirit, if not in exact syntax, to these formats. As we will see, the exact syntax is not the main point. An excerpt of our sample file is shown below:

TSINK, 99, 2.7 # Thermal Sink 99
NODE, 10, 300, 12150 # Node 10
…
NODE, 80, 300, 3220 # Node 80
LOAD, 60, 750 # Load on node 60
CON, 1020, 10, 20, 0.71 # Conductor 10 - 20
…
CON, 5020, 50, 20, 0.71 # Conductor 50 - 20
RAD, 2070, 20, 70, 8.9e-9 # Radiation 20 - 70
…
RAD, 8099, 80, 99, 1.8e-7 # Radiation 80 - 99

There are five types of data in this file: TSINK, NODE, LOAD, CON, and RAD. These represent temperature sinks, temperature nodes with mass, thermal loads on particular nodes, conductive connections between nodes, and radiative connections between nodes or between a node and a sink. The information after the # symbol is a comment. If you are familiar with similar types of files, the above format can be described as a new dialect in the computer science sense.

To borrow a few more phrases from computer science, the information in this file describes a connected graph, where the records describing conduction and radiation are edges between the nodes. This is particularly useful, since AI tools are good at dealing with these types of data structures. Let’s keep that in mind for later, but now let’s turn our attention to getting this data into COMSOL^®.

The COMSOL Equivalent

Although COMSOL Multiphysics^® does include a Lumped Thermal System interface as part of the Heat Transfer Module add-on, we want something with a little less overhead for larger models. The simplest equivalent for representing the data above is to use the Global Equations interface. Since the data we are trying to import represents a three-dimensional structure, but the input file does not contain enough information about the shape and dimensions of the geometry, representing the model in an abstract format is justified.

The Global Equations interface enables us to represent the graph network as a coupled system of ordinary differential equations, so equations of the form:

 
These equations can be represented within the user interface as shown in the screenshot below. We can also write out, and read in, a whole set of these global equations using the Save to File and Load from File buttons.

Screenshot showing how to implement a set of ordinary differential equations. Note also the Save and Load buttons.

The model that we’re trying to reproduce here can be set up with nine entries for the nodes and the sink. I put these together by hand. It was a bit tedious, but I knew that I would only need to do it once, so it was worth the effort. Note that the original file was written with an assumed set of units, so all equations were nondimensionalized. After entering this into the user interface, I wrote the model out to a text file. A few lines are shown below:

NODE10 12150*d(NODE10,t)[s]-((0.71*(NODE20-NODE10)+0.71*(NODE30-NODE10)+0.71*(NODE40-NODE10)+0.71*(NODE50-NODE10))+(4.5e-8*(TSINK99^4-NODE10^4))) 300 0 "Node 10"
...
NODE60 12150*d(NODE60,t)[s]-((750)+(0.71*(NODE20-NODE60)+0.71*(NODE30-NODE60)+0.71*(NODE40-NODE60)+0.71*(NODE50-NODE60))+(4.5e-8*(TSINK99^4-NODE60^4))) 300 0 "Node 60"
...
TSINK99 TSINK99-2.7 2.7 0 "Thermal Sink 99"

Using AI to Build the Translator

At this point we should have a good understanding of these two files. We should also understand that a file with hundreds or thousands of entries will require some level of automation. What may not be so obvious is that we are converting graph data from an edge-list representation to an incidence-list representation that exploits sparsity. Readers with a little background in computer science will recognize that this type of conversion algorithm is not trivial to implement.

Fortunately, these are exactly the kinds of tasks that AI can help automate. To be clear, we are not going to ask AI to do the conversion; we are going to ask AI to write a general-purpose conversion tool. But what do we need to ask for? How do we instruct the AI tool to write the conversion tool for us?

We already have almost all of the data that we need within these two text files. All we need to do is upload them to the AI tool of choice and use a prompt along these lines:

Here is a file that contains nodes and sinks, representing nodes on a graph network. The connections between the nodes are defined by CON and RAD lines. Everything after a # is a comment. I need to convert the first file into the format of the second file. Please write out the transformation rules between them.

In this case, I used ChatGPT. After about a minute, it presented a full human-readable description of how to do the conversion. It was particularly impressive to me that the AI tool recognized that the file referred to a thermal problem and identified the nature of the nonlinear radiative term.

Excerpt of the ChatGPT output.

As I spent some time reviewing this output in detail, I gained confidence that the tool understood the task. Keep in mind that the input file uses a unique, undocumented dialect for describing a lumped thermal model. I needed only one more step: asking the AI tool to generate code that I could run independently. This required one more prompt:


Please write a monolithic piece of Java code that I can use in COMSOL's Method Editor to convert any file of the original type into this format.


I asked for Java code here because I wanted to run this entirely within COMSOL Multiphysics^®. Just as ChatGPT had no difficulty understanding the unique dialect of the input file, it also had no difficulty writing code in the programming language I requested.

The resulting code was several hundred lines, so I will not show it here. I did, however, review it briefly and noticed that ChatGPT added comments for readability and even included error checking so that it would fail with an informative message if given an invalid file. The code compiled and ran without issues, and I was able to do some preliminary checking by generating results.

In general, it is important to verify and validate translation code, and the right approach depends on the form of the data. Sometimes a visual spot check of a few cases may be sufficient; other times, you may need to be more rigorous. One approach is to write a back-converter as well, essentially round-tripping the data and verifying that it still matches the original. We can use AI for this step as well.

Sample output computed from the input file.

At this point, some readers may be asking: Why not just incorporate this translator into the COMSOL product suite? Let me emphasize that this example is a minimal thermal network file in a unique dialect. More practically, lumped thermal model files can incorporate many more record types, require collections of multiple files, and even contain customized code that implements, for example, a specific type of thermostatic control algorithm. Writing a general-purpose translator could involve covering many edge cases that might never arise in the same file. Furthermore, one-to-one translation is not always the best path. Sometimes you need to take a step back and understand the specific modeling intent of the person who created the original file. That is particularly true of these types of thermal network models, but that is a topic for another day.

Remarks on File Translation and AI in General

It is worth reemphasizing how little effort was required to get to this point. We already had the sample input file. We did spend some time figuring out how best to represent this within COMSOL^® and generating sample syntax, but the interaction with the AI tool was minimal: Two files were uploaded, two prompts were given, and the results were applicable to any file using the same format. The total interaction time with the AI tool was a few minutes.

If you need to perform this type of structured data or file translation, you can become familiar with this workflow quickly. What if you need to incorporate more features? Try adding a few more lines of sample input and output, and then prompt the AI tool again. What about other data that you want to automatically bring into COMSOL? Your company may have a large proprietary material database in a custom format that you need to import. Use this workflow to write the translator. Keep in mind that the proprietary data itself does not need to be shared with the AI tool; only the data format is shared. What you get back is an algorithm and source code.

More broadly, as discussed in a previous COMSOL blog post, AI tools continue to improve rapidly. For the kinds of problems that COMSOL users need to solve, these capabilities are becoming increasingly useful.

For now, I will leave you with the results of one final prompt, where I asked the AI to draw its interpretation of the hardware represented by this file. The result is shown below. It is not quite there yet, but I will be sure to revisit this in a few months to see how things have changed.

How AI interprets the global equations.

Next Step

How might you use this workflow to help with your COMSOL modeling? Leave your thoughts below, or contact us!