Pulse-echo measurements are a standard application in the oil industry for detecting material properties behind pipes. The measurement setup is simple, but modeling is challenging due to the high frequencies of ultrasonic pulses and the complexities of 3D modeling for this type of application. A 2D axisymmetric simplification can be a valuable step that reduces computation time and allows for parametric studies. In this blog post, we discuss how we used this type of simplification to improve our pipe simulations.

Pulse-Echo Measurements in Pipes

In the oil industry, pulse-echo ultrasonic measurements are important for obtaining information about the material properties behind an oil pipe and about the bonding of the material to the pipe. Relevant considerations are, for example:

• What the quality of the cement behind the pipe is
• If the shale behind the pipe is bonded to the pipe
• If there is a fluid-filled gap between the pipe and the solid surrounding material

These considerations are important before and during production, but also during plug and abandonment operation when closing an oil field.

Pulse-echo measurements are used to derive the material properties outside of a pipe by measurements from inside the pipe.

From the circular transducer, a short Gaussian pulse is sent with normal incidence toward the pipe wall (Figure 1). Once the signal reaches the pipe wall, part of it is reflected back into the fluid toward the transducer while the other part is transmitted further through the pipe wall, to the outer material. The pulse is then reflected back and forth inside the pipe wall. Every time the pulse is reflected on the inner pipe wall part of the signal is transmitted toward the transducer.

These returning signals are recorded with the same transducer that is used to send the initial pulse. The decaying strength of the signal inside the pipe wall depends on the material parameters, specifically the impedances of the material inside the pipe, the pipe itself, and the material outside the pipe. The decay rate of this signal measured at the transducer can be used to estimate the material properties outside of the pipe if the properties of the pipe and the material inside the pipe are known.

Figure 1. Concept of pulse-echo measurements in pipelines (left). An animation showing the transmitted and reflected signals (right).

Modeling this type of pulse-echo measurement in pipes is challenging due to the distance between the transducer and pipe, the high frequencies of the ultrasonic pulse, and the 3D geometric setup of a circular transducer within an elongated pipe. A 2D axisymmetric simplification is carried out with respect to the transducer symmetry axis, which proved to be a valuable step for the parametric study we were carrying out.

Through this model, we were able to reduce the computation time significantly, allowing us to build a database with about 1400 simulations including variations in pipe thickness and curvature, material parameters, and transducer–pipe distance, among other factors. These simulations allow for further research of the pulse-echo method, its sensitivities, and possible improvements to the interpretation of this type of data.

The COMSOL Model

To model this pulse-echo measurements setup correctly requires a 3D model. The transducer is circular, making it axisymmetric with the x-axis, while the pipe is axisymmetric with the z-axis (Figure 1). Hence, we can make use of two symmetry planes for this model, reducing the model domain to a quarter. Figure 2 shows the 3D model with the two implemented symmetry planes at four time steps. The ultrasonic pulse has propagated from the transducer surface (12 µs) to the pipe (36 µs), exciting acoustic waves in the pipe (52 µs), and the reflected pulse has traveled back to the transducer, where the signal is recorded (73 µs). Energy from the vibration within the pipe transmitted toward the transducer is visible behind the initial pulse in the snapshot at 73 µs. The modeling domain is surrounded by absorbing layers to ensure that no reflections of the waves occur at the domain boundaries. The xy-plane and the xz-plane are symmetric planes.

Figure 2. 3D modeling results at four time steps, making use of two symmetry axes.

Parameters of a standard pipe and transducer geometry, common in industry applications, were implemented for the model shown in Figure 2. These parameters can be found in the table below:

The time-dependent study step is computed from 0 to 140 µs to allow the propagation of the ultrasonic pulse to the pipe and back as well as the recording of a significant part of the pipe’s reverberations.

Use of Time-Explicit Domain

The pipe in this model is normally filled with some kind of fluid. Here, we used oil-based mud with longitudinal wave velocities of 1301 m/s. The pipe itself is of steel, with a longitudinal wave speed of 5800 m/s. We calculated the element size for the mesh using:

with the element size , the material’s wave speed , and the maximum frequency . Due to the large difference in the speed of sound in steel and oil-based mud, there are significant differences in the required mesh size for these domains. Figure 3 shows the mesh size for a 2D plane. The most efficient way to model this was by using the time-explicit mode. The acoustic domain of the fluid-filled pipe and the elastic wave domain of the pipe and surrounding material were then coupled by identifying the respective surfaces as identity pairs and using the Pair Acoustic–Structure Boundary multiphysics coupling.

The domain is surrounded by absorbing layers to prevent reflections from the domain boundaries. Here, two absorbing layers were defined, one for the acoustic domain and one for the elastic wave domain. 

Geometry Reduction

To be able to carry out parametric studies with hundreds of variations, it is important to have a model with a relatively short computation time. Hence, using the 3D model with a computation time of more than 7 hours on our machine for each model is unrealistic. Therefore, we explored the possibility of reducing the geometry dimensions.

A standard way to reduce the computation time is by going from a 3D simulation to a 2D simulation. Taking a slice in the xy-plane and making use of the pipe’s mirror symmetry makes it possible to reduce the model significantly. In this 2D case, the pipe is modeled correctly, as the 2D model assumes infinite extend in the third direction. However, the source also becomes infinite in the third direction and is therefore not modeled correctly. This leads to significant deviations between the model results from 2D models compared to 3D models for this measurement geometry (Ref. 1).

Another possibility to reduce the geometry dimension is the use of a 2D axisymmetric model. As pointed out, the transducer is axisymmetric with respect to the x-axis, while the pipe is axisymmetric with respect to the z-axis (Figure 1). Here, we chose to carry out the axisymmetric model so that the symmetry axis aligns with that of the transducer (Figure 3). Thus, the geometry of the circular transducer is modeled correctly. Choosing again the xy-plane for the modeling of the 2D-axisymmetric part means we are modeling the curvature of the pipe. However, applying the symmetry about the x-axis means that the pipe curvature is assumed to be axisymmetric, hence modeled as a part of a spherical shell (Figure 4).

Figure 3. Building the model and grid for the 2D axisymmetric model. The time-explicit study step is used. The acoustic and solid domains are coupled by an identity pair.

Figure 4 shows the results of the wave propagation for the 2D axisymmetric model. The modeling is carried out for the domain shown in Figure 3, and the results in Figure 4 are plotted with the radial extension. The four presented time steps are the same as for the 3D model. The measured signal integrated over the transducer surface from the 2D axisymmetric and 3D model are plotted in Figure 5. (A detailed discussion of the comparison of 3D, 2D, and 2D axisymmetric models and the justification for using 2D axisymmetric models can be found in Ref. 1.)

Figure 4. 2D axisymmetric modeling results in the revolved geometry for four time steps. The modeling domain is surrounded by absorbing layers.

Building a Database

The solution time for the 2D axisymmetric model was around 13 minutes, a significant improvement from the multiple hours of computation time for the 3D model. This speedup made it possible to build a database with hundreds of variations in the model (Ref. 2). Beside geometrical variations, we also introduced a fluid-filled annulus between the outside of the pipe and the surrounding solid material in the range of 10 to 1000 µm. We did this by making use of the advantages of the time-explicit implementation, using additional Pair Acoustic–Structure Boundary couplings for the transition between the fluid and solid domains. For each calculated model, the pressure at the transducer was exported and integrated over the transducer surface, giving the results of the measured signal from the pulse-echo modeling for further analysis and research (Figure 5).

Figure 5. COMSOL model results of the pressure over the transducer surface.

To simplify the calculation of all these models, we used LiveLink™ for MATLAB^®, which enables the integration of COMSOL Multiphysics^® with MATLAB^® By doing so, we were able to drive the variation of all the calculations we were interested in automatically and export the pressure over the transducer surface. The input information and the averaged pressure over the transducer surface were then written into a JSON file. These results make up the database, which can be used for further analysis.

Access the Data and Models

To further explore the model discussed in this blog post, download it via the following Application Exchange entry: Modeling pulse-echo ultrasonic data.

The data generated within this project, and the 2D axisymmetric and 3D COMSOL^® models (v6.2), are available on Mendeley data (DOI:10.17632/3bs65nzpv2.1).

References

1. A. Diez, T.F. Johansen, E.M. Viggen, “From 3D to 1D: Effective numerical modelling of pulse-echo measurements in pipes,” Proc. 46th Scandinavian Symposium on Physical Acoustics, pp. 1–23, 2023; ISBN 978-82-8123-023-1.
2. A. Diez, E.M. Viggen, T.F. Johansen, “Ultrasonic pulse-echo dataset from numerical modelling for oil and gas well integrity investigations,” Sci Data 12, 544, 2025; https://doi.org/10.1038/s41597-025-04851-x

Acknowledgement

This work was a collaboration between A. Diez, T.F. Johansen (SINTEF), and E.M. Viggen (Norwegian University of Science and Technology) for the Centre for Innovative Ultrasound Solutions, funded by the Research Council of Norway under grant no. 237887.

About the Author

Anja Diez is a researcher in the Acoustics Group at SINTEF. She has a background in geophysics and worked the first years of her career using seismics and ground-penetrating radar to investigate the ice sheets in Antarctica and Svalbard, determining ice properties and conditions at the glacier bed. In recent years, she has worked on acoustic projects related to industrial applications and nondestructive testing, combining signal processing and data analysis with COMSOL^® modeling for wave propagation in fluids and solids.

What AMPP 2026 Was All About

If there’s one thing that was immediately clear when walking the floor at AMPP 2026, it’s that corrosion is a much bigger deal than most people outside the field realize. We’re talking about a challenge that touches oil and gas pipelines, power systems, transportation infrastructure — pretty much anything made of metal that has to survive harsh conditions over a long period of time. What struck me most was how diverse the crowd was. You had senior engineers deep in the technical weeds right alongside procurement teams, environmental managers, and operations leaders, all of them there because corrosion affects their work, their budgets, and their safety responsibilities. The conversations happening on the show floor weren’t just technical; they were about risk, cost, and how to do things smarter. In this vein, we noticed a theme that kept coming up: People want better ways to predict and prevent corrosion before it becomes a problem, not just ways to respond to it after the fact.

COMSOL at the Exhibit Hall

The crowd stopping by our booth was genuinely all over the map: engineers, operations managers, procurement folks, and everyone in between. But what connected all of them was the same underlying challenge: Physical testing alone is either expensive or slow and doesn’t always give you the full picture. Physics-based simulation, on the other hand, enables engineers to virtually analyze real problems — from how quickly a pipeline might corrode to what protection strategy works best — before prototypes are built or physical tests get underway. That’s a value proposition that resonates whether you’re an engineer solving a technical problem or a manager trying to reduce costs and risk. At our booth, we were able to demonstrate our simulation capabilities live and have meaningful conversations about how virtual testing can complement, or even replace, costly physical tests.

Chen You, Bertil Nistad, and myself were on hand at the COMSOL booth to answer questions about the COMSOL^® software and demonstrate its capabilities.

COMSOL’s Technical Presentation

One of the things I was most excited about was getting to watch two of my colleagues present in the technical program. COMSOL Senior Applications Engineer Chen You and Development Manager Bertil Nistad delivered a session together on modeling corrosion prediction and cathodic protection in the COMSOL Multiphysics^® software, and the room was engaged from start to finish. Even coming from a nontechnical background, I found myself genuinely following along. They broke down how simulation can be used to model the conditions that lead to corrosion and how products like the Corrosion Module add-on to COMSOL Multiphysics^® help engineers design better cathodic protection systems. The Q&A at the end reflected this engagement, with attendees asking about how the models are applied in practice, how cathodic protection systems are simulated, and how the capabilities of the COMSOL^® software can be used in different scenarios.

As part of the AMPP 2026 technical program, Chen You and Bertil Nistad gave a presentation on using the COMSOL^® software to model corrosion prediction and cathodic protection.

“The AMPP Annual Conference is a great place to meet our existing customers and showcase our modeling software to a wide audience of possible prospects,” Nistad said. “For me, it was a week full of great presentations and lively discussions about corrosion modeling.”

Nistad also had the chance to be interviewed at the conference by Inspenet. You can check that out here.

Why It Matters

Here’s the thing about corrosion: It costs the global economy trillions of dollars every year. That’s not a typo. Whether they stem from a leaking pipeline, a degraded coating on a ship hull, or a power component failing ahead of schedule, the downstream effects of corrosion are real: safety risks, unplanned downtime, and major costs that hit both operators and the environment.

Being at the AMPP Annual Conference + Expo made it easy to see why simulation is becoming such an important part of how the industry works. When businesses can model a problem and virtually test different materials, protection strategies, or designs, they make better decisions before anything is built. That’s good for engineers, good for business, and good for the environment.

Model Corrosion Processes

Want to learn more about how modeling and simulation can be used for understanding corrosion and designing and optimizing corrosion protection systems? Check out the Corrosion Module via the button below.

Why Induction Heating is a Modeling Challenge

Induction heating is widely used in materials processing applications for many reasons, including:

• The workpiece is indirectly heated.
• The zone being affected can extend into the volume of the material.
• The heating profile over time can be very precisely controlled by altering the magnitude and frequency of the driving currents.

Precise control over the heating process, however, does require some kind of process feedback loop and experimental verification. This loop should ideally be augmented with a numerical model that incorporates experimental data and can predict the temperature evolution over time. A numerical model is especially useful since it can predict behavior in regions that are not directly measurable in real time.

Practically speaking, numerical models do pose several challenges. You need good input data, starting with an accurate description of the geometry of the workpiece and induction heating coil, as well as knowledge of the driving currents and frequencies. From there, you also need to know the material properties and how they vary.

All material properties vary with temperature, and since induction heating usually starts with the workpiece at around room temperature and raises it up to around the melting temperature, these nonlinearities can never be neglected. Both the electrical and thermal conductivity of metals usually decrease with increasing temperature, but not always.

Representative B–H curves at different temperatures. Above the Curie temperature, the slope of the curve equals the magnetic permeability of free space.

Furthermore, there are some materials, specifically soft magnetic materials, that have magnetic permeability that is nonlinear with magnetic field but with neglectable hysteresis. These materials can be experimentally characterized by a B–H curve, relating the magnitude of the -field to the -field. These curves also exhibit a change with temperature, decreasing gradually in magnitude up to the Curie temperature, at which point the material becomes nonmagnetic. Although the form of such curves will usually be quite similar to what is shown above, the actual experimental data can be difficult to acquire and will often represent the biggest hurdle to implementing a numerical model. 

Implementing the Model: A Few Conceptual Preliminaries

Induction heating processes usually occur over a timespan lasting at least a few seconds, while the excitation frequencies can range from as low as 50 Hz up to 450 kHz or even beyond. This means that there are two temporal scales to the problem. The temperature fields can be assumed to vary quite slowly relative to one period of oscillation. In other words, the electromagnetic fields do not observe a variation in temperature over one cycle. This motivates a solution method termed the Frequency-Transient study approach, where the electromagnetic fields are solved in the frequency domain, while the thermal problem is solved in the time domain.

The electromagnetic fields do get recomputed but only when the material properties vary with respect to a change in temperature. The explicit assumption of a frequency domain analysis is that the material properties are constant over one cycle. However, in practice, the electromagnetic fields usually do vary nonlinearly over one cycle. We are usually trying to heat the material as quickly as possible, meaning that the excitation currents and resultant fields will be driving the material into the nonlinear part of the B–H curve.

To resolve this apparent contradiction, we use the concept of the effective B–H curve, which uses an energy-based method to compute a curve that approximates the nonlinear behavior. These curves can be computed by starting with the experimental tabulated B–H curve. The workflow begins with acquiring the B–H data, cleaning it up using the B–H curve checker, and then using the smoothed data within the effective B–H curve calculator.

Following the above concepts, you should now have at least three sets of data for the workpiece material:

• temperature
• nonlinear electrical and thermal conductivity
• a set of effective B–H curves at different temperatures for each material being heated

To set up a complete thermal model, you will also need the density and the specific heat, as well as the thermal emissivity of the surface since we will always need to consider at least radiative cooling. The specific heat also always varies with temperature. Density is always held constant — any change in volume has to be modeled by also solving for thermal expansion of the solid material. The surface emissivity is held constant in the example we will look at, although in practice it may vary.

The inductive coils are almost always water-cooled copper pipes and thus are much simpler to model since their temperature is fixed. We are usually not interested in the heat distribution within the coil itself. We only need to know how the currents in the coil are heating up the workpiece. This leads to a simpler modeling approach of using the Impedance Boundary Condition to model only the coil surface, thus avoiding the need to mesh the interior volume of the coil. This simplification can be relaxed by modeling the volume of the coil as well, if desired.

The workpiece, on the other hand, must always be modeled as a volume, since the variations in the properties cannot be accurately modeled via a surface condition. Furthermore, the fields will often vary quite sharply in the direction normal to the surface but quite gradually in the tangential direction. This motivates the usage of a boundary layer mesh. When working in 3D, it can even be useful to partition the geometry by introducing layers at the surface.

A Quick Aside: Using Analytic Equations for the Effective B–H Curve

Although experimental data is the ground truth, it is not always easy to acquire this data to high accuracy, especially over a wide temperature range. So, sometimes it motivates to use simpler expressions for the B–H curve, especially when we can use these to subsequently come up with analytic expressions for the effective B–H curve.

Two convenient expressions for the magnitude of the -field, , as a function of the magnitude of the -field, , are:

These have the same limiting behavior: at low magnetic field, the slope (the differential relative permeability) is , and in the high-field limit, the slope is , equal to free space. The saturation flux density, , is the point at which the slope approaches the free space limit.

By using this form, we can easily introduce a temperature dependency to approximate the behavior as we approach the Curie temperature, above which materials can be approximated as being nonmagnetic. We do so by making the saturation flux density a function of temperature, where the function can have any form as long as it drops to zero above the Curie temperature. Usually, it will be monotonically decaying, so we could posit a curve of the form shown in the plot below, with polynomial decrease up to the Curie temperature. Such a curve can be defined using the Piecewise function in COMSOL^®, and a smoothing term can be added for numerical stability.

Approximation of the change of the magnetic behavior with temperature. Smoothing is applied over a small region about the Curie temperature.

With these expressions, we can now derive an analytic expression for the effective B–H curve using the Simple Energy method, where:

Plugging in our two expressions from above gives:

Note that we could also introduce a temperature dependency to the term, altering both magnitude and slope of the B–H curve. It should be emphasized that these are solely presented as convenient approximations, and you could repeat this derivation for other functions and other equations for computing the effective B–H curve. In the absence of high-quality experimental data, these can be a reasonable starting point for your modeling.

Implementing a Model Using the Magnetic Field Formulation

The case we will consider here is of a steel square channel heated by a three-turn coil. We will assume that the coil pitch does not strongly affect the solution and will use this assumption to exploit the symmetry of the workpiece and coil to reduce model size.

A three-turn inductive coil around a section of a steel square channel (left) can be reduced to a one-sixteenth model (right) of the workpiece and the surrounding free space.

The approach that we will take here is to use the Magnetic Field Formulation interface, which solves for the -field directly. This makes it slightly more suitable for modeling of materials where the nonlinearity depends explicitly on the -field itself. We’ve already covered how to introduce excitations into this formulation and will use a global constraint to fix the driving current on the coil.

When entering the material nonlinearity, you can reduce the computational cost by avoiding symbolic differentiation with respect to magnetic field in the effective B–H material relationship. This is done using the nojac () operator, and although this can occasionally lead to an increased number of solver iterations, each iteration takes less time and memory, resulting in an overall improvement. Note that this holds, regardless of whether you’re using the analytic expressions derived above or using tabulated data.

The solution approach uses the Frequency-Transient solver, which will default to using a segregated approach for 3D models. The -field is solved for using a direct solver, and this is a great place to try out the new NVIDIA CUDA^® direct sparse solver (NVIDIA CuDSS) if you have an appropriate graphics card. 

The results are animated below, and show the effects of the rise in temperature. The effective relative permeability changes as the part heats up, and eventually drops to unity once the Curie temperature is exceeded.

Animation of the effective relative permeability, which is affected by field strength and temperature, over time.

With regards to the meshing, you can use the maximum operating frequency along with the conductivity and maximum effective permeability to compute the minimum skin depth. Along with that, you can investigate different element order, although when you have highly nonlinear materials, the linear order is usually preferred.

The computational costs of a three-dimensional model can of course be significant, so reducing the model by exploiting symmetry is always encouraged. To take it one step further, you can often reduce the model to two dimensions. Such models are particularly useful since they can help guide you towards an efficient meshing strategy in your three-dimensional models.

Two-dimensional models can often provide reasonable approximations of a three-dimensional model at much lower computational cost. The 2D results are mapped onto their corresponding 3D faces.

Comparison of the two-dimensional and three-dimensional model results.

The importance of building 2D models cannot be overemphasized, as they can often provide a very good estimate of the 3D solution, and it is much easier to investigate meshing, solver, and discretization strategies. As you expand these types of models, you will likely also want to include more physics, such as solving for the mechanical properties: Induction Hardening of a Cylindrical Pin.

Closing Remarks

We have shown here that it is easy to use the -field formulation to set up and solve induction heating of materials that are nonlinear with both magnetic field and temperature. This formulation can be more stable, and hence, converge faster than other approaches, particularly when you have strong nonlinearities. 

Interested in getting started in this area? Download the tutorial model from the link below.

A Quick Overview of the Different Formulations

The five interfaces that can be used to solve for time-varying magnetic fields within the AC/DC Module are:

• Magnetic Fields interface
• Magnetic and Electric Fields interface
• Magnetic Fields, Currents Only interface
• Rotating Machinery, Magnetic interface
• Magnetic Field Formulation interface

The first of these, the Magnetic Fields interface, is by far the most commonly used and is the workhorse of most inductive modeling. It solves for the magnetic vector potential, the -field, and has convenient ways to excite structures such as: lumped port boundary conditions, coil domain conditions, and background field excitations. Most people do, and should, start their modeling with this interface.

Next, the Magnetic and Electric Fields interface solves for both the -field and the electric potential, the -field, making it an – formulation and introducing additional complexity. This added cost is warranted primarily in two cases:

• When solving magnetohydrodynamic models
• When solving at line frequencies, where you need to very precisely model the losses within dielectrics and at interfaces

The Magnetic Fields, Currents Only interface is useful strictly when there are no magnetic materials present. It is primarily meant for computing partial inductances. It also solves for the -field but does so using Lagrange elements rather than Curl elements, making it less computationally expensive.

The Rotating Machinery, Magnetic interface, as its name implies, is for modeling of motors and generators. This formulation allows parts to rotate and move relative to each other by solving for both the magnetic vector potential and the magnetic scalar potential, the -field, which introduces additional cost and complexity. For an introduction, see our Motor Tutorial Series.

These formulations are all derived from Maxwell’s equations and all solve for the -field, either on its own or with an additional field. Contrast this with the Magnetic Field Formulation interface, which is, of course, also derived from Maxwell’s equations but instead solves for the magnetic field, the -field, in the following form in the time domain:

and in the frequency domain:

The solutions to these equations are also solutions to Maxwell’s equations, so why would we use this instead of the -field formulations? There are (at least) two reasons:

• When dealing with nonlinear materials, it is possible to introduce specific types of nonlinearities with respect to currents that arise in superconductors, as well as admitting simpler ways of entering nonlinearities with respect to the magnetic field.
• The -field formulation admits different ways of exciting the system. Depending upon your modeling needs, these can simplify the model setup.

It is this second point that we will examine in more detail here within the context of frequency-domain modeling. Most of the frequency-domain excitations that we will introduce will also work in the time domain, with a few exceptions that will be explicitly addressed.

Understanding the Common Boundary Conditions

The Magnetic Field Formulation interface, or the -field interface, includes four boundary conditions that are in common with the -field formulations and have the same interpretations:

• Magnetic Insulation boundary condition: This is primarily interpreted as a boundary to a perfectly conductive material. For a discussion about the various other ways in which this condition can be interpreted, see our Learning Center article: “Understanding the Magnetic Insulation Boundary Condition”.
• Perfect Magnetic Conductor boundary condition: This condition is primarily interpreted either as a boundary to a nonconductive media or as a type of symmetry condition. For more details, see our blog post: “Exploiting Symmetry to Simplify Magnetic Field Modeling”.
• Periodic boundary condition: This condition imposes periodicity, such as for structures with a twist (like submarine cables).
• Impedance boundary condition: This condition models the presence of a lossy material adjacent to the modeled space. It is a frequency-domain condition and should only be used when the skin depth is much smaller than the part size. For more details, see section 4 of our Learning Center article: “Understanding the Options for Meshing and Geometric Modeling of Conductors in Electromagnetic Fields”.

These conditions do not, on their own, introduce excitations, but they can form part of the solenoidal current path that is necessary to set up a solvable model. That is, even though they do not excite a current, they can define a path along which current can flow. 

Understanding the Excitation Boundary Conditions Specific to the MFH Interface

The three boundary conditions unique to the MFH interface that can be used to introduce excitations are:

• Tangential Magnetic Field boundary condition: This imposes the nonhomogeneous Dirichlet condition.
• Tangential Electric Field boundary condition: This introduces a nonhomogeneous Neumann condition.
• Surface Magnetic Current Density boundary condition: This is also a nonhomogeneous Neumann condition but can introduce a jump in the electric field across an interior boundary within the model space.

As their names imply, these conditions involve the tangential components of the fields on the boundaries where they are applied. If you were to enter a field that is purely normal to the boundary, that would be the same as a zero contribution. That is, the Tangential Magnetic Field boundary condition with a zero tangential component would be equivalent to the Perfect Magnetic Conductor boundary condition, and the Tangential Electric Field boundary condition without a tangential component would reduce to the Magnetic Insulation boundary condition. These conditions can only excite the system if either the electric or magnetic field has a component tangential to the boundary.

To help build up our understanding of these excitations, we will be using a model of a coaxial cable. It is worth remarking that despite its apparent simplicity, the coaxial cable is probably one of the best ways of learning electromagnetics, especially since there exist analytic solutions for the electric field and magnetic field variation with radius:

where the electric field points radially outwards, and the magnetic field circulates about the center. We’ll start by using these expressions to help excite our model and work forward from there.

Introducing an Electric Potential Difference with the Tangential Electric Field

We’ll begin by imposing an electric potential difference between the inner and outer conductors at one end of the model and short the other end. The short is imposed via the Magnetic Insulation boundary condition at the far end, and the insulating space around the exterior of the coax is modeled via the Perfect Magnetic Conductor boundary condition. To impose the electric potential difference, the Tangential Electric Field boundary condition is applied over the cross-sectional faces at one end, with the electric field being set to zero within the conductors, and the above equation is used to define the field over the annular face of the dielectric. 

The results are shown below. Currents flow through the inner and outer conductors, and at the operating frequency, the skin depth is apparent in the plot of the losses within the volume. The arrow plots of the electric and magnetic fields and the power flow match what we should expect from a shorted section of a coaxial cable. Power flow is primarily along the axis of the cable and dissipates into heat uniformly along the length. The short at the far ends means that there are zero electric fields. 

Plot of the heating (left) and the electric field (red), magnetic field (green), and power flow (blue) in a coaxial cable excited with an electric field at one end and shorted at the other (right).

Note that although the losses within the dielectric are small, they are nonzero. This is because the -field formulation does require a nonzero conductivity in all domains, unless solving in the frequency domain at very high frequencies, where the displacement currents are significant. This conductivity is often referred to as a numerical stabilization term. Using a direct solver can make this stabilization conductivity quite small. We use here an expression for the electric conductivity such that the skin depth in the insulators is ~1000 times greater than the domain size. Furthermore, we do not need to use gauge fixing for models that include numerically significant inductive effects, which describes most time- and frequency-domain cases.

It is important to note that, even though a coaxial cable is a type of transmission line, there is no transmission of power through this model. The fixed electric field injects power, and this is entirely dissipated into heat within the modeled space. There is no transmission or reflection. To introduce that, we will move on to the next boundary condition type.

Using the Surface Magnetic Current Density Boundary Condition

To introduce a boundary to the nonmodeled space on either end of the coaxial cable, we need to introduce the concept of impedance — the relationship between the electric and magnetic fields. Within the coax, we know the impedance is:

and we can use this to define the surface magnetic current density as:

where is the tangential component of the -field along the boundary. This condition can be used to model an impedance on either end of the cable. We can take this one step further and introduce a source term as well:

The factor of two arises because the source term will excite a signal propagating into the modeling domain, as well as a signal propagating into the boundary. This factor of two assumes that the impedance of the modeled space matches the impedance of the transmission line represented by the boundary, so any mismatch of impedance will lead to some reflection. This behavior is somewhat similar to the Cable-type Lumped Port boundary condition. However, the Lumped Port boundary condition in the Magnetic Fields interface excites via a surface current, while the MFH interface excites via an electric field, so the two are not exactly analogous. Note that the Surface Magnetic Current Density boundary condition acts in contribution with other Neumann conditions, so a separate feature can be used to contribute the source term to the matched impedance term.

Plot of the heating (left) and the electric and magnetic fields and power flow (right) in a coaxial cable excited with a matched impedance at both ends and an excitation at one end.

From the above plots, we can see that the power is flowing through the cable, so we now have a model of a two-port device. However, it is important to remark that only the frequency-domain form of the -field formulation includes the displacement current contribution. The time-domain form considers solely the inductive and conductive currents, so it is not a full-wave formulation but rather a quasi-static formulation. If you do want to model transient wave phenomena, use the Electromagnetic Waves, Transient interface, an -field formulation, as covered in the Learning Center article: “Modeling Capacitive Discharge”.

Sometimes, though, we are not interested in power flow but only the losses within the conductors. In a coaxial cable, for example, most of the losses are in the inner conductor, so we might want to analyze just that conductor and ignore the losses in the outer conductor. We also often want to assume a fixed current. For this, we use the last boundary condition, the Tangential Magnetic Field boundary condition.

Using the Tangential Magnetic Field to Impose a Current

Looking back at our analytic expressions for the magnetic field, we see that we know exactly what it will be between the inner and outer conductor, and that this expression is for the tangential component to a coaxial cylinder. With this, we can reduce our model to a cylindrical space around the inner conductor and apply the Tangential Magnetic Field boundary condition along its exterior. Conceptually, this is equivalent to imposing a tangential surface current density, so we use it in conjunction with the Magnetic Insulation boundary condition at the top and bottom to provide a solenoidal current path through the inner conductor. By imposing a magnetic field along the outside, we induce current to flow within the inner conductor. We do need to be a bit careful with this, though, as it has an interesting consequence.

Plot of the surface currents and heating (left) and the electric and magnetic fields and power flow (right) in a coaxial cable excited with a specified tangential magnetic field along the outside. The mesh on the excitation boundary is aligned with the direction of the current flow implied by the imposed tangential magnetic field.

The magnetic field that we are applying is exact, but the finite element mesh on which it is being applied is an imperfect description of a cylinder. There can always be some geometric discretization error that leads to the modeled geometry having sharp corners and edges in the finite element space. These lead to local singularities that can lead to what looks like quite noisy fields near the boundary. We can avoid this issue by ensuring that the mesh is aligned with the direction of current flow, as shown in the image above.

There is one other feature of this method that is worth remarking on. Observe from the plot that the electric field is now aligned with the axis of the cable, rather than going between the inner and outer conductors, and that the power flow is pointing inwards towards the center. This represents a completely different kind of excitation, and it could be argued that such an excitation is not physically realizable. In terms of the current distribution, magnetic fields, and heating within the conduction, it is still a useful modeling approach.

At this point, we’ve looked at three different ways of exciting a structure. To briefly summarize:

1. The Tangential Electric Field boundary condition can be used to impose a potential difference between the conductors of a device.
2. The Surface Magnetic Current Density boundary condition can model a boundary with specified impedance along with an electric potential source term.
3. The Tangential Magnetic Field boundary condition can be used to impose a known current.

Although we’ve looked here at a simple geometry, exactly the same techniques can be applied on more complicated parts, as long as you keep in mind that:

• The electric field should be specified at the cross section of a TEM-type transmission line or some portion of the model where it is reasonable to impose the assumption of curl-free fields.
• Displacement currents are being considered in the frequency domain but not in the time domain.
• The magnetic field excitation boundary has a different physical interpretation but an equivalent result in terms of the heating and requires careful meshing.

Let’s now introduce two other ways of exciting such models that are quite unique to the -field formulation.

Enforcing a Constraint on Current Within a Closed Path

One of the core capabilities of the modeling architecture of COMSOL Multiphysics^® is the ability to automatically interpret the symbolic expressions that are entered within the user interface (UI). This includes not just symbolic differentiation; it also extends to constraint elimination. Although a complete description of how this works is a bit beyond this short article, we will introduce how to implement this within the UI.

Returning to our previous model of only the inner conductor of the coaxial cable and the surrounding dielectric, we will now introduce an edge midway along the inner conductor. We want to impose the condition that a specified total current flows along the conductor, and from Ampere’s law, we know that the line integral of the magnetic field along a closed path equals the current flowing through it:

or, for nonmagnetic materials:

We can use an integration coupling to compute this current, and then use the Global Constraint feature to impose the condition that this integral must equal the desired current. From there, COMSOL^® takes care of everything else, and we can go ahead and look at the solution.

Screenshot of the Global Constraint feature that enforces current flow through an integration loop over a set of edges within the Magnetic Field Formulation interface.

Plot of the heating (left) and the electric and magnetic fields and power flow (right) in a coaxial cable excited with a global constraint on the current through a loop.

Observe that the current distribution and heating profile match the previous model. The magnetic field also matches the analytic expectation. The electric field and power flow, however, should again give us pause. The power flow implies that the edge is acting as a kind of an interior line source, but this again represents something that is not physically realizable, so the electric field does not have meaning. Despite that, the current distribution and the heating from such a model can be used. Note also that the heating distribution is smooth around the excitation edge.

The advantage of this approach is that we can excite current flow through any kind of structure. We do not need to impose any excitation at the boundaries of the modeling space. That is, we can excite a conductive loop entirely within the model space or different currents in different conductors.

Coupling the Magnetic Field Formulation to the Magnetic Fields, No Currents Interface

The last type of excitation that we will cover involves using another physics interface, the Magnetic Fields, No Currents interface. This interface solves for the magnetic scalar potential, the -field, and, as its name implies, is classically meant for the modeling of systems where there are no currents present.

It turns out, though, that when used in combination with the -field formulation, the -field interface can introduce an excitation using the Magnetic Scalar Potential Discontinuity boundary condition. This condition has to be applied to a boundary representing the cross section of a closed loop of current. The current can flow either through the volume of a conductive domain modeled using the -field or implicitly along the Magnetic Insulation boundaries of the -field interface. 

Schematic of a coaxial cable excited using the – approach. The Magnetic Scalar Potential Discontinuity boundary (green) excites current to flow along a solenoidal path through the two conductive domains (grey), modeled using the -field formulation, and along the Magnetic Insulation boundaries (yellow) at the top and bottom of the modeling domain. The arrows provide a visualization of this loop around the excitation boundary.

The – formulation will have less degrees of freedom (DOFs) than the -field formulation but can have greater computational costs since the resultant system matrix is nonsymmetric and usually requires a direct solver. It is conceptually, and sometimes geometrically, complex. It has the advantage that free space regions can be modeled as having zero loss, but there is no capacitive coupling though these domains. In fact, only the magnetic field is computed in the nonconductive domains, so the electric field and, hence, power flow are not defined. Only the magnetic fields, currents, and losses within the conductor can be extracted.

Plot of the heating (left) and the magnetic field (right) in a coaxial cable modeled using the – formulation.

A Quick Overview of Other Possibilities

There is a near-infinite complexity to the types of low-frequency inductive modeling that can be addressed using these excitation methods of the -field formulation. It should be emphasized that just because they can does not necessarily mean that they should be done using these approaches. It is very much something that needs to be determined on a case-by-case basis. With that being said, here are a few quick examples of other ways in which these approaches can be used.

Modeling a Unit Cell of a Twisted Cable

The Periodic Condition boundary condition can be imposed on faces that are identical but with a twist relative to each other. This boundary condition is demonstrated in the example of a twisted submarine cable and can be combined with the approach of exciting the current via a global constraint. 

A unit cell model of a one-quarter section of a twisted three-phase cable. The results on the modeled section can be patterned along the length, showing the heating and magnetic field.

Modeling a Symmetric Coil

When exploiting symmetry, it is often impractical to use the Lumped Port excitation within the Magnetic Fields interface; a Coil domain condition usually has to be used. This has a drawback at high frequencies, since the skin depth of the coil will need to be meshed very carefully. Using the MFH interface instead, the driving coil can instead be modeled using the Impedance Boundary Condition and the current imposed via the method of a global constraint.

A symmetric model of an inductively heated billet of material, using the -field formulation (left) and the -field formulation (right). The deposited heat is very similar.

Modeling a Solid Conductive Part within a Background AC Field

When a conductive part is exposed to an AC magnetic field, there will be induced currents within the volume of the material. As long as there are no holes in this part, which can intercept some fraction of the field, it is possible to model this using the coupled – formulation. It is strictly required that the conductive material does not form a closed loop.

When using the – formulation, the resultant system matrix is nonsymmetric, so computational costs will grow more rapidly than when using the -field formulation in all domains or the -field formulation, despite the lower number of DOFs. It thus is primarily useful if the ratio of the part volume to the entire modeling volume is very small.

A spiral piece of wire exposed to the background AC magnetic field. Note the wire does not form a closed loop. The -field is solved within the wire and excited by the boundary conditions applied to the -field solved for within the surrounding space.

Closing Remarks

We have shown here a few ways in which to excite the -field formulation to solve inductive problems, and of course there are many more variations on these techniques. Although these examples all are solved in the frequency domain, they are also applicable in the time domain with the caveats noted earlier. With these examples in our toolbelt, we’re ready to further explore the electromagnetic modeling possibilities with the AC/DC Module.

Note also that many of these examples are solved using direct solvers. It is worth remarking that with the NVIDIA CUDA^® direct sparse solver (NVIDIA cuDSS), available in COMSOL^® version 6.4, one can address larger and larger problems in less time than ever before.

To gain hands-on experience with the models discussed in this blog post, click the button below.

NVIDIA and CUDA are trademarks and/or registered trademarks of NVIDIA Corporation in the U.S. and/or other countries.

This transition is especially visible in enhanced geothermal systems (EGS). The question is no longer whether engineered geothermal reservoirs can work, but how to optimize them for scalable, long-term performance.

From Niche Resource to Multiuse Subsurface Platform

The growth of the workshop mirrors a broader global trend. Geothermal energy is increasingly being recognized for its potential to provide reliable baseload power with a relatively small land footprint compared to solar and wind. But the real transformation lies in the shift from conventional hydrothermal systems, limited to specific geologic settings, toward “geothermal anywhere”.

EGS enable this shift by creating engineered reservoirs in hot, dry rock deep underground. By hydraulically fracturing these rocks and circulating water through the system, EGS turns previously inaccessible heat into usable energy.

Conceptual evolution of geothermal technologies from conventional hydrothermal to next-generation EGS and superhot systems. Image source: U.S. Department of Energy.

At the workshop, keynote talks highlighted that the industry is moving beyond EGS proof-of-concept toward efficient execution by refining well placement, stimulation design, and multiwell connectivity, improving performance and repeatability of this method.

Furthermore, geothermal is evolving beyond standalone power generation. Several sessions at the event explored opportunities for critical minerals and lithium coproduction, subsurface thermal energy storage, and geological hydrogen pathways. Increasingly, geothermal sites are turning into multiresource platforms.

Superhot Rock: The Next Frontier

One of the most ambitious directions discussed was superhot rock (SHR) — the industry’s moonshot, targeting temperatures above 400°C at depths greater than 10 kilometers. Under these conditions, water enters a supercritical phase, carrying significantly more energy than conventional steam.

Under optimal conditions, a single supercritical well could produce vastly more power than today’s standard geothermal well, but accessing that potential introduces new engineering challenges the industry is racing to address: Materials behave differently at extreme temperatures, well integrity becomes more complex, and brittle rock shifts to a ductile state under higher pressure and heat.

Recent field progress suggests that this frontier is closer than many anticipated. At the workshop, Mazama Energy reported reaching 331°C in Oregon. Fervo Energy demonstrated rapid drilling performance in a 290°C sedimentary basin well using an AI-enabled drilling system, showing that EGS is not confined to crystalline rock. Quaise Energy shared progress on millimeter-wave drilling technology aimed at reaching greater depths.

As geothermal pushes into extreme environments, tighter integration across disciplines becomes essential — a recurring theme throughout the workshop.

Conceptual geothermal system showing heat extraction from a geothermal reservoir (< 250^°C) and the deeper superhot rock zone (> 400^°C). Conventional systems typically operate at shallower depths, while superhot resources target depths of several kilometers (> 5 km), where higher temperatures increase energy density and can enable up to ~10× greater power output per well compared to conventional geothermal systems.

The Future of Geothermal Depends on Coupled Systems

In the discussion of moving the industry forward, something the presenters emphasized was optimizing the “3 Cs”:

1. Connectivity (well communication)
2. Conductivity (heat uptake)
3. Conformance (uniform flow)

Achieving optimization depends on understanding how these factors evolve together over time.

Today’s EGS projects already involve fracture creation, heat extraction, fluid flow, and chemical transport. Temperature influences stress. Stress alters permeability. Permeability reshapes flow paths and heat recovery, while chemical processes modify fracture performance. This interconnected behavior — thermal–hydraulic–mechanical–chemical (THMC) coupling — is central to long-term geothermal performance, and moving toward 400°C increases that interdependence.

The thermal–hydraulic–mechanical–chemical (THMC) framework for geothermal systems. Reservoir performance depends on the continuous interaction of heat transfer, fluid flow, mechanical deformation, and chemical processes.

In that sense, the geothermal industry increasingly requires a fourth C: coupled modeling, or the ability to evaluate interacting physical processes. Modeling multiphysics within a unified framework can reduce uncertainty and support scalable development.

Geothermal is not just getting hotter; it’s becoming more integrated and data-driven, with multiphysics thinking at the core of system design. The projects that succeed will connect diagnostics, optimization, and coupled physics into repeatable strategies.

Further Reading

Learn more about how multiphysics simulation can be used within the geothermal industry by checking out these blog posts:

• Modeling Geothermal Processes with COMSOL^® Software
• Coupling Heat Transfer with Subsurface Porous Media Flow
• Fully Coupled Hydromechanical Modeling of Fractured Media

What Is Sputtering?

Sputtering is a physical process where ionized particles of a plasma or gas are bombarded on a solid surface, or target, to displace the target’s atoms, which can then land on another surface, or substrate. Sputtering occurs naturally in outer space and can be an unwelcome source of wear in precision components. However, it is also useful for thin-film deposition and for physical surface removal in semiconductor manufacturing, optics, and MEMS applications. A schematic of a sputtering setup can be seen in the figure below.

A typical sputter deposition setup. Image in the public domain, via Wikimedia Commons.

Various types of sputtering mechanisms, namely magnetron, reactive, and ion-beam sputtering, are used in the semiconductor industry to deposit thin films, which are crucial for creating the complex multilayer structure required in today’s integrated circuits. These films may form barrier layers, adhesion layers, or conductive layers within chips. Sputtering offers several advantages in this regard. For instance, it enables excellent control over the thickness of the deposited film while being applicable for different materials, such as metals, alloys, and compounds. Moreover, sputtering is the ideal choice for substrates that cannot withstand higher temperatures.

In addition to thin-film deposition, sputtering is also used for the selective removal of material from a wafer in dry etching by employing ionized gas particles guided by a protective mask. Inert gases such as xenon and argon are typically preferred as sputtering gases because of their minimal reactivity and ability to cause higher displacement of the target material due to their significant molecular weight.

A sputtering device, with plasma in a vacuum glass tube.

Understanding Complexities with Simulation

While sputtering is commonly used in chip manufacturing and thin-film processing, complexities arise mainly due to various phenomena simultaneously at play, such as plasma generation and ion bombardment. Hence, maintaining meticulous balance and control over various electrical, chemical, and physical parameters, such as gas pressure, ionization potential, and substrate bias voltage, is critical to achieve the desired outcomes.

Numerical simulation can be used to analyze factors such as pressure, geometry, and material properties, making it beneficial in the optimization of thin-film growth, substrate erosion, and device fabrication. Furthermore, virtual analysis and optimization can help teams save on costly physical experiments.

Let’s look at an example of modeling argon ion sputtering on a silicon surface in the COMSOL Multiphysics^® software.

How Is Argon Sputtering Modeled in COMSOL Multiphysics^®?

The Particle Tracing Module, an add-on to COMSOL Multiphysics^®, provides physics interfaces that can be used in sputtering-related models. In this example, the Charged Particle Tracing interface is used to simulate argon ion trajectories near a silicon surface. It can be used to compute ion and electron trajectories in electric and magnetic fields. The Deformed Geometry interface is used to visualize the surface evolution caused by sputtering and capture the resulting morphological progression.

Representation of the model setup.

As seen in the image above, the sputtering process involves two ion beams incident at 45^° and −45^°, respectively. When argon ions hit the silicon surface, silicon atoms are sputtered. Each particle is assumed to represent 10¹² actual argon ions. The mask height is taken as 0.2 μm, and the etching window width is 0.5 μm. Mask sputtering is ignored in this simplified example.

What the Results Depict

The figure below shows the sputtering yield in relation to the incident angle of the beam. Sputtering yield refers to the number of target atoms sputtered per incident ion. As can be seen, the number of secondary particles varies with respect to the angle of incidence. In this simplified model, the direction of secondary particle emission is assumed to follow a cosine distribution, and elastic-collision-based assumptions are used to estimate the emitted particle velocity. As observed in the graph, sputtering yield increases almost linearly with respect to incident angle between 30^° and 70^°. Accordingly, the sputtering yield can be calibrated by controlling the incident angle of the beam. For example, a lower angle of incidence can be employed for lesser material removal, whereas increasing the angle of the incident beam will result in more removal of target material within this range.

Number of secondary sputter particles with respect to the incident angle.

Particle trajectories and sputtering progression with beams incident at 45^° and −45^°.

The above figure depicts particle trajectories with beams incident at 45^° and −45^° and the resulting sputtering progression. More material is observed to be sputtered at the corners, resulting in faster height reduction.

How Can These Findings Be Improved?

This model can be further extended to incorporate other factors that affect sputtering. For example, in addition to the incident angle, sputtering yield also depends on the energy of the incident ion as well as the absolute and relative mass of the incident ion with respect to the target atoms. Furthermore, the bombardment process is significantly more complex in the real world, with some ions penetrating the substrate to generate cascading collisions and atoms being ejected only after multiple scattering events. These effects are not captured in this simplified model and generally require more specialized sputtering or collision-cascade descriptions.

Try It Yourself

Want to try modeling the sputtering process yourself? Download the related MPH-file in the Application Gallery:

Inducing Streaming with Surface Acoustic Waves

When alternating voltage is applied on the surface of a piezoelectric material, it generates a strain determined by the electrical field, and waves start to propagate on the surface. These waves are called surface acoustic waves (SAWs), and they are distinguishable by how the material deforms relative to the propagation and normal directions. Two types of SAWs include Rayleigh waves and Love waves. This blog post will focus on Rayleigh waves, which make the surface deform in the normal direction. For the generation of the SAWs on the substrate, normally a set of comb-like terminals, or interdigital transducers (IDTs), is used to impose alternating electric potential. IDTs can both generate and receive the SAWs. When used as a filter in electric components, one more set of terminals is placed on the path of generated SAWs. The two terminals of the IDT on the receiver side will have different electric potentials according to the strain the substrate is experiencing, and information about the oscillation can be determined.

Instead of placing the second IDT on the surface as a receiver, we will put a droplet on the propagation path. The droplet will start interacting with the SAWs and absorb their energy. The SAWs attenuate as they travel under the droplet and are called “leaky SAWs” for this behavior. Energy is radiated to the droplet in the form of bulk waves incident at an angle called the Rayleigh angle. In the droplet, the incident waves reflect on the free surface of the droplet while also losing energy due to viscous dissipation, eventually resulting in a steady circulating flow component called acoustic streaming. We can induce steady flow just by oscillations. This plays an important role in the microfluidics area — we can enhance mixing inside the droplet without needing to physically put something in the fluid to stir it; the method is noninvasive. The resulting streaming may have different circulation patterns inside the droplet depending on the energy of the waves, the dimensions of the system, the material properties of the droplet, and so on.

Schematic of SAW-induced streaming. The SAWs (wavy blue line) depart from the IDT. Once they reach the droplet, the energy is transferred into the droplet (solid yellow arrow), and finally, the streaming occurs (dashed yellow arrow).

As we have seen so far, the droplet streaming has multiple physics areas involved. Due to their complexity, these contributing factors are frequently modeled separately by dividing the process into several steps. Note that the analysis of the acoustic field in a droplet takes up considerable RAM (under the current conditions), while the analysis of the other fields usually does not. We will use COMSOL Multiphysics^® to handle the complexity and create a model covering the whole process of transferring energy in the system. In this model, the drop has a wetting diameter of about 2 mm and a contact angle of 78° with the solid surface. The SAW device is excited with a frequency of 20.37 MHz. The droplet is assumed to be a glycerol–water mixture, and the properties of the piezo crystal are taken from the lithium niobate material in the Piezoelectric material library. As a reference, the values and setup are similar to the ones used in the paper: “On the Influence of Viscosity and Caustics on Acoustic Streaming in Sessile Droplets: An Experimental and a Numerical Study with a Cost-Effective Method” (Ref. 1).

Streaming Model Setup in the 2D Configuration

First, we will check what the model looks like in 2D. We know that the streaming has a 3D structure because of the hemispherical shape of a droplet, but it is always a good starting point to create a 2D model to see if the necessary settings and physics are complete and that the phenomenon we hope to model occurs under the 2D assumption. We will use the same cut plane as the schematic shown above for the 2D simulation. We will use the Electrostatics, Solid Mechanics, Pressure Acoustics, and Creeping Flow interfaces to excite the SAWs, and to capture the streaming, the Acoustic Streaming Domain Coupling multiphysics is used. Given that the time scales are quite different between the oscillation and the streaming, the simulation is done in two steps: a Frequency Domain study and a Stationary study.

In piezoelectric analysis, we need to pay attention to the crystal cut of the piezoelectric material. In this model, 128^° YX-cut lithium niobate (LiNbO₃) is used for the substrate; thus, the rotation should be reflected on the material properties. Functionality in COMSOL Multiphysics^® enables us to consider the angle of the crystal cut by defining a coordinate system using, for example, the Rotated System feature and specifying the coordinate system in the Piezoelectric Material node of the Solid Mechanics interface. We also have an entry in the Application Gallery that explains the coordinate settings: Euler Angles Rotation in SAW Modeling. Note how the angles are set differently in the Application Gallery model for a 2D component (XY sagittal plane) vs. a 3D component (XZ sagittal plane). Further below, we will create a 3D model with the XZ plane as the sagittal plane.

Settings of the Rotated System feature in the 2D model. In the 3D model, the value for β is set to -38 [deg].

We also need to make sure that an appropriate loss mechanism is taken into account in the Pressure Acoustics node. In the current setup, it is mostly the Eckart streaming that drives the flow inside the droplet. Therefore, the bulk attenuation of the sound wave should be modeled to capture it. The Fluid model in the Pressure Acoustics node specifies what type of attenuation the acoustic waves will experience. Here, we simply choose Viscous. If we left Linear elastic selected (the default), we would see no streaming in the results.

One last point to check is the Values of variables not solved for settings in the Stationary step. In the streaming analysis, the Acoustic Streaming Domain Coupling feature couples the Frequency Domain and Stationary studies. When the coupling is activated in a Stationary study, it refers to the variables solved in a Frequency Domain study to calculate the terms contributing to the streaming. Since we use multiple study nodes, the coupling does not know which solution has the frequency domain data of interest, so we need to specify it in the settings of the study step.

Now, let’s run the Frequency Domain and Stationary studies in order. The result should show distributions like the images below. The IDT is placed to the left of the droplet, outside of the images. There are SAWs generated on the piezoelectric material, moving to the right. The wave propagation direction can be checked more clearly using the Animation feature with the Dynamic data extension. The surface waves turn almost invisible after they have traveled over half of the contacting area. In exchange, the bulk waves in the droplet propagate in the upper-right direction, resulting in a complex pressure pattern. As might be expected, the circulating flow field is also confirmed from the stationary result. It has a large vortex across the whole domain area, but note that this might be 2D-specific. In a 2D configuration, we cannot simulate a vortex whose axis is not normal to the screen. However, it is a good start that we have set up a SAW model that induces streaming similar to what we have expected. Let’s move on to the 3D model.

Results computed in the 2D model. Stress in the substrate and acoustic pressure in the droplet (left); displacement of the substrate and velocity distribution of the streaming in the droplet (right).

Tackling the RAM Issue in the 3D Configuration

Despite the difference in the dimension, from 2D to 3D, the principles of the simulation do not change. We use the same interfaces with the same multiphysics couplings. What does differ from the 2D model is the requirement for the computational resources. The 3D model’s geometry dimension is much larger than the wavelength, and it is highly likely to face memory issues. Not only would we experience longer computation, but we would also need to reduce the memory consumption so that the model can fit into the RAM of the computer.

First, we can simplify the geometry to an extent that will not significantly affect the result. In the 3D model, the IDT is modeled as several sets of simple rectangular 2D terminals that are placed in parallel. Moreover, the droplet and the substrate are cut in half by the middle plane to halve the degrees of freedom (DOFs). A periodic condition is used for the transverse direction of the piezo while a symmetry boundary condition is applied to the drop. This is a good approximation for the physics, especially when we are interested in the flow field in the droplet. We would need to widen the substrate domain if the model exhibited a strong dependency on the transversal length.

Second, considering that the resolution levels required to capture the waves are different in the droplet than in the substrate, we will use different mesh and mesh sizes for them. This is achieved by using the Form Assembly method instead of the Form Union method in the geometry sequence. This feature will allow the model to have multiple geometry objects in a component, mesh each object separately, and connect them using the so-called Pair feature that couples nonconforming meshes. Remember to use the Union operation where applicable so that the software can recognize that some of the geometry instances belong to a single object. The domains that belong to the same object will have a conformed mesh, and no pair is used between them.

Magnified view of the nonconforming mesh near the leading edge of the droplet.

Lastly, we need to use an iterative solver to reduce the RAM usage in the Frequency Domain study. A general guideline on the settings of iterative solvers for acoustics problems is available in the Solving Large Acoustics Problems Using Iterative Solvers section of the Acoustic Module’s User Guide documentation. In this model, the number of DOFs for the Electrostatics and Solid Mechanics interfaces is much smaller than for the Pressure Acoustics interface, and direct solvers are expected to work well for solving for the electric potential and the solid displacement fields. Therefore, we activate hybrid preconditioners and use a Direct Preconditioner for the dependent variables of the Electrostatics and Solid Mechanics interfaces. With this feature, we can use a direct solver for small fields while applying efficient solvers to the other large fields. If the model did not have a piezoelectric domain, we could use a Segregated solver to reduce the RAM usage even further. However, as explained in Solving Large Acoustic–Structure Interaction Models section of the Acoustic Module’s User Guide documentation, acoustic–structure interaction models with piezoelectricity need to use the Fully Coupled solver, so the linear solver is the only part we can tweak. The acoustics part of the equations uses the shifted Laplace approach as an efficient form of the multigrid preconditioner.

The Settings window for Direct Preconditioner in the 3D model. Note that hybridization is activated by choosing Multi preconditioner in the Hybridization section.

Now, the Frequency Domain study solves with 130 GB RAM in about 1 hour 10 minutes in our environment. However, even after running the study, there is something we need to keep in mind: rendering. In a large model, rendering a result can take a long time. To simplify working with large models, it is recommended that you enable the Only plot when requested checkbox in the Settings window for the Results node.

The results below clearly show the influence of the 3D geometry. Now the pressure pattern has a more defined peak in the middle slice than the 2D result. In addition, although it is a bit difficult to discern from a single picture, a small vortex is formed near the leading edge, and the rest of the droplet is occupied by a large vortex. These show good agreement with the reference, where a 3D simulation was conducted by decomposing the computation into 2D subproblems. In our model, we have been taking a rather straightforward approach; thus, we can simply set up multiphysics modeling with some modifications in the settings. The strategy to use the nonconforming mesh and the iterative solver would also work for other large problems, such as complicated MEMS devices.

Results computed in the 3D model. Stress on the substrate surface and acoustic pressure in the droplet (left); streamlines of the streaming colored by the velocity magnitude (right).

Now It’s Your Turn

In this blog post, we have dealt with the complexity of modeling streaming in a droplet and the multiphysics involved using 2D and 3D models. We observed the interactions between physics interfaces, even with the 2D model, that solved quickly. Extending a 2D model to a 3D one can sometimes be challenging; it might require trial and error to see what configuration works for each model. We hope the current strategy will be useful to your multiphysics problem. The models are available from the following links:

• SAW-Induced Streaming in a Droplet – 2D Setup
• SAW-Induced Streaming in a Droplet – 3D Setup

As mentioned above, piezoelectric materials are frequently used to excite SAWs, but the definition of the properties requires special attention in terms of the coordinate system. Moreover, the design of the IDTs is closely related to the wave speed to excite; therefore, it is also important to see if SAWs are generated as desired. The following models work as useful references to test your SAW setup before building a complex model:

• Euler Angle Rotation in Surface Acoustic Wave Modeling
• Surface Acoustic Wave Velocity Calculations from a Unit Cell

This blog post does not mention particle movement in the streaming flow. To model acoustic manipulation of such things as particle movement, it might also be necessary to consider acoustic radiation force exerted on the particles. The Particle Tracing for Fluid Flow interface includes this functionality in the Acoustophoretic Radiation Force node, which can be used along with the Drag Force node. This combined use will consider the force due to the steady flow component at the same time. The following models would be useful as the starting point of such applications:

• Acoustic Streaming in a Microchannel Cross Section
• 3D Acoustic Trap and Thermoacoustic Streaming in a Glass Capillary

Reference

1. A. Riaud et al., “On the Influence of Viscosity and Caustics on Acoustic Streaming in Sessile Droplets: An Experimental and a Numerical Study with a Cost-Effective Method,” Journal of Fluid Mechanics, vol. 821, pp. 384–420, 2017. DOI: https://doi.org/10.1017/jfm.2017.178

Cupping Therapy Explored

Suction pressure elevates skin, fat, and muscle to varying extents depending on the pressure level and contact area. Along with cupping therapy (Figure 1), suction is used in medical applications such as cryolipolysis for fat freezing and reduction, lymphatic therapy to promote fluid drainage, and laser therapy, such as in port-wine stain (PWS) laser treatment, where lasers are used to lower blood vessel visibility.

Etim and her team are working to increase understanding of how suction pressure and applicator size influence the stress distribution across tissue layers. Computational modeling can be used to study the stress patterns and specific responses in each layer of skin during suction. By understanding how skin deforms under suction pressure, researchers can better predict the effects of devices and assess extreme loading conditions. Moreover, these insights can help medical device designers develop more reliable models and optimize device design.

Before diving into how Etim and her team used simulation, let’s consider the layers of skin and their relationship with medical applicators.

Figure 1. Cupping therapy, which results in deformation of the skin and underlying tissue due to suction.

Layers of Skin & Aperture Sizes

Skin is made up of three primary layers:

1. Epidermis
2. Dermis
3. Hypodermis

The epidermis is the protective outer layer of skin, consisting of four layers over most of the body and five layers on the palms and soles. The dermis contains blood vessels, nerves, hair follicles, and sweat glands. This layer deforms anisotropically due to the preferential orientation of collagen fibers. The foam-like hypodermis insulates the body, cushions organs, and stores subcutaneous fat, behaving as an energy-absorbing layer because of its collagen matrix, which encloses lipid-filled adipocytes.

Suction devices have typically been made with apertures under 10 mm; however, clinical applications such as cupping therapy, laser-assisted PWS removal, and some cryolipolysis applicators can use larger apertures ranging from 30–65 mm, highlighting the need to evaluate a broader range of loading geometries. Smaller apertures primarily engage the epidermis and dermis layers, where the fiber alignment may play a stronger role, whereas larger apertures engage the hypodermis and muscle underneath.

Now, we’ll dive into how simulation was used to analyze tissue.

Setting Up the Model

Etim and her team used a phenomenological approach to estimate the properties of skin, fat, and muscle for different participants, also considering different aperture sizes. They began with a single-phase solid mechanics model in COMSOL Multiphysics^®, with the governing momentum balance equation

where ρ is the density, u is the displacement field, f_v represents the volume forces, F is the deformation gradient, and S is the second Piola–Kirchhoff stress tensor. The Green–Lagrange strain tensor is defined as

Stress and strain are connected by the strain energy density function W as

Skin () was modeled using the polynomial model:

and fat () was modeled using the Mooney–Rivlin formulation:

The muscle was modeled with an Ogden formulation:

The team implemented a multilayered, axisymmetric finite element model of skin, fat, and muscle in COMSOL Multiphysics^®. Skin thickness was set to 2 mm, muscle thickness was set to 10 mm, and fat thickness was set to values specific to the results of participant ultrasound measurements gathered in a study. The simulation represented suction loading for applicators with aperture sizes of 50, 30, and 16 mm. The layers are assumed to be bonded without friction. Etim explained that a uniform suction pressure was applied at the skin surface over each aperture radius, and a boundary condition was imposed to enable the skin to slide along the rim of the cup.

Etim and her team used the Optimization Module, an add-on to COMSOL Multiphysics^®, to estimate the properties of skin and fat while keeping the properties of muscle constant. was assumed to be zero, while , , and were optimized by minimizing the deviation between the deformations seen in experiments and those in the model. A time-dependent model involving a wide pressure range was used, but the material properties in this optimization were initialized with the results of an optimization for a single high pressure. This methodology enabled the team to capture nonlinear pressure-displacement behavior.

Simulation Results

The results of the simulation showed that tissue deformation increased with suction pressure, with the magnitude and distribution influenced by aperture size and fat thickness (Figure 2). As expected, the simulation showed larger deformation for both larger suction pressure and larger aperture size, with the smaller apertures increasing the engagement of the upper skin layers.

Since the results were based on real data from individual participants with varying levels of fat thickness, the results also showed that differences in fat thickness of approximately 2 mm between participants led to distinct shifts in deformation profiles.

Figure 2. Displacement of skin and fat at 27091.109 pascals (Pa) for a participant with a 4.0-mm fat thickness when aperture sizes of 50 mm (a), 30 mm (b), and 16 mm (c) were used.

The effect of different material property combinations was investigated, and it was found that variations in the fat stiffness of the participant had a stronger influence on deformation with the 50-mm aperture than the 30-mm and 16-mm apertures. These results highlight the importance of taking into account the applicator geometry as well as the composition of the tissue in the individual participant.

Let’s take a closer look at a stress analysis in one individual participant who has a 6.3-mm fat thickness. The skin reacted with the highest stresses near the center of the suction cup, particularly in the circumferential and axial directions when a pressure of 27091.109 Pa (or 203.2 mmHg) was applied. In contrast, the fat layer showed low circumferential and axial stresses, despite large deformation. Shear stress remained minimal across conditions. The muscle experienced less stress, although it was more widely spread out. Other participants who had thinner fat thickness resulted in higher and more widely distributed values. In the muscle layer, the largest stress component is the axial stress (Figure 3). This stress analysis suggested to Etim and the team that fat thickness modulates load transfer and acts as a protective layer, providing a buffer to deeper tissue during suction loading.

Figure 3. Three components of the stress tensor are plotted at the mid-depth of the skin, fat, and muscle layers. The results are for a pressure of 27091.109 Pa in a cup with a 50-mm aperture.

Insights from Simulating Tissue

The results of this study indicate that both aperture size and fat thickness significantly influence the deformation caused by suction. The larger aperture, 50 mm, engaged more of the fat layer, deeper into the tissue, and produced more consistent results. The smaller apertures, 30 mm and 16 mm, engaged primarily the upper skin layers and showed higher variability. As Etim explained, these results highlight the importance of characterizing tissue under multiple loading geometries and pressures. Accurate representation of both skin and fat properties is crucial in the design and performance of medical devices that use suction.

Further Learning

For more information on this work, read the UML team’s full paper, which won a Best Paper award at the COMSOL Conference 2025 Boston! The paper outlines the team’s modeling approach and results.

Creating Surrogate Models: The Workflow

How do you create a surrogate model in COMSOL^®? The workflow can be described in a few steps as follows, starting from a completed parametric single physics or multiphysics model:

1. Add and run a Surrogate Model Training study, which is based on design of experiments (DOE) methods to sample the model parameter space.
2. Add a suitable surrogate model and train it on the simulation data stored in a Design Data table. Optionally train the surrogate model on, for example, experimental data.
3. Use the surrogate model in an app or digital twin, or for other purposes.

This workflow is also illustrated in the figure below.

The workflow for using surrogate models.

Using Surrogate Models in Simulation Apps

One practical use of a surrogate model is to speed up a simulation app created with the Application Builder. In the example of the microstrip patch antenna app shown below, a function call to a surrogate model replaces the need to solve the full finite element model, resulting in near-instantaneous response times when varying the antenna’s dimensions or material properties. In the image below, we can see how the app displays the antenna gain pattern from the trained surrogate model, shown alongside the result from the original model for comparison.

A simulation app comparing results from a surrogate model with those from a full simulation model.

The image below shows another example of a simulation app that has been accelerated using surrogate models. In this case, a set of surrogate models is used to reconstruct the electric potential, temperature, and stress in a MEMS actuator. The user can interactively adjust four geometric dimensions of the CAD model, along with the applied voltage, using sliders. Thanks to the surrogate models, the app responds quickly, enabling a much more interactive experience than would be possible with a full simulation model.

A simulation app of a MEMS actuator that uses multiple surrogate models to quickly visualize and evaluate physical field quantities such as the temperature, stress, and electric field.

In the MEMS app shown above, the visualizations and result evaluations are generated from a set of surrogate model functions that are called upon behind the scenes. A function corresponding to the temperature field is shown in the image below.

A function call to a deep neural network (DNN) surrogate model function, used behind the scenes in the MEMS app.

The syntax dnn1_T(x, y, z, dw, gap, wv, L, DV) calls a DNN function named dnn1_T, with the eight input arguments listed in parentheses:

• Three spatial coordinates: x, y, and z
• Four CAD dimensions: dw, gap, wv, and L
• The applied voltage: DV

This type of function call replaces calls to the field quantities defined by the full simulation model, which in this case is a finite element model.

Each surrogate model can define multiple functions, each with any number of input arguments. These functions typically represent physical quantities such as an electric field, temperature, or stress. Since surrogate model functions can be differentiated, they are well suited for use in gradient-based optimization workflows, such as inverse modeling, where sensitivities with respect to input parameters are required.

Types of Surrogate Models

There are three types of surrogate models available in COMSOL^®: DNN, Gaussian process (GP), and polynomial chaos expansion (PCE). The DNN surrogate model is included in the platform product and does not require any add-on products. The GP and PCE surrogate models are part of the Uncertainty Quantification Module, where they are automatically created and trained using dedicated solvers or studies for uncertainty quantification. However, any of the three surrogate model types can be trained on any kind of simulation or experimental data.

A response surface for a GP surrogate model showing the uncertainty estimate (standard deviation).

Once trained, surrogate models are available as functions under the Global Definitions node, ready to be used throughout the model. Each type of surrogate model comes with its own advantages. The choice of surrogate model depends on the problem at hand: DNNs are powerful for complex, high-dimensional problems with large training sets, while GP and PCE models are better suited when you need access to the confidence or uncertainty in a prediction.

A DNN surrogate model definition that includes six functions in eight input arguments. The user interface makes it possible to customize the number of layers as well as the number of nodes per layer.

Note that for small datasets, GP models may be easier to create, and they perform better than DNN models.

Now let’s take a closer look at how the data used to train these models is generated.

Generating the Training Data

The Surrogate Model Training study is used to generate the data needed to train surrogate models. It performs a parametric sweep using methods based on DOE, and it can be configured to sweep over virtually any combination of input and output parameters. The result is a table of simulation data that serves as the basis for training. Surrogate models are not limited by physics; they can be used in applications across electromagnetics, structural mechanics, acoustics, fluid flow, heat transfer, chemical engineering, or any combination of multiphysics.

The first few rows of a data file generated by the Surrogate Model Training study.

Training Surrogate Models

Training is the step where the collected data is fitted to a surrogate model. Once trained, the surrogate model can be used in place of the original simulation, achieving significant speedup while maintaining sufficient accuracy in many cases. Surrogate models can either be trained automatically after the data generation step or added and trained manually in a separate step.

The fidelity of a surrogate model is controlled by the amount and quality of training data. A higher-fidelity model generally requires more data, which can be obtained from simulations, physical experiments, or a combination of both.

When working with surrogate models, the data-generation step is typically more time-consuming than the training step. However, both steps can be accelerated. Data generation can be sped up by running simulations on a cluster, allowing multiple design points to be computed in parallel. Training DNN surrogate models can also be accelerated using GPUs, which can significantly reduce training time for large datasets.

An app for analyzing the test cycle of a battery, highlighting how a DNN surrogate model can be used to reconstruct time-varying physical quantities.

Dive Deeper into Surrogate Models

In this blog post, we have provided a brief introduction of the functionality in COMSOL^® for creating and building surrogate models. To get a more comprehensive overview of this functionality, check out our Learning Center course on surrogate modeling, which is an 8-part self-guided course that covers an introduction to creating surrogate models, fitting data with a DNN, evaluating model uncertainties, geometry sampling, and more.

Tip: To learn about the theoretical background of surrogate modeling, check out this Learning Center course: “Surrogate Modeling Theory”.

MM MaschinenMarkt “Design Engineering” Award: Best of Industry 2025

The Best of Industry Award, presented annually by MM MaschinenMarkt and determined by a public vote from their readers, celebrates outstanding achievements that move the development of industrial technologies forward. As a leading trade publication for mechanical engineering and manufacturing, the MM MaschinenMarkt community is comprised of engineering experts, designers, and business leaders who are up to date on the latest innovations, rendering their recognition of the Electric Discharge Module an even greater honor.

At left: Phillip Oberdorfer, technology communication manager at COMSOL (second from right), receiving the Best of Industry Award 2025 on behalf of the company, among various other award recipients. Photo courtesy of the Vogel Communications Group. At right: A close-up of the Best of Industry Award 2025 for the “Design Engineering” category.

MM MaschinenMarkt presents this award in association with trade media outlets such as Automobil Industrie, PROCESS, Labor Praxis, DeviceMed, and konstruktionspraxis, covering a comprehensive list of achievements in manufacturing and construction, process manufacturing, laboratory and medical technology, and automotive development.

Speakers presenting on stage at the MM MaschinenMarkt Best of Industry Awards 2025. Photo courtesy of the Vogel Communications Group.

We are grateful to be recognized alongside so many leading voices in technical advancement and are incredibly excited about the role that the Electric Discharge Module will continue to play in the development of next-generation electrical systems.

Understand, Analyze, and Predict Electric Discharges and Breakdown

The Electric Discharge Module extends the capabilities of COMSOL Multiphysics^® for assessing electrical insulation performance, adding specialized functionality for modeling electric discharges in gases, liquids, and solid dielectrics. Engineers and researchers can simulate phenomena such as streamer, corona, dielectric barrier, and arc discharges; lightning-induced electromagnetic pulses; surface charge accumulation; and more — all within a unified multiphysics environment. The Electric Discharge Module seamlessly integrates with other products in the COMSOL product suite, including those for electromagnetics, structural mechanics, and fluid dynamics, enabling users to explore the full spectrum of multiphysics effects often associated with electric discharges.

“The combination of specialized functionality and multiphysics integration reflects COMSOL’s long-standing focus on bridging the gap between physics understanding and practical simulation,” said Lipeng Liu, technical product manager for the Electric Discharge Module at COMSOL. “Built on the COMSOL Multiphysics^® platform, the Electric Discharge Module gives engineers and researchers robust, validated capabilities for modeling electrical insulation and discharge phenomena while helping reduce reliance on costly experiments.”

That focus is especially important in industries like electronics, aerospace & defense, automotive, healthcare, and manufacturing, with applications ranging from consumer electronics to high-voltage power system components. For many teams, the Electric Discharge Module serves as an important tool for product development, helping to reduce costs associated with experimental testing and prototyping while garnering a better understanding of electrical insulation design.

A streamer discharge propagating in transformer oil simulated with the Electric Discharge Module.

We would like to thank the MM MaschinenMarkt community and the other award recipients for continuously pushing technical boundaries in the spirit of developing better products and more efficient manufacturing solutions.

Next Steps

Learn more about the award-winning Electric Discharge Module add-on to COMSOL Multiphysics^® via the button below.

History of Estimating Pi

Historically, pi has been calculated through a variety of geometric methods by several ancient cultures. Later, during the 17^th and 18^th centuries, infinite series like the Gregory–Leibniz and Machin-like formulas advanced pi calculations, allowing mathematicians to compute up to hundreds of digits by hand.

Using modern technology, the calculation of pi has reached extraordinary precision using powerful computers and advanced algorithms. In a previous blog post, we discussed how the Monte Carlo method can also be used to estimate the value of pi. Today, state-of-the-art calculation methods involving rapidly converging infinite series, such as those derived from Srinivasa Ramanujan and the Chudnovsky brothers (Ref. 1), have allowed us to compute trillions of digits of pi.

The Simple Pendulum

In this blog post, however, we are moving away from sophisticated technology and will instead rely on an experiment you can easily replicate in real life. A simple pendulum consists of a mass (ideally point-like) called a bob, attached to one end of a string (ideally massless). The time period T of such a pendulum’s swing, which is the time it takes to complete one oscillation back and forth, can be formulated as

where L is the length of the pendulum and g is the acceleration due to gravity. If the time period of a pendulum’s swing is known, the value of pi can be estimated as

By measuring the time period of a pendulum’s oscillations, you can estimate the value of pi from direct observation. The equation assumes that the motion is simple harmonic — which is only true when the angle is small enough that sin⁡θ ≈ θ — so it works best for small release angles. You may notice that the time period is evidently independent of the mass of the bob. In realistic systems, however, the string will have a nonzero mass and the bob will have a nonzero radius. This means that the accuracy of the estimated value of pi improves when the pendulum mimics an ideal case, i.e., when the string is massless, when all of the mass is concentrated at a point-like bob, and when the initial angle is small (<15 degrees). A pendulum comprising a string of length L and a bob of radius R_b being released at an initial angle θ.

Building the Simulation App

A simulation app based on a COMSOL model can serve as a predictive tool, allowing you to explore how different design parameters, such as pendulum length or mass, affect the results before performing the actual experiment. To build a model of a simple pendulum, we’ll use the Multibody Dynamics interface, available in the Multibody Dynamics Module add-on product, as shown in this example model of a double pendulum. In the model the Events interface is used to track the number of oscillations in the simulation using a Discrete State variable called count. The process occurs in two stages, represented by two implicit events. When the bob’s velocity in the x direction becomes negative, count is incremented by a positive fractional value less than one. When the x-component of the velocity becomes positive again, count is rounded up to the next integer using the ceiling function. A stop condition is implemented as another implicit event that is triggered when count reaches the specified number of oscillations. Once the simulation terminates, the total simulation time (T_sim) and the number of oscillations (N) are used to calculate pi according to the equation

Now let’s use the Application Builder to create an app that estimates pi for different configurations of the simple pendulum. You can easily build custom simulation-based apps without needing extensive programming by using the drag-and-drop form design and the Method Editor for short code snippets. The app serves as an interactive tool to estimate pi for different pendulum parameters by computing the time period. This approach makes it especially useful in educational settings. For instance, a teacher could use the simulation to design a pendulum of appropriate size and timing for classroom demonstrations, helping students connect theoretical predictions with hands-on measurement. 

The app’s UI.

The app contains sliders to control the length of the pendulum, the radius of the bob, the initial angle of release, and the number of oscillations to solve. It also contains a field where you can provide a numeric value of the ratio of the bob’s density to that of the string. 

Once you have chosen the desired parameters for your pendulum, clicking the Plot Pendulum button displays what your pendulum looks like based on the chosen values. The Compute button can then be used to simulate the pendulum, during which the Graphics window and the energy plot are updated in real time. After the solve, the estimated value of pi and the error from the true value are displayed to reflect the results obtained from the current solve. Feel free to share the set of input values that gave you the best estimate of pi in the comment section below!

A screen recording of the app in use.

Building an app like this in the COMSOL Multiphysics^® software using the Application Builder can be done by leveraging templates and without extensive programming experience. This way, engineers and scientists can link theory and simulation practically to make their models accessible to users who are not familiar with modeling.

Next Steps

You are welcome to download the MPH file with the app design and related files from the Application Gallery via the button below!

Further Reading

Want to see how simulation apps are being used in the real world? Check out a few examples below:

• Read about a smartphone app powered by multiphysics simulation: Forecasting Fruit Freshness with Simulation Apps
• Take a look at how a university professor uses apps within the classroom: Bringing Lab Courses to Remote Learning Students with Simulation Applications

Reference

1. Borwein, J.M. and Borwein, P.B., 2004. Ramanujan and Pi. In Pi: A Source Book (pp. 588-595). New York, NY: Springer New York.