Table of Contents

1. Formats for Importing Surface Meshes
2. Repairing and Editing Imported Surface Meshes
3. Demonstrating with a Skull Mesh:
   1. Moving and Rotating Imported Meshes
   2. Repairing Holes and Problems Like Intersecting Elements
   3. Making a Cutout of the Mesh
   4. Removing Elements from Meshes
   5. Remeshing Faces to Control the Element Size
   6. Merging Meshes, Uniting Meshes, and Creating Computational Domains
   7. Joining Domains and Creating Selections
4. Concluding Thoughts

Formats for Importing Surface Meshes

The 3D surface meshes I’ll focus on in this blog post often originate from 3D imaging. In a medical application, for example, a 3D surface mesh can be used to describe the geometry of an internal organ or even the whole body. Common file formats supported for import into COMSOL Multiphysics^® are STL, 3MF, and PLY. STL is perhaps the most widely recognized format since it is also used for 3D printing, so I’ll use it for discussion here. However, note that the mesh editing methods I’ll present in this post will work for any mesh file in a supported format imported into a Mesh-Based Geometry sequence.

STL files do not contain any boundary or domain information, so this will need to be created either during import, which is done automatically, or manually after import. Boundaries are always created during import. However, domains — i.e., volume regions you can assign materials to and fill with a volumetric mesh — cannot be created during import if the imported mesh has intersecting elements due to overlapping regions (which you’ll see an example of shortly). Most physics interfaces in the software require domains, so it is important to create them. For meshes with overlapping regions, you need to first resolve the intersections so that there are edges tracing all of them. Then, you can manually create the domains. Filling the domains with a volume mesh is done when the computational mesh is generated, which I will demonstrate in Part 2 of this series.

Note: For designs created in CAD software, it is recommended to export them in any of the formats supported for import with the CAD Import Module.

Repairing and Editing Imported Surface Meshes

In addition to the challenges of working with the STL format, another common challenge that arises when importing a surface mesh is that the quality of the triangular mesh is often too poor to be readily used for finite element simulation. From a simulation perspective, the quality of a mesh is poor if it contains:

• Holes or self-intersecting boundaries preventing domain creation
• Thin spikes or other unwanted irregularities in shape
• Triangles with very acute or obtuse angles
• Triangles that differ greatly in size

COMSOL Multiphysics^® supports a range of operations for handling all of these challenges, as well as for making other sorts of common mesh edits. You can use these capabilities to:

• Move, scale, and rotate the imported mesh
• Modify the repair tolerance during import
• Fill small or large holes in the mesh
• Remove elements from the mesh, replacing them with new meshed surfaces
• Combine and merge meshes
• Combine imported mesh with parametric CAD designs (or geometry drawn in COMSOL Multiphysics^®) to run parametric sweeps
• Intersect multiple imported meshes with each other or with a plane
• Remesh the surfaces to control element size and the shape
• Use a curved surface representation of nonplanar boundaries

Editing done to the surface mesh of a skull.

In this blog post, I’ll use a mesh of a human skull to demonstrate the workflows for several of these actions.

Demonstrating with a Skull Mesh

The mesh for the following demonstrations was obtained from the NIH website. (Ref. 1) Here, I’ve imported the mesh into a Mesh-Based Geometry sequence, which is a sequence that defines the geometric model with the domains and boundaries needed to set up and solve any physics. Alternatively, you can do the import, repair, and editing in a Mesh Part sequence, as is demonstrated in Part 2 of this blog series (coming soon) and in the STL Import Tutorial Series. Moving portions of the work to one or several mesh parts is good for organizing it but also for when you want to reuse the same mesh in several sequences or use mesh parts as construction geometry (also shown in Part 2).

For this example application, I’ll need to rotate and move the mesh to get it into position, remove some elements that are not needed, and calculate the intersection between the meshes of the bone and teeth. I will also touch on some other useful operations and workflows that can help in other cases. To aid with visualization, the mesh surfaces are colored as in the image above, unless stated otherwise. In the example model that can be downloaded at the end of this post, the boundaries remain gray until the colored selections are created at the end.

Moving and Rotating Imported Meshes

Many 3D imaging techniques are performed on the bodies of humans who are lying down, which means that an imported mesh may need to be rotated to, for example, put a body model in an upright position. Moving and rotating an imported mesh is done by adding one or several Transform attributes to the Import node. Rotation around multiple axes requires several Transform attributes.

For this example, I implemented a rotation by first aligning the z-axis to the axis specified by the Axis type setting; then, the mesh was rotated 83 degrees around the specified x-axis, and I added a displacement to place the mesh at the origin.

The orientation of the mesh in the file (left). By setting the Axis type to x-axis, the mesh is rotated to align the z-axis with the x-axis (middle). Lastly, the rotation angle of 83 degrees around the x-axis is applied (right).

Note that it is recommended to add any Transform attributes before doing any other editing or repair, as reimporting the mesh with new Transform settings will likely result in lost input selections in downstream features in the sequence.

Repairing Holes and Problems Like Intersecting Elements

During import or when building other operations, Information nodes may appear to notify you that there are holes in the mesh and provide you with a selection of edges that you can zoom in on.

To repair smaller holes relative to the overall size of the mesh (typically the size of one or a couple of mesh elements), it is often sufficient to increase the Repair tolerance of the import, as this can collapse elements or align a mismatch in vertex coordinates. For larger holes, or if increasing the repair tolerance doesn’t seal off all holes, try using the Fill Holes operation. If doing so fills all of the holes and the surfaces thus form watertight regions, the Fill Holes operation will also create domains inside the regions.

For problems like intersecting elements, as seen below at right, isolate the problematic region, delete the elements, and fill in the hole. For more details about this workflow, refer to the Removing Elements from Meshes section of this blog post.

The skull mesh at hand doesn’t require any repair of this kind, but the STL Import Tutorial Series includes several examples of smaller holes and intersecting elements, as can be seen in the following images. Download the application files to learn how to fix these.

A hole with zero area (i.e., a slit) is bounded by edges highlighted in blue (left). A hole in the mesh caused by a mismatch in vertex coordinates (middle). Intersecting elements in an imported mesh (right).

Making a Cutout of the Mesh

What if you only want to include a smaller region of the imported mesh in the simulation? For example, with the skull mesh, I ultimately wanted to create a model to simulate a dental implant, so I didn’t need the full skull mesh. Assuming that parts of the mesh are not needed for the simulation, you can delete the parts of the mesh that you are not interested in. The skull mesh was imported as several mesh surfaces, which made it possible to extract the surfaces of the jaw bone and teeth using a Delete Entities operation, as shown in the following image.

I used a Delete Entities operation to delete the parts of the mesh that I was not interested in; I only kept the upper jaw and two teeth (highlighted in light yellow).

Next, I used Intersect with Plane operations to cut the upper jaw bone into smaller pieces and then delete the parts that weren’t needed for the simulation. You can specify the planes using coordinate planes, points, coordinates, the normal vector, and more. Once a plane has been specified, one or several parallel planes can be added, as you can see in the bottom two of the following images. The Intersect with Plane operation is used to create selections of what is above and below the plane. You can use these selections to easily select what to delete.

Clockwise from top left: I used multiple Intersect with Plane operations to cut the upper jaw into smaller pieces, deleting what was above and/or below the planes as I went.

Removing Elements from Meshes

For the skull mesh, I wanted to create domains inside the bone and teeth to assign materials to when adding physics to the model, and to do this I needed to resolve the intersection between the bone and teeth surfaces with a Union operation. As you can see in the following image, the surfaces intersect in a complex way, with small overlaps and gaps between the meshes. I also needed to resolve the intersection to remove the indents from the missing tooth because I wanted to simulate that the dental implant (to be added to the mesh in Part 2) had fully healed into the bone. It is easier to unite meshes where the intersection is reasonably close to 90 degrees rather than the small angles that the skull mesh had. Therefore, I removed a rather large portion of the bone surface and replaced it with a flatter new meshed surface.

Using a Clip Plane operation to see inside the mesh. The surface of the bone (light yellow) contains depressions that intersect the teeth (white) and the future dental implant.

The workflow for removing elements typically consists of isolating the irregularity, deleting it, creating a new surface, and joining it with an adjacent boundary. Some STL meshes contain a large number of problematic regions, so isolating the surface can been automated by using the Mesh Partition with Ball add-in, which is available from the Developer tab in the software. The partitioning operations can quickly partition a mesh but are restricted to specific shapes or logical expressions by their respective settings. On the other hand, the Create Edges operation is more flexible, as you select each mesh edge manually. For practical reasons, the Create Edges operation is typically only used when isolating a smaller number of mesh elements.

I’ll demonstrate a workflow here that uses both types of operations and starts with loading the Mesh Partition with Ball add-in. To load the add-in into a model, click the Add-in Libraries button on the Developer tab. In the list of add-ins, select the checkbox for the Mesh Partition with Ball add-in, and click Done. Now, click the Add-ins button on the Developer tab and select Mesh Partition with Ball.

When the add-in is loaded into the model, middle-click with the mouse to place the center of rotation, which will be used as the center for the ball. Click the Create button in the Settings window for the add-in to add a Partition with Ball operation to the sequence with its center and radius filled in by the add-in. The settings can be changed if needed. In the example, I opted to only partition the largest bone boundary (colored blue at left in the following set of images), as that was the only boundary that contained elements I wanted to replace. For other use cases of the add-in, see the Spray Particle Deposition in Human Airways tutorial model.

Looking at the mesh from below with the teeth hidden for visibility. In the Graphics window, middle-click the mesh to set the center of the ball. The radius of the ball (light pink) is set automatically by the add-in. The teeth are hidden for visibility (left). The resulting mesh after a Partition with Ball operation has been added to the sequence to isolate a region enclosed by the ball (right).

For the rest of the depressions I used the Create Edges operation, as the bone wall around the depressions in the middle was quite thin and clicking the edges manually gives you more control over which elements to isolate. With the Create Edges operation, you click the mesh edges in the Graphics window to create geometrical edges that will isolate a surface. These edges are selected to attach to the edges created by the add-in (the thicker black edges in the following image at left). Use of the Create Edges operation is also demonstrated in the STL Import Tutorial Series.

Click the mesh edges (highlighted in blue) in the Graphics window to form the bounding edges of the surface containing the rest of the depressions (left). Looking at the mesh from the side using a Clip Plane operation: The highlighted surfaces (blue) are selected to be deleted and replaced (right).

Once the depressions were isolated, I deleted the surfaces and created a new meshed surface using the Create Faces operation, as shown in the following images. The generated surface is a minimal surface, which means that it aims to be as planar as possible with a coarse mesh size. If the hole had been small compared to the bounding box of the mesh, I could have used the Fill Holes operation instead.

The highlighted edges (blue) are used as input to the Create Faces operation (left). In gray, the new meshed boundary (right).

Finally, I used the Join Entities operation to join the larger boundary of the bone with the newly created boundary.

Remeshing Faces to Control the Element Size

To help the Union operation calculate the intersections between the meshes of the teeth and bone, it’s ideal to have mesh elements of more equal size, as well as high-quality elements on the intersecting faces. Achieving these conditions can be accomplished by using the Remesh Faces operation. To gain more control over the element sizing, change the element size in the Size subnode. To get a similar size on all surfaces, regardless of the curvature and any narrow regions, set the Minimum element size and Maximum element size to similar or equal values. Note that the mesh generated when using the Remesh Faces operation is only to help the Union operation and to get a smoother representation of the surfaces. The Remesh Faces operation accomplishes the latter by placing new nodes on a curved representation of the surface derived from the linear mesh. The computational mesh is generated later in a Mesh sequence.

Remeshing the teeth and the joined surface of the bone to prepare them for the Union operation.

The Remesh Faces and Remesh Edges operations can also be used to control the shape of surfaces and edges. For an example where the Remesh Faces operation is used to smooth the mesh surfaces, see Generating a Simulation Mesh of a Femur From 3D Data.

Merging Meshes, Uniting Meshes, and Creating Computational Domains

Going back to the example of the dental implant, I used the Union operation to calculate the intersection between the overlapping meshes of the teeth and bone. The Union operation calculates the intersection edges, partitions the surfaces, and splits the mesh elements as necessary, as you can see in the following image at right.

If the Union operation reports problems with intersecting elements, keeping domains, or the like, try another mesh size in the Remesh Faces operation, try lowering the Absolute repair tolerance of the Union operation, and/or switch to the Linear option in the Placement of mesh vertices setting. When the Curved placement of vertices option is used (as shown in the middle image below), new mesh vertices are placed on a curved representation of the mesh surface that is derived from the linear mesh. The Linear option simplifies the problem, with the drawback of having more angular intersection edges.

For the example model, I used the Linear placement of mesh vertices option for the Union operation, which created edges tracing the intersections which the mesh conformed over (image at right).

The mesh before building a Union operation: The triangles of the tooth are intersecting the triangles of the bone (left). When the Curved placement of mesh vertices option is used, the intersection edge is smoother (middle). When the Linear placement of mesh vertices option is used — as I opted to do in the example model — the intersection edge is more angular (right).

Until this point, domains had not been created inside the mesh because of the intersecting elements between the teeth and bone. Now that the intersections in the mesh were resolved, I used the Create Domains operation to form the computational domains inside the watertight regions. To get a visual of the domains, you can right-click anywhere on the mesh in the Graphics window, select Clipping and then Clip Plane. Then, in the Graphics toolbar, click the Clipping Active button and select Show Cross Section.

An alternative way to see whether the domains have been created is to check if the Selection List reports any domains. After the Create Domains operation has been built, the Graphics window should show cross-sectional faces for the domains (following image at right), the Selection List should be populated on the domain level, and the Information section of the Create Domains operation’s settings will list the number of created domains as well as the total number of domains.

Before building the Create Domains operation, there are no domains in the mesh — only boundaries, edges, and vertices (left). After the Create Domains operation has been built, the domains are visualized with cross-sectional faces in the Clip Plane if the Show Cross Section setting has been enabled (right).

Meshes imported using several Import nodes are treated the same way as different objects are treated in the Geometry sequence when using the Form Assembly method. To connect touching boundaries in such meshes, it is recommended to use the Merge Entities operation. This operation can also be used to collapse gaps in a mesh, as shown in the following example, which is taken from the STL Import Tutorial Series. Download the application files to learn more.

Merging two faces (yellow) to seal a gap between a vertebra and disc (left). The resulting mesh after the gap has been sealed (right).

Joining Domains and Creating Selections

As expected, the domains of the teeth in the example mesh were partitioned with the surface of the bone (following image at left). I could have kept the domain partitioning as it was, but instead I used the Join Entities operation to join the domains of the teeth into two domains (image at right). After that, I could define domain and boundary selections to be used when setting up materials and physics. For this example, I created domain selections to facilitate assigning material properties.

Joining the domains of the teeth (left). The domains with assigned colored selections: Now that the domains have assigned colors, the cross section also shows these colors (right).

The final mesh-based geometry.

This mesh-based geometry can very well be used for simulation at this point; all that’s needed is to build the computational mesh in the Mesh 1 node.

Concluding Thoughts

COMSOL Multiphysics^® offers multiple features to repair and edit imported meshes. In this blog post, I have demonstrated how to rotate an imported mesh, remove elements, merge entities, unite meshes, and create computational domains. If you haven’t already, check out the model file of the example featured in this post. Note that all edits done for this mesh may not be necessary for other meshes. The order in which you add the operations can also differ, depending on the mesh at hand. I have prepared this mesh to later add a dental implant where a tooth is missing, and the idea is to generate a computational mesh that could be used for a solid mechanics simulation. The topic of combining the mesh prepared in this blog post with imported CAD geometry is covered in Part 2 of this series (coming soon).

Learn More

If you are interested in learning more about working with imported meshes, check out these resources:

• Watch this video on editing and repairing imported STL meshes and combining them with CAD geometry.
• STL Import Tutorial Series, a tutorial model series showcasing how to repair, edit, and combine STL meshes with each other and with geometry created in the software
• Analyzing Porous Structures on the Microscopic Scale, a tutorial model that showcases import of an STL file and the generation of a physics-controlled mesh
• Spray Particle Deposition in Human Airways, a tutorial model featuring the Mesh Partition with Ball add-in to repair the STL mesh
• Generating a Simulation Mesh of a Femur From 3D Data, a blog post showing two ways of creating a mesh from data in a file and smoothing the shape of the meshed surfaces

Reference

1. kbrowne, “Skull and Eyes – Visible Human Male. Version 1.03,” NIH 3D, 12, Apr. 2026; https://3d.nih.gov/entries/3DPX-020591.

This is the second blog post in a two-part series on defining load cycles in battery models. Read Part 1 here.

Introduction

Battery cycling scenarios can involve complex, multistep protocols and are not always as straightforward as applying a constant current for charge and discharge with switching based on a voltage cutoff. In practice, batteries may operate under current, voltage, or power control, or a mix of all three, specified either by constant or variable inputs or by tabulated data. Transitions between modes can be governed by simple duration conditions or more custom, user-defined criteria based on a variety of performance outputs. Additionally, load profiles may consist of simple repeated sequences or include nested loops within a broader protocol.

Specifying load profiles is critical for accurately modeling real-world battery system performance.

When incorporating these profiles into a battery model, it is equally important to consider their numerical implications beyond defining and customizing them. Abrupt changes and complex switching conditions in the load profile can lead to numerical stability issues during simulation. As discussed in Part 1, having an event-based profile, either using the Events interface directly or through the predefined Charge–Discharge Cycling feature, is often the best approach for solver behavior during sudden load transitions. However, the Charge–Discharge Cycling feature is limited to constant current–constant voltage (CCCV) profiles, with or without rest periods, where switching is based on current and voltage thresholds. While the Events interface does not impose restrictions on the types of cycling scenarios that can be defined, implementing load profiles with it requires a certain level of expertise and can become cumbersome for more advanced cycling cases. To address this, the Load Cycle feature, an event-based functionality, was introduced in version 6.4 to simplify the process. It is designed to be both flexible and comprehensive, allowing users to define a wide range of cycling scenarios while maintaining numerical robustness.

Where to Find the Load Cycle Node in the Model Tree

Where to find the Load Cycle model tree node depends on the physics interface used in the model. In detailed battery models, such as the Lithium-Ion Battery and Battery with Binary Electrolyte interfaces shown in the screenshot below, as well as the Current Distribution interfaces, the Load Cycle node is available as a boundary condition. In such models, since the electrodes are explicitly represented, the negative side is grounded, and the load profile is assigned to the positive side to reflect the operating conditions. In cases where Electrode Surface, Thin Porous Electrode, Perforated Electrode Surface, or Highly Conductive Porous Electrode nodes are present within these interfaces, Load Cycle has been included as an option in the Electrode Phase Potential drop-down menu of these nodes. For simplified models, Load Cycle is included as an Operation Mode in the Lumped Battery interface and as a boundary condition for current conductors in the Battery Pack interface.

The Load Cycle feature can be found as a boundary condition by right-clicking on the Lithium-Ion Battery and Battery with Binary Electrolyte interfaces or selecting from the ribbon menu.

What the Load Cycle Feature Offers

In the Load Cycle Settings window, the load type can be selected. In the Load Type section, you can choose between Galvanostatic, Potentiostatic, or a combination of both, as shown in the screenshot below. Selecting either galvanostatic or potentiostatic limits the available operation modes to Current and C-rate, or to Voltage, respectively, along with the Rest and Subloop child nodes. Choosing the combined option Potentiostatic and galvanostatic enables access to all modes, including power.

Transitions between added operation modes can be based on multiple criteria, including elapsed time, cutoff values (such as voltage or current), or any user-defined expression. If the Use elapsed time only option is selected at the top level, transitions are restricted to elapsed-time-based (explicit) switching. Otherwise, additional options are also available to control transitions between modes.

The Load Cycle Settings window includes options for terminating the load cycle, such as limits based on total cycling time, number of cycles, voltage thresholds, or user-defined criteria. This built-in functionality allows users to define the end of cycling without manually introducing variables or adding stop conditions in the solver configuration. These conditions can be selected from the Cycling Stop Condition drop-down menu in the Load Cycle Settings window.

The settings also allow users to enable built-in global probes, which automatically monitor voltage and current. Using probes in the model, regardless of the physics involved, offers advantages, such as allowing you to view results without waiting for the simulation to finish, to use probe variables (which are global) in the results section, or to define expressions within the model. In battery modeling, in particular, monitoring voltage and current during the simulation is highly recommended, as it helps understand how the battery behaves under the load cycle and assists with troubleshooting.

The Load Cycle Settings window shows the load type selection, the option to restrict the continuation method to elapsed time only, and settings for cycling stop conditions and voltage and current probes.

Using the Load Cycle Feature

Any load profile can be defined by specifying its steps, the characteristics of each step, the transitions between them, and the conditions for stopping the cycling sequence. With the Load Cycle node added to the Lithium-Ion Battery interface (or to any Electrochemistry interface) and the load type set in the Load Cycle Settings window, construction of the load profile begins by adding the profile steps in sequence: right-click the node and add the corresponding step. Once the steps have been added and arranged as desired, each step can be customized by applying a constant value or assigning a function to capture that step’s characteristics. The input type for each step can be selected in the corresponding Settings window from a drop-down menu, with options for value, function, or step sequence to define the input values. Each step of the load cycle is followed by the next when the switching condition, defined through the Continuation Condition, is met and enforced using solver events.

Although using functions to define the load profile was discussed in Part 1, incorporating them into this event-based framework improves numerical stability when the load is applied to the model since function smoothing is not required. Alternatively, you can use a step sequence, which enables importing or locally define a table of times at which the set value is updated. This option is particularly useful when importing a user-defined load profile, such as a drive cycle, from a text file into COMSOL^®. The screenshots below represent two different load profile examples, each constructed with different steps and input types.

A 1C constant-current charge and discharge profile with known durations, with a rest period in between ad final rest step is defined using the Load Cycle feature.

A drive-cycle based on a table of time versus C-rate is defined by selecting “Step sequence” as the input type for the C-rate.

To cover scenarios that require a nested loop within a main sequence, the Subloop child node can be added. This feature is only accessible if Use elapsed time only has not been selected in the Load Cycle settings. The Subloop contains its own nodes for different modes of operation, similar to the main sequence, representing a loop within the overall protocol. The Subloop iterates over these modes, and its duration is controlled by a selected break condition. This condition can be based on elapsed time, number of cycles, or a user-defined criterion. Once the condition is met, the simulation returns to the main sequence or, if the subloop is the final step, proceeds to the first node under the Load Cycle node. The screenshot below shows a simple profile with a subloop. Such a profile can be easily extended to include multiple subloops with completely different modes of operation and termination criteria.

The Subloop feature defines a profile consisting of two cycles of charge and discharge pulses, followed by a charge pulse and a rest period.

Conclusion

Although all the methods for defining load cycles covered in the Part 1 remain valid, the Load Cycle feature introduced in version 6.4 offers a more straightforward and numerically robust approach that also supports complex, custom profiles. Key aspects of this functionality have been covered in this blog, and users can now start using it in their battery simulations. Most examples in the Application Library have been updated to use this feature and can serve as demonstrations of different scenarios implemented with it.

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”.