Editor’s note: Read the previous blog post, “Minimizing the Charging Time of a Battery”, to learn more about the relationship between charging time and battery longevity for a lithium-ion battery.

Modeling the Energy Density of a Lithium-Ion Battery

Let’s consider a lithium-ion battery model constructed using the Lithium-Ion Battery interface. The cell is simulated under a constant-current 1C discharge, and the model can be used to analyze the battery’s electrochemical behavior and performance. The cell consists of a graphite negative electrode and an NMC111 positive electrode, separated by an electrolyte-filled separator. The electrode thicknesses and active material volume fractions are specified as part of the porous electrode properties. For a given set of these parameters, the model predicts the cell voltage throughout the discharge and, consequently, the amount of energy the cell can deliver.

The energy delivered during discharge can be calculated by integrating the product of the cell voltage and applied current density over the discharge time:

The time integration can be implemented as an ordinary differential equation (ODE) in COMSOL Multiphysics version 6.4 and earlier, while in COMSOL Multiphysics version 2027, the functionality will be built into the Time Integral feature.

In this example, the cell current density, , is set to correspond to a constant C rate, but, in principle, the maximization could be performed for any given discharge load profile.

Optimizing Electrode Dimensions

The electrode dimensions can be varied by changing the geometry and remeshing in every optimization iteration. However, by using the Deformed Geometry interface instead, the geometry can be modified by deforming the computational mesh, which allows the electrode dimensions to vary continuously without remeshing for each design configuration, helping maintain a smooth dependence of the model response on the geometric design variables. This condition makes the problem well suited for gradient-based optimization, which is particularly advantageous when working with many design variables.

In practice, the charges in the two electrodes should be balanced, and therefore the negative electrode thickness, , is expressed in terms of the thickness of the positive electrode, :

where:

• = lithium inventory loss factor during formation
• , = negative and positive (initial) electrode volume fractions
• , = maximum negative and positive concentrations
• = maximum lithiation level in the negative electrode
• = minimum lithiation level in the positive electrode

To gain a deeper understanding of how to balance electrodes and see other example models, view this blog post.

The battery design variables are limited to the thickness of the positive electrode, , and the initial volume fractions of the negative and positive electrodes, and .

Solver and Optimization Settings

The partial differential equations (PDEs) for the electrolyte potential, electrolyte salt concentration, and the electric potential are solved with a time-dependent solver. Here, the backwards differentiation formula (BDF) is used, then terminated when the insertion particle concentration at the boundary between the negative electrode and the separator falls below the minimum lithiation level. The method of moving asymptotes (MMA) is used as the optimization solver. Robustness is improved using the Move limits option, which limits the change of design variables between optimization iterations.

Results

The objective history is plotted with and without shape optimization in the figure below. The shape optimization improves the objective from 2.5 to 3.2 J/m³, although this improvement depends on the initial values of the design variables.

The objective history shows that 50 iterations of the MMA solver are sufficient.

The thickness of the positive electrode and the positive and negative electrode volume fractions are plotted below. Compared to the initial values, the optimized battery has significantly denser electrodes (that is, with larger electrode volume fractions). In this way, the energy density gets larger since the electrodes can store more lithium per volume. The electrode volume fractions may however not become too high, since this at the same time reduces the electrolyte volume fraction of the porous electrodes. This in turn reduces the effective ion transport capabilities of the electrolyte.

Similarly, slightly thicker electrodes can contribute to a higher energy density in the optimized battery by increasing the amount of active material. However, the electrodes should not be made too thick, since increasing the electrode thickness also increases the ion transport distance through the electrolyte phase. The resulting transport limitations can lead to larger concentration gradients and increased polarization, reducing the usable capacity and energy of the cell.

The optimization results in higher active material volume fractions and slightly thicker electrodes.

Takeaways

In this example, optimization is used to improve the battery energy density by varying the active material volume fractions and the thickness of the electrodes.

More generally, COMSOL Multiphysics facilitates gradient-based optimization of fully coupled multiphysics problems. Different optimization approaches can also be combined within the same model, including parameter, shape, and topology optimization as well as optimal control.

Next Steps

• Interested in the details of this model’s setup in COMSOL Multiphysics?
  • Download and explore it for yourself: Maximizing the Energy in a Lithium-Ion Battery
• Want to learn more about optimization?
  • View this tutorial model from the COMSOL Application Gallery about minimizing a lithium-ion battery’s charging time.
• Check out these previous blog posts.

1. Grain Packing

Industries such as agriculture, food & beverage, and pharmaceuticals involve grain packing. Packing fraction — the volume occupied by the grains divided by the total sample volume — depends on grain shape, size, and interaction forces. Understanding how these factors affect packing efficiency after the container is sealed is important for the safety of the product and the optimization of packaging processes.

The first tutorial model we’ll highlight demonstrates how grains settle under gravity to form a packed bed. The model geometry consists of a rectangular box with a length and width of 0.2 mm and a height of 0.3 mm. In the box, 5000 grains are released, with the diameters of each grain being drawn from a normal distribution with a mean of 10 μm and a standard deviation of 2.5 μm. The grains are allowed to settle under the influence of gravity and the study’s specified contact forces for 8 milliseconds.

Three contact force models are used to study the effects of different forces on the packing efficiency of the grain:

1. Hertz–MD, one of the most common force models used in discrete element method (DEM) simulations
2. Hertz–MD with adhesion, which extends the Hertz–MD model with surface adhesion, allowing tensile forces to prolong contact after a collision
3. Hertz–MD with adhesion and van der Waals forces, which extends the previous force model by adding the (noncontact) van der Waals forces between grains, which is important when grain sizes approach the micron scale

Comparing these force models illustrates how adhesive interactions affect the predicted packing fraction.

To reduce the computational cost of these simulations, the Granular Flow interface uses a non binary search (NBS) algorithm to identify potential grain–grain contacts and an octree search algorithm to identify potential grain–wall contacts in 3D. These algorithms limit the grains and wall elements checked for each grain during the time-dependent study. The search radius defaults to the grain’s physical radius and can be increased using the contact search expansion ratio in the Grain Properties feature.

Left: the packed grains, using the Hertz–MD model. Middle: the packed grains, using the Hertz–MD with adhesion model. right: the packed grains, using the Hertz–MD with adhesion model with van der Waals forces enabled.

Including adhesion lowers the final packing fraction in this example, while adding van der Waals forces produces a smaller, additional decrease.

In this example, the Hertz–MD model predicts a significantly higher packing fraction than the models that include adhesive interactions. The Granular Flow Module can account for adhesive and van der Waals forces when these interactions significantly influence grain packing.

2. Hopper Flow

Hoppers that store and discharge granular materials are used in the agriculture, food processing, manufacturing, mining, and construction industries. The controlled material discharge from a typical funnel-shaped hopper often exhibits some combination of two flow patterns: mass flow and funnel flow.

In mass flow, the entire granular bed moves when material is discharged. This pattern typically follows first-in, first-out behavior, with the lower layers discharging before the upper layers and limited mixing between them.

In funnel flow, grains move through a flow channel while surrounding regions remain comparatively stagnant. Material from different layers can mix during discharge, and some material may remain near the walls. Hoppers that exhibit this type of flow have the advantage of having a smaller footprint and thus may be preferred in applications where space is a constraint.

Left: mass flow. Right: funnel flow.

This example model considers how the shape of a hopper design influences the flow pattern and thus the results. Two hopper designs are modeled, with hopper walls inclined at 15^° and 75^° from the vertical. For each design, a total of 150 grains with equal diameters of 10 mm are released at the inlet boundary every 20 ms up to 0.58 s, leading to a total of 4500 grains. The grains are allowed to settle under gravity for the first 1.0 s, after which the outlet is opened to allow the discharge of the stored material for a total of 4.0 s.

The time-dependent simulations show mass flow for the hopper with walls inclined at 15^° from the vertical and funnel flow for the hopper with walls inclined at 75^° from the vertical.

A Parametric Sweep study compares the two hopper designs, while coloring grains by release time helps visualize the discharge sequence and mixing between layers. The Bounding Box feature removes grains that have left the region of interest to improve computational efficiency.

The mass remaining in each grain layer for the hopper with walls inclined at 15^° from the vertical.

The corresponding results for the hopper with walls inclined at 75^° from the vertical.

3. Angle of Repose

The angle of repose is the angle between the surface of a stable granular heap and the horizontal plane. This angle often depends on a variety of factors, including the size and shape of the particles, surface roughness, friction coefficients, and density. The angle of repose is therefore a complex property of the granular system, and it is often used in material characterization.

This model uses a similar hopper geometry as in the Hopper Flow model. Here, a total of 5000 grains, 20 mm in diameter, are deposited through the hopper outlet and onto the surface below. Once the grain heap stabilizes, the grain positions are used to compute the angle of repose. The grains are colored by speed in m/s, with their very low speeds indicating that the heap has stabilized.

The heap of grains formed after emptying the hopper.

The outline of the left and right sections of the free surface, along with their least-squares fit, are shown in the figures below. The least-squares fit obtained for each section shows that the heap’s outline is fairly linear.

Left: the outline and linear fit for the left section of the heap. Right: the outline and linear fit for the right section of the heap.

4. Mixing Grains in a Rotating Drum

Many industrial applications necessitate uniform mixing of different granular materials. To achieve this result, organizations commonly employ rotating drums that include baffles to disrupt bulk flow. In this model example, we will look at how the number of baffles in a rotating drum influences the flow of materials.

The model uses the Granular Flow interface to simulate the mixing of two different grain types within a rotating drum. The extent of mixing is quantified using four mixing indices — Kramer, Lacey, Beaudry, and Valentin — calculated from the grain positions using a model method.

The mixing process with four baffles.

Mixing process with eight baffles.

The evolution of the four mixing indices provides a quantitative comparison of the two drum designs.

Evolution of the mixing indices with four baffles.

Evolution of the mixing indices with eight baffles.

For the conditions studied, the drum with eight baffles produces faster mixing and higher final mixing-index values than the drum with four baffles. These results can help teams assess drum designs before manufacturing.

5. Thermal Conduction in Granular Materials

Granular materials are a common insulator in industrial systems. The arrangement of grains influences heat transfer paths, making the insular properties of the granular mass dependent on their packing density and orientation.

Here, we will look at an example study involving a two-stage DEM simulation of thermal conduction in compressed granular materials. In the first stage, a top wall moves downward under an applied pressure of 20 kPa, compressing an assembly of spherical grains. The resulting grain arrangement provides the starting point for the second stage, which examines heat transfer through the bed. This process reflects manufacturing applications in which compaction is used to increase packing density and modify a material’s mechanical and thermal properties.

Grains undergoing compression.

In the second stage, the compacted granular bed is subjected to thermal conduction with a fixed temperature boundary condition of 300 K set for the bottom wall. Grain motion is frozen during this stage, assuming that the compacted bed remains mechanically stable. Only conductive heat transfer through grain–grain and grain–wall contacts is considered.

Solving the thermal stage separately allows larger time steps than those needed to resolve grain contact dynamics, reducing computational cost. In this example, 8 hours of thermal evolution are simulated in approximately 30 seconds of computation time.

Thermal conduction in the granular material over eight hours.

6. Powder Spreading Using a Moving Roller

Additive manufacturing techniques such as powder bed fusion (PBF) rely on obtaining thin uniform layers of powdered materials. The quality of the powder bed depends on its packing density, uniformity, and thickness. These properties are influenced by the roller geometry and motion. This process involves heaping the powder in front of the bed and employing the moving roller to flatten the powder evenly over the area. Different roller speeds and geometries can be used to achieve different qualities in the bed.

In this model example, 20,000 grains (all 1 mm in size) are released onto the bed and allowed to settle. Then, the roller is set in motion with a translational speed of 100 mm/s for a total linear traversal length of 163 mm. The roller also rotates with an angular speed of 6 rad/s around its axis, resulting in the formation of a compact powder bed. The roller’s combined translational and rotational motion is prescribed using the Wall Movement settings. Sequential studies simulate grain settling, heap formation, and spreading.

The moving roller flattening the grain heap into a uniform powder bed.

The resulting bed is analyzed using packing fraction, grain-height distribution, and coordination number — the number of contacts each grain has with other grains. The packing fraction remains approximately constant at 0.47 along the bed, indicating relatively uniform packing.

7. Screw Conveyor

Screw conveyors can transport bulk granular materials horizontally or over an incline while also incorporating other specialized tasks, such as heating, cooling, and drying. These conveyors are widely used in food processing, pharmaceuticals, chemical industries, agriculture, cosmetics, and plastics, among other areas.

In this example, we’ll look at a screw conveyor model that highlights how the Granular Flow Module can be used to study the behavior of grains depending on how they are released into a screw conveyor. Moreover, this model highlights how the module can be used to build complex geometries for granular flow applications.

The screw conveyor transporting grain.

In this example, the grains have a density of 4200 kg/m3 and a diameter of 10 mm. They are released in 10 batches of 2500 grains, in intervals of 0.5 seconds. The screw is rotated at a constant angular speed of 120 rpm, and the simulation lasts 20 seconds. The screw rotation is prescribed using the Wall Movement settings, and a Force Accumulator subnode evaluates the forces exerted by grains on the screw. A model method calculates the mass flow rate through a plane within the conveyor.

Grain positions at the end of 10 seconds.

The image shows grain positions at 10 seconds, colored by release time. Most grains from each batch remain clustered, although partial mixing occurs during transport. Some grains remain near the inlet because of the conveyor’s inclination and the clearance between the screw blade and casing. The throughput and screw-load results can support comparisons of rotational speed, conveyor inclination, and screw geometry.

Next Steps

Learn more about the Granular Flow Module and how it can be used to simulate the behavior of granular materials:

Explore the Tutorial Models

Ready to explore the models featured in this blog post? Find them in the Application Gallery, along with many other models for simulating granular flow:

• Grain Packing
• Hopper Flow
• Angle of Repose
• Mixing Grains in a Rotating Drum
• Thermal Conduction in Granular Materials
• Power Spreading Using a Moving Roller
• Screw Conveyor

In this blog post, we’ll look at several features and options in COMSOL Multiphysics that support accessible workflows and alternative ways that users can interact with the software.

Adjusting the Appearance of the COMSOL Desktop

The appearance of the COMSOL Desktop can be adjusted to accommodate different viewing preferences. COMSOL Multiphysics includes default, light, and dark Color Themes, which change the appearance of the user interface.

The Default, Light, and Dark themes available for the COMSOL Desktop, shown in the Preferences window.

Default mode (left) and dark mode (right).

The Graphics window also offers several additional color themes that can be selected independently of the desktop theme to help distinguish geometry, selections, and other graphical elements.

COMSOL Multiphysics also supports display scaling on all supported operating systems, enabling users to enlarge the elements of the user interface. On Windows^®, users can also take advantage of accessibility tools such as Magnifier to enlarge portions of the screen.

The COMSOL Desktop viewed using Magnifier at 200%.

Also on Windows, COMSOL Multiphysics works with screen readers such as Windows Narrator, which can announce text and information about many controls as the user navigates the interface.

The blue frame indicates the user interface element currently selected by Windows Narrator.

Keyboard Navigation

The COMSOL Desktop provides several ways to navigate the user interface and access commands using the keyboard instead of a mouse. For example, pressing the Alt key displays a KeyTip by each ribbon command that indicates the corresponding key to press to execute the command. The Tab and arrow keys can also be used to navigate between controls and through lists or trees, while Ctrl+Tab enables navigation between windows.

KeyTips shown for the Home ribbon tab, displayed by pressing the Alt key followed by the H key.

Additional options for personalizing the COMSOL Desktop are covered in the Learning Center article Personalize the COMSOL Desktop, which is part of the How to Navigate the COMSOL Multiphysics User Interface course.

The Learning Center also provides information on using keyboard shortcuts, and a complete list of available keyboard shortcuts can be found in the COMSOL documentation.

Syntax Highlighting

Syntax highlighting makes it easier to read and edit complex expressions and functions when, for example, creating parameters and variables. In expression text areas in the Settings window, variables, math functions, units, numbers, and strings are visually distinguished using different font colors and formatting (i.e., italics, bold, etc.).

Matching parentheses are also highlighted when the cursor is placed next to one parenthesis in a pair. This can be especially useful when working with expressions that contain multiple sets of parentheses.

Together, these features can make complex expressions easier to read and navigate.

The text field below the Variables table shows an expression with syntax highlighting.

Accessible Visualizations

Choosing an appropriate color table can make simulation results easier to interpret. COMSOL Multiphysics includes perceptually uniform color tables, which are designed so that equal changes in data values correspond to approximately equal perceived differences in color.

As an alternative to the Rainbow color table, which uses a broad range of colors to represent variations in simulation results, the perceptually uniform Cividis color table uses gradual changes in color and lightness and is designed to remain distinguishable for many users with color vision deficiency (CVD).

The same simulation results visualized using the Rainbow (left) and Cividis (right) color tables.

Users can also create custom color tables to control how colors are distributed across a visualization. This provides additional flexibility for creating visualizations based on the data being presented and the intended audience.

For more information about accessible color tables, see the blog post Accessible Color Tables to Aid in Accurate Perception of Results.

Alternative Ways to Access Simulation Results

In addition to graphical visualizations, COMSOL Multiphysics provides several ways to explore and analyze simulation results numerically. Results can be evaluated and displayed in tables, and numerical data, including the data used to generate plots, can be exported for further analysis. The Report Generator provides another way to present simulation results, combining text, tables, and images in structured reports.

The comprehensive COMSOL API for use with Java also enables users to interact with models programmatically, including when accessing and processing simulation results.

Temperature distribution in an electronic enclosure cooling model, with average temperatures of individual components evaluated and displayed in a table.

Accessibility Resources

For additional information about accessibility options in COMSOL products, see:

• COMSOL Accessibility
• COMSOL Accessibility Conformance Report

 
Microsoft and Windows are trademarks of the Microsoft group of companies.

Understanding Thin Conductors

The term thin conductor requires a precise meaning in the current context. We begin with the thickness of the material: It must be small compared to the other dimensions. A PCB trace or conductive tape is a typical example. The material must also be sufficiently thin that the tangential current density can be approximated as uniform through the thickness. This approximation is appropriate when the thickness is much smaller than both the radius of curvature and the skin depth. The latter can be estimated from the material resistivity and characteristic operating frequency. We further assume that conduction currents dominate over displacement currents and therefore neglect displacement-current contributions. Finally, we assume that the magnetic response of the conductor can be neglected, so no magnetic material model is needed.

A thin conductor can be modeled as a sheet current.

With the term thin conductor so defined, we can move on to addressing how such a conductor can be modeled. Starting with Maxwell’s equations, we can derive the magnetic field formulation, which uses the finite element method to solve for Faraday’s law within domains:

where the electric field is computed from:

and where the material resistivity, , is the inverse of conductivity and can depend nonlinearly on the magnetic field, the current, or both. Solving these equations also requires a set of boundary conditions and excitations, several of which are addressed in our blog post Modeling of Excitations with the Magnetic Field Formulation.

In addition, in the absence of a surface current, the tangential component of the magnetic field is continuous across an interior boundary:

where the subscripts refer to the up and down sides of the interior boundary. This condition is implicitly satisfied within the Magnetic Field Formulation interface because the tangential magnetic field is represented using curl-conforming vector elements, which maintain tangential continuity across interior boundaries. That is, by default, the interface has a single unknown at the boundary and thus enforces:

However, if there is a nonzero current flowing along a surface, then we have:

This condition is not available within the existing Magnetic Field Formulation interface, but as we are about to see, it’s very easy to add. Along with the current flowing along the surface, we also need the thickness-averaged current density in the conductor and the electric field tangential to the surface:

where is the geometric thickness. So, our objective is twofold: Introduce into the existing Magnetic Field Formulation the possibility of a difference in the magnetic field across an interior boundary due to a sheet current and introduce a constitutive relationship of the conductor carrying that current.

Here, we present an intuitive explanation that assumes some familiarity with the finite element method. Since we already have expressions for the tangential current and electric field, we can write the electrical power dissipated over the surface:

In the weak formulation, the corresponding virtual-power contribution is obtained by replacing the surface current with its test variation:

Since the surface current is determined directly from the jump in the magnetic field, its variation is likewise determined by the magnetic-field test functions. This contribution is therefore added to the variational form solved by the finite element method. Several references are provided at the end of this article for readers interested in the details. It’s also worth pointing out that AI tools can often provide specific step-by-step derivations of these types of equations.

Implementation Within COMSOL^®

The first requirement for implementing a thin conductor is to allow different magnetic fields on the two sides of a zero-thickness interior boundary. This setup is accomplished by adding an Auxiliary Slit feature to the desired boundaries. The feature disconnects the tangential magnetic-field degrees of freedom on the two sides of the boundary, allowing the fields to differ. We therefore need to add an equation relating these degrees of freedom. The fields on the two sides can be accessed using the up() and down() operators.

Adding an Auxiliary Slit feature to any interior boundaries of an interface.

As soon as this slit is introduced, we can define a set of variables using the equations developed in the section above. The screenshot below shows a sample of such a definition. Note that the normal vector is defined in terms of <nxmesh,nymesh,nzmesh>, which is the normal as evaluated on the discrete finite element mesh — so, for example, the z-component of the surface current is:

nxmesh*(up(Hy)-down(Hy))-nymesh*(up(Hx)-down(Hx))

Variables can be defined on the boundaries.

Finally, we need to add a Weak Contribution feature to impose the relationship between the electric field and the surface current. The following contribution is added to the residual:

Es_x*test(K_x)+Es_y*test(K_y)+Es_z*test(K_z)

Keep in mind that the electric field can depend nonlinearly on the surface current via a variable defining the resistivity. For a detailed introduction to weak-form modeling, see this Learning Center series: Modeling with Partial Differential Equations in COMSOL Multiphysics^®.

The Weak Contribution feature implements the relationship between the electric field and surface current.

With the Auxiliary Slit and Weak Contribution features added, along with the appropriate variables defined, the implementation is complete. The selected boundaries can now carry surface current, which can also enter or leave adjacent conductive domains. The surface current can be excited using any of the standard approaches for working with the magnetic field formulation discussed in the previously mentioned blog post on magnetic-field excitations.

An Example From the Field of Superconductivity

As an example, we consider an application from the field of superconductivity: two thin superconducting tapes in free space joined by a relatively large metallic conductor, as illustrated below. The superconducting material itself consists of a very thin layer deposited on a relatively insulating substrate. In this example, we assume that currents in the substrate are negligible and therefore omit the substrate entirely. Current can then flow from one superconducting tape, through the metallic joint, and into the other tape.

The situation under consideration: Two thin superconducting tapes joined with a metallic joint are modeled within a cylindrical free-space region.

A superconductor has a highly nonlinear relationship between electric field and current density. A commonly used constitutive relationship is the power-law expression for the resistivity:

where is the critical electric field criterion, is the critical current density, and is the power-law exponent. In general, can depend on magnetic-field magnitude and orientation, as well as on other fields such as temperature; here, we use the simplified form shown above.

A numerical requirement of this power-law material model within the magnetic field formulation is that the resistivity can never be exactly zero, so we introduce into the above expression a constant, . This constant is chosen to be sufficiently small so that its effect on the solution is negligible.

Similarly, the resistivity of the surrounding insulating domains cannot be infinite. As a practical starting point, we choose the stabilization resistivity so that the skin depth in these domains is approximately one thousand times the characteristic model length, , to get:

This provides a useful balance between numerical conditioning and artificial dissipation. Nevertheless, the resulting artificial losses should be checked after solving. Because of the large contrast in resistivity and the strongly nonlinear constitutive relationship, the problem is well suited for a direct solver. With a compatible GPU, the solution can also benefit from GPU acceleration.

The excitation is introduced via the Tangential Magnetic Field feature. The prescribed field corresponds to a uniform return current along the exterior of the cylindrical modeling domain. The model is solved in the time domain, with a ramped sinusoidal excitation of current such that the local current density in some parts of the tape periodically exceeds the critical current density. The losses within the thin layer, the metallic domain, and the free-space domain are integrated and compared.

The results below show current flowing along one superconducting tape, through the metallic conductor, and into the other tape. The second plot shows the losses over time and illustrates the effect of the nonlinear constitutive relationship. It also shows that artificial losses associated with the stabilization resistivity are small compared with the conductor losses. It’s also useful to vary the stabilization resistivity and confirm that the losses in the physical conductors remain essentially unchanged. The example also demonstrates an important feature of the formulation: Surface current can naturally enter or leave adjacent volumetric conductors.

Current flow (blue arrows) along the thin conductive surfaces and through the adjacent volumetric conductor.

Losses in the thin conductive surfaces, metallic joint, and surrounding free-space domain as a function of time.

Closing Remarks

We have shown how the capabilities of an existing COMSOL Multiphysics^® physics interface can be extended by introducing jump conditions and additional equations on interior boundaries. Although the example presented here is from the area of electromagnetic field modeling, this workflow can be useful in other areas of physics.

It’s worth emphasizing that AI tools can be a great help for this kind of model development work. The priming data for most large language models (LLMs) includes the COMSOL documentation and many technical references, so synthesizing this data into a concrete suggestion for a new feature or formulation can be done in only a few minutes.

The formulation developed here enables thin conductors to be modeled as boundaries within the Magnetic Field Formulation interface and is particularly well suited to superconductors with highly nonlinear resistivity. This approach can be applied to other classes of problems, and if you are interested in learning more, please contact us! The model is available via the button below:

References

1. J.-M. Jin, Theory and Computation of Electromagnetic Fields, 2nd ed., Wiley-IEEE Press, 2015.
2. C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Dover, 2009.

Accelerating Multiphysics Simulation with NVIDIA cuDSS

The NVIDIA CUDA Direct Sparse Solver (cuDSS) is fully integrated into the standard solver framework of COMSOL Multiphysics^®. With this integration, the software can leverage the high memory bandwidth and massive parallel architecture of NVIDIA GPUs for multiphysics analysis. GPU acceleration can substantially speed up the solution of large sparse linear systems, including those that arise repeatedly in nonlinear iterations, implicit time-stepping methods, and continuation methods used in parametric studies. By speeding up the linear system solves at each iteration or time step, NVIDIA cuDSS can improve overall simulation performance.

Nonlinear and time-dependent studies are common in semiconductor modeling, including modeling of applications involving electrokinetics, electromagnetics, and fluid dynamics. GPU acceleration can reduce the runtime of such studies from overnight to just a few hours. For example, for the mechanical polishing benchmark model discussed here, using two GPUs accelerated the runtime from 14–16 hours on the tested CPU configurations to just 2 hours and 9 minutes. Similarly, early benchmark tests on the photoresist spin-coating model presented here showed an approximately 92% lower runtime compared with the tested Intel^® Xeon Gold configuration. Below is a more detailed look at these and other benchmark semiconductor models that show similar speedups thanks to GPU acceleration.

NVIDIA cuDSS Semiconductor Benchmarks in COMSOL

Model: Simulation of the Flow Field and Abrasive Particles in Chemical Mechanical Polishing Equipment

Model geometry of the chemical mechanical polishing (CMP) process, where the wafer is positioned on the polishing head and pressed onto the larger disc of fluid.

CMP is a crucial technique for surface planarization that is widely used in semiconductor fabrication. This model simulates the flow field around the wafer and tracks the motion of abrasive particles carried by the flow. The strongly coupled physics and many degrees of freedom make this a very computationally demanding model. With NVIDIA cuDSS, the model ran more than 6× faster than on the tested CPU configurations, allowing it to be solved and adjusted several times throughout a workday, whereas the CPU-based solvers would need more than a full workday just to run once.

¹ NVIDIA RTX^™ Pro 6000 Blackwell Max-Q 96 GB x 2
² Reported CPU peak memory usage 92 GB

Model: Electrodeposition of a Microconnector Bump with Deforming Geometry in 3D

Model geometry of a microconnector bump, a common feature of chips that are bonded directly onto substrates using bumps and underfill.

This model simulates the evolution of a microconnector bump over time as copper deposits onto an electrode surface. Cupric ions are transported through the electrolyte by convection and diffusion, and the electrode kinetics are described by a concentration-dependent Butler–Volmer expression. The coupling of electrochemistry, fluid flow, and deformed geometry makes for a sizeable computational load, which NVIDIA cuDSS handled more than 2× faster than the conventional CPU solvers.

¹ NVIDIA RTX^™ Pro 6000 Blackwell Max-Q 96 GB
² Reported CPU peak memory usage 8 GB

Model: Copper Deposition in a Trench Using the Phase-Field Method

2D model geometry of a “trench” that causes uneven copper deposition.

During copper deposition, narrow regions can lead to nonuniform plating and the formation of voids, as shown by the blue area in the 2D geometry. The model captures this effect by coupling ion transport and electrode reactions with a phase field formulation based on the Cahn–Hilliard equation, where it tracks the evolving shape of the copper deposit. The benchmark demonstrates that even for a more lightweight model, NVIDIA cuDSS can significantly outperform conventional CPU solvers, in this case running more than 2× faster.

¹ NVIDIA RTX^™ Pro 6000 Blackwell Max-Q 96 GB
² Reported CPU peak memory usage 130 GB

Model: Photoresist Spin Coating

The model geometry of a photoresist coating on top of a substrate.

Photoresist coating is an important process in semiconductor manufacturing where the thickness of the photoresist layer needs to be precisely controlled. Usually, the photoresist layer is thinned by spin coating, where centrifugal force is used to throw excess photoresist away from the wafer surface until a specific thickness is achieved. This model simulates the process using swirl flow to calculate the film thickness at different viscosities and rotational speeds. NVIDIA cuDSS substantially reduced the runtime for this model, achieving a speedup of up to 13× compared with the CPU-based solvers.

¹ NVIDIA RTX^™ Pro 6000 Blackwell Max-Q 96 GB
² Reported CPU peak memory usage 24 GB

Model: Electrostatic Discharge Test of a PCB

A PCB mounted on a copper ground plane being subjected to electrostatic discharge (ESD) events.

PCBs are subject to ESD events that can cause malfunctions or logical errors. ESD testing helps engineers assess whether electronic assemblies can withstand specified discharge conditions. Numerical simulation can identify ESD problems during the design stage while providing detailed electromagnetic field distributions that are difficult or even impossible to measure accurately in the real world. ESD currents are highly nonlinear, which makes the ESD testing model a good candidate for GPU acceleration. In this case, NVIDIA cuDSS achieved a speedup of about 3× compared to the tested CPU configurations.

¹ NVIDIA RTX^™ Pro 6000 Blackwell Max-Q 96 GB x 2
² Reported CPU peak memory usage 22 GB

Model: Residual Stress and Warpage Due to the Curing of the Underfill in a Chip Package

Model geometry of a flip chip (FC), a widely used form of semiconductor packaging where the chip is directly bonded to the substrate using bumps and a polymer underfill that fills the gap between bumps and substrate.

FC packaging allows for high interconnection density as well as a reduction in package size, but the different thermal expansion properties of the materials involved may cause warpage of the package during underfill curing. This model simulates the warpage of the chip package throughout the curing process as the package is exposed to different thermal conditions. NVIDIA cuDSS reduced the runtime substantially, with a speedup of more than 6× compared to the tested CPU configurations.

¹ NVIDIA RTX^™ Pro 6000 Blackwell Max-Q 96 GB
² Reported CPU peak memory usage 62 GB

Learn More About GPU Acceleration in COMSOL^®

By accelerating multiphysics simulations with NVIDIA cuDSS, the COMSOL^® software helps engineers evaluate more process and design options in less time — from copper deposition and photoresist coating to package warpage and polishing equipment. Faster simulation supports shorter development cycles and more informed engineering decisions before physical testing. To hear more about the capabilities of NVIDIA-accelerated computing support for semiconductor applications, reach out to your local COMSOL representative.

Table of Contents

• Selecting an Electrochemical Interface
• Primary, Secondary, and Tertiary Current Distribution Theory
  • Primary Current Distribution Theory
  • Secondary Current Distribution Theory
  • Tertiary Current Distribution Theory
• Physics Interfaces in the Electrochemistry Branch
  • Primary Current Distribution Interface
  • Electrode, Shell Interface
  • Secondary Current Distribution Interface
  • Tertiary Current Distribution, Nernst–Planck Interface
  • Aqueous Electrolyte Transport Interface
  • Electrophoretic Transport Interface
  • Concentrated Electrolyte Transport Interface
• Summary of Physics Interface Assumptions

Selecting an Electrochemical Interface

The flow chart below can be used to provide guidance for selecting an interface. By considering which physical phenomena need to be resolved, including electrolyte mass transport, migration, electrode kinetics, and dissociation equilibria, as well as the assumptions that can be made about the electrolyte and electrode geometry, a suitable physics interface can be selected as a starting point. For each interface, the underlying theory is discussed in more detail, explaining the assumptions behind it.

Flowchart for selecting an electrochemical interface.

In the following sections, the current distribution underlying each physics interface is introduced, along with representative examples of how to use the interface.

Primary, Secondary, and Tertiary Current Distribution Theory

Primary, secondary, and tertiary current distribution represent increasing levels of physical detail. Current distribution theory starts from charge conservation in the electronic conducting and electrolyte phases:

where and are ​the current densities in the electrode and electrolyte phases, and and are volumetric source or sink terms.

In primary and secondary current distributions, current conduction in both the electrode and electrolyte phases is described by Ohm’s law. The electrolyte composition is assumed to be uniform, and a constant electrolyte conductivity can be used. Current densities in the electrode and electrolyte phases are then written as:

where and are the conductivities of the electrode and electrolyte, and and are potentials in the electrode and electrolyte, respectively.

Primary Current Distribution Theory

In the primary current distribution, only ohmic losses in the electrochemical cell are considered, while electrode kinetics are neglected. Thus, no activation overpotential is considered, and the potential difference between the electrode and electrolyte is equal to the equilibrium potential, which can be expressed as:

Secondary Current Distribution Theory

In the secondary current distribution, electrode kinetics is included in addition to the ohmic effect. In this case, the local current density at the electrode–electrolyte interface is described as a function of the overpotential, which is the deviation of the electrode and electrolyte potential difference from the equilibrium potential:

The local current density for an electrochemical reaction can be described using an appropriate kinetic expression, such as the Butler–Volmer equation:

where is the exchange current density, and are charge transfer coefficients, is the universal gas constant, is Faraday’s constant, and is the temperature.

Tertiary Current Distribution Theory

In the tertiary current distribution, mass transport in the electrolyte is considered, and concentration gradients are no longer assumed to be negligible. The mass balance equation for each species i in the dilute electrolyte is given by the mass conservation equation:

where , , and denote the total flux, concentration, and reaction source of species i, respectively. In the electrolyte, the total flux of species i can be expressed as:

where is the charge number, is the diffusion coefficient, is the ​mobility, and is the velocity field. The first, second, and third terms on the right-hand side represent diffusion, convection, and migration, respectively.

The electrolyte current density is then obtained by summing the contributions from all ionic species fluxes:

As a result, the electrolyte potential is solved for together with the species concentrations, rather than from Ohm’s law alone. This allows the effects of concentration gradients on ionic species transport to be accounted for in the current distribution. In tertiary current distribution, the electrolyte current is therefore influenced not only by conductivity and electrode kinetics, but also by changes in electrolyte composition.

Physics Interfaces in the Electrochemistry Branch

The assumptions underlying primary, secondary, and tertiary current distributions were introduced in the previous section. The interfaces in the Electrochemistry branch reflects the same progression, with different assumptions about current conduction, electrode kinetics, and mass transport in the electrolyte. In this section, these interfaces are introduced together with representative examples using each interface.

Physics interfaces in the Electrochemistry branch of the Select Physics window.

Primary Current Distribution Interface

The Primary Current Distribution interface is based on the theory of primary current distribution. It is used when the current distribution is governed mainly by ohmic effects and the cell geometry, and the influence of activation losses is neglected. The following model examples use this interface:

• Electrocoating of a Car Door
• Electrochemical Machining of a Microbore
• Primary Current Distribution in a Lead–Acid Battery Grid Electrode

Electrode, Shell Interface

The Electrode, Shell interface can be used for thin electrodes when the potential variations are in the direction normal to the electrode surface. Rather than explicitly resolving electrode thickness, the domain is replaced by a boundary equation that describes tangential current conduction. Current conduction in the shell electrode is described by Ohm’s law together with charge conservation. The following model example uses this interface:

• Electrodeposition on a Resistive Patterned Wafer

Secondary Current Distribution Interface

The Secondary Current Distribution interface extends the Primary Current Distribution interface by including activation overpotentials through electrode kinetics. It can be used when electrode kinetics affect the current distribution, while mass transport in the electrolyte can be neglected. The following model examples use this interface:

• Current Distribution in a Chlor–Alkali Membrane Cell
• Decorative Plating
• Electroplating of a Printed Circuit Board
• Aluminum Anodization
• Anode Film Resistance Effect on Cathodic Corrosion Protection
• Galvanic Corrosion of a Magnesium Alloy in Contact with Steel
• Jelly Roll
• Jelly Roll Using a Flattened Geometry

Tertiary Current Distribution, Nernst–Planck Interface

As discussed, tertiary current distribution is used when mass transport in the electrolyte is no longer negligible and when concentration gradients influence the current distribution. The Tertiary Current Distribution, Nernst–Planck interface accounts for these effects by coupling ionic species transport to electrolyte charge conservation. The interface includes five electrolyte charge-conservation models that describe coupled charge and mass transport in the electrolyte.

Charge conservation models.

Charge Conservation Model, Electroneutrality

The Electroneutrality option is used for systems with significant concentration gradients of ionic species. Under this assumption, the concentrations of the major ionic species are resolved explicitly, while local electroneutrality is imposed according to:

This is the general tertiary formulation for systems in which the major charge-carrying ions are included explicitly, and electroneutrality is assumed in the electrolyte. By combining the ionic flux expressions with the electroneutrality condition, the electrolyte current density can be expressed as:

where electrolyte conductivity can be defined as:

It is commonly used when concentration gradients are significant, and the main ionic species must be resolved directly. The following model examples use this option:

• Electrochemical Treatment of Tumors
• Capacitive Deionization of Saline Water
• Desalination in an Electrodialysis Cell (Ion-exchange membrane)
• Corrosion Inhibitor Transport
• Lithium Sulfur
• Vanadium Redox Flow Battery
• Electrochemical capacitor with porous electrode

Charge Conservation Model, Water-Based with Electroneutrality

The Water-based with electroneutrality option extends the electroneutrality option by incorporating the water auto-ionization equilibrium reaction for aqueous electrolyte. It accounts for proton and hydroxide ion transport, allowing the local pH to be calculated from proton concentration. The following model examples use this option:

• Alloy Deposition
• Parasitic reaction in an electrochemical capacitor

Charge Conservation Model, Supporting Electrolyte

The Supporting electrolyte option can be used when the system contains a strong or well-mixed supporting electrolyte so that most of the concentration of charge-carrying species does not vary significantly. In this case, the electrolyte conductivity may also be treated as constant. The following model examples use this option:

• Cathodic Protection of Steel in Reinforced Concrete
• Zinc–Bromine Redox Flow Battery

An equivalent implementation of the Tertiary, Supporting electrolyte option can be obtained by coupling the Secondary Current Distribution interface with the Transport of Diluted Species interface from the Chemical Species Transport branch. The following model examples use this implementation:

• Wire electrode
• Orange battery
• Rotating cylinder hull cell
• Galvanized nail

Charge Conservation Model, Poisson

The Poisson option can be used when local electroneutrality is not assumed, and the potential in the electrolyte is obtained from Poisson’s equation. It can be used in cases where charge-separation effects are of interest, typically within a few nanometers of an electrode surface. The following example models use this option:

• Diffuse Double Layer with Charge Transfer
• Ion-Exchange Membranes and Donnan Potentials

Charge Conservation Model, Electroanalysis

The Electroanalysis option can be used for electroanalytical problems in electrolytes containing a large excess of inert supporting electrolyte. Under this assumption, the conductivity is high enough that ohmic losses can be neglected. In addition, the electrolyte potential is not solved, and the migration effect is neglected. The transport of ionic species is described mainly by diffusion, with convection included when needed. The following model examples use this option:

• Cyclic Voltammetry at a Macroelectrode in 1D
• Electrochemical Impedance Spectroscopy
• Voltammetry at a Microdisk Electrode
• Alternating Current-Induced Corrosion
• Adsorption–Desorption Voltammetry
• Glucose Sensor
• Thin Layer Chronoamperometry

The Electroanalysis option can be either chosen from the Tertiary Current Distribution, Nernst–Planck interface or directly chosen as a physics interface from the Electrochemistry branch when adding physics.

Aqueous Electrolyte Transport Interface

The Aqueous Electrolyte Transport interface can be used for aqueous systems containing weak acids, weak bases, ampholytes, and complex species. It is based on the same assumptions as the Tertiary Current Distribution, Nernst–Planck interface when the Water-based with electroneutrality option is selected. Due to its more efficient handling of multiple equilibrium reactions and easier model setup, the Aqueous Electrolyte Transport interface can be a more convenient choice in some cases than the Tertiary Current Distribution, Nernst–Planck interface, especially for applications such as atmospheric corrosion, mechanistic corrosion modeling, biological electrochemistry, and electrochemical sensor modeling. The following model examples use this interface:

• Accelerated Corrosion Test of a Scratched Galvanized Steel Sample
• Atmospheric Corrosion with Mass Transport
• Corrosion Under an Evans Droplet
• Crevice Corrosion of Iron in an Acetic Acid/Sodium Acetate Solution
• Carbon Dioxide Corrosion in Steel Pipes
• Pitting Corrosion in Stainless Steel

Electrophoretic Transport Interface

The Electrophoretic Transport interface can be used for water-based systems in which species transport is driven by electrolyte potential gradients, and electrode kinetics are not considered. This interface supports weak and strong acids and bases, ampholytes, proteins, and uncharged species. Transport is described by the Nernst–Planck equations together with electroneutrality, dissociation equilibria, and water auto-ionization. Compared with the Aqueous Electrolyte Transport interface, it is tailored to electrophoresis applications. The following model examples use this interface:

• Zone Electrophoresis
• Isoelectric Focusing
• Isoelectric Separation

This interface can be added from the Chemical Species Transport branch.

Concentrated Electrolyte Transport Interface

The electrolyte transport formulations discussed so far are based on dilute-solution theory. The Concentrated Electrolyte Transport interface can be used when the electrolyte needs to be modeled as a concentrated solution, with transport governed by the Maxwell–Stefan diffusion equation assuming local electroneutrality. In contrast to the Nernst–Planck equations, the concentrated solution theory does not assume that the electrolyte species are diluted in a neutral solvent of constant concentration. Typical electrolytes that can be modeled include ionic liquids, molten salts, and highly concentrated solutions. More theory is discussed in the blog post Introduction to the Concentrated Electrolyte Transport Interface. The following model example uses this interface:

• Molten Carbonate Transport

The Electrochemistry branch includes interfaces at different levels of electrochemical modeling, from generic current-distribution interfaces to specialized interfaces for different types of electrochemical applications.

Summary of Physics Interface Assumptions

The interfaces discussed above use different assumptions for current distribution, electrolyte transport, and electrolyte description. The table below summarizes the main assumptions and formulations associated with each interface. A check mark indicates that the corresponding physical effect, assumption, or transport formulation is included in the interface. A blank indicates that it is not part of the default formulation or is not applicable.

Specification chart of generic electrochemical interfaces.

The interfaces discussed in this blog post provide a starting point for building electrochemical models in COMSOL Multiphysics^®. They span different assumptions, levels of detail, and application areas. The interface that is chosen helps define the scope of the model from the beginning and determines how current distribution, electrode kinetics, and species transport are represented. It also provides the basis for extending the model to include additional physical effects or multiphysics couplings.

With Formula E debuting the revolutionary GEN4 race car for its 2026–2027 season and COMSOL Multiphysics^® v2027 set to be released later this fall, the agreement comes at a time of major technological advancement for both organizations.

What will this partnership mean for the future of electric racing and simulation-led development? Let’s take a look under the hood at each industry leader’s latest innovations to find out.

Formula E’s New Era: GEN4

Formula E’s 2026–27 season (Season 13) ushers in GEN4, setting new benchmarks for speed and innovation in motorsport. Capable of hitting speeds of 208 mph (335 km/h) and reaching 124 mph (200 km/h) in just 4.4 seconds, the GEN4 car represents a huge technical leap for the world’s leading all-electric motorsport series. 

GEN4 also sets a new global benchmark as the world’s most sustainable race car, becoming the first to be 100% recyclable or reusable, with a minimum of 20% of recycled materials used for the tires, batteries, chassis, and front powertrain kit (FPK).

The GEN4 Formula E car seen during the GEN4 Unleashed launch event held April 21–22, 2026, at Circuit Paul Ricard, Le Castellet, France. (Photo by Simon Galloway/LAT Images)

A Performance Jump

The leap to GEN4 starts with a major performance improvement over previous generations. Redefining high-performance electric motorsport, the GEN4 race car will deliver 50% more power in race mode than the GEN3 Evo. With a power output of 600 kW in ATTACK MODE, GEN4 represents a 71% increase in base power output compared to GEN3 Evo.

GEN4 cars also feature permanently active all-wheel drive (AWD), which improves traction and torque management and helps them accelerate from 0 to 100 km/h (0 to 62 mph) in around 1.8 seconds. Teams, drivers, and fans will see a marked difference in performance between these cars and any from previous seasons.

Energy Management and Efficiency Remain Key

The GEN4 championship fight will continue Formula E’s traditional challenge to competitors: to find the most energy-efficient way to be the fastest over the course of a lap, a race, and a season.

The race energy available to each GEN4 car has increased by 43% compared to the GEN3 Evo (to 55kWh). However, this additional energy capacity comes alongside the increase in power and performance, making how efficiently teams deploy and recover energy more important than ever. These changes will lead to faster races without compromising on Formula E’s emphasis on energy efficiency.

The GEN4 car will showcase its speed across two major race formats: the E-Prix and the E-PrixUnleashed. The former maintains the racing series’ customary focus, challenging competitors to manage their energy better than the rest of the field. However, the latter will show off the car being driven at its ultimate limit, highlighting the raw speed of the GEN4 as it laps some of the world’s most recognizable circuits.

COMSOL Multiphysics^® v2027: Innovating with Speed, Efficiency, and Accuracy

With each generation of the ABB FIA Formula E World Championship, greater performance is unlocked by solving complex engineering challenges. Similarly, COMSOL Multiphysics^® was designed to help users tackle a wide variety of engineering challenges quickly and accurately. The software provides built-in functionality for analyzing how the interactions between different physical phenomena affect the way individual parts and the system perform. This fall, the software will become even more capable with its own leap forward: v2027. The new version will introduce faster computation, system-level modeling, and expanded tools for building and deploying simulations.

Unveiling the Electric Motor Module

The Electric Motor Module, a new add-on product to the software, will provide a streamlined environment for modeling electric motors and generators in a wide variety of commercial or industrial applications, including road cars. Using the module’s dedicated functionality for evaluating performance, operating ranges, losses, material effects, and faults in both 2D and 3D, engineers can optimize the performance, efficiency, and material utilization of machine designs. Just as GEN4 capabilities like permanently active AWD and greater energy capacity have already begun to do, electric motor designs developed using this new module will pave the way for future commercial vehicle technologies.

A 3D model of an axial-flux motor for an electric vehicle, illustrating the detailed electric machine modeling capabilities of the Electric Motor Module.

Connecting Models Through COMSOL Systems™

COMSOL Multiphysics^® v2027 will also expand the scale at which engineers can approach complex systems with the introduction of COMSOL Systems™. This add-on product provides a new modeling environment for analyzing complete systems while retaining access to detailed multiphysics modeling wherever local physical effects need to be resolved.

Using a library of nearly 5000 ready-to-use components, engineers can build electrical circuit diagrams, process flowsheets, and custom equation-based system models. These system-level representations can then be directly connected to 1D, 2D, or 3D multiphysics models created in COMSOL Multiphysics^®.

An electrical circuit model created in COMSOL Systems™ coupled to a detailed 3D multiphysics model created in COMSOL Multiphysics^®, enabling circuit-level and component-level physics to be analyzed together.

Such technologies could also be used in the context of an electric race car, like Formula E’s GEN4. Important aspects of the electric powertrain could be modeled at the system level using COMSOL Systems™ while incorporating detailed multiphysics models of components such as the motor, battery, or regenerative braking system. This multilevel approach would enable them to zoom in on, for example, the electric motor to capture the effects of electromagnetic losses, heat generation, and temperature on motor operation as needed while efficiently analyzing the behavior of the vehicle as a whole.

Engineering a Faster Future, Together

Being at the forefront of electric motorsport competition and engineering a successful product both require an emphasis on optimization and efficiency. This synergistic intention makes the partnership between Formula E and COMSOL mutually beneficial as each party reaches pivotal milestones in its technological journey.

With GEN4, Formula E is introducing its fastest car yet while maintaining the focus on efficiency and energy management that has defined the Championship. At the same time, COMSOL Multiphysics^® version 2027 will offer engineers powerful new capabilities for modeling electric machines, analyzing interconnected systems, and investigating the multiphysics interactions that influence real-world performance across many industries.

Separately, these two technology leaders are driving innovation. Through this partnership, they can accelerate it.

Accelerating Electric Racing with Simulation

To learn more about COMSOL Multiphysics^®, the upcoming release of v2027, and how the software is used in motorsport, check out the links below:

• Modeling Drag Reduction Systems for Motorsports Using CFD (blog post)
• COMSOL in Sports and Racing (ebook)
• System-Level Modeling and Agentic AI Take the Spotlight in COMSOL Multiphysics^® Version 2027 (press release)

Building the COMSOL^® Model

To get a better idea of the interference these sensors may introduce, a collaborative team from GEPEA–Oniris (Process Engineering–Environment–Agri–Food Joint Research Unit, National College of Veterinary Medicine, Food Science and Engineering, Nantes, France) and CTCPA (agri-food technical center in France) ran a study that evaluated these in-package microwave heating processes by simulating the effects of different sensor positions within a single-mode microwave applicator. The study focused on a lab-scale 915-MHz microwave heating process for 300 g of mashed potato in a polypropylene tray. The 3D model geometry of the microwave equipment was built in the COMSOL Multiphysics^® software, while the experimentally measured geometry of the polypropylene tray was imported to COMSOL^® from a custom-made STL file.

Figure 1. Description of the 915-MHz single-mode microwave apparatus.

The model the team built was robust, taking into account all of the geometrical details of the microwave apparatus. It included the complete design of a waveguide transition, an impedance matching element (i.e., an iris), a 915-MHz single-mode microwave applicator, and a sliding short circuit. The microwave input power source was modeled with a coaxial port at the microwave antenna. All surfaces of the waveguide, including the copper sliding short circuit, were considered as perfect electric conductors. The dielectric properties of the mashed potato were experimentally measured and assumed to correspond to a constant moisture content of 76% on a wet basis.

To avoid microwave interference with metallic parts, the wireless sensor was surrounded by an ogive-shaped 316-L stainless steel microwave shield and placed at the upper surface of the mashed potato sample. Additionally, a polypropylene film was heat-sealed at the top surface of the polypropylene tray to close the system. Each of these components came together to create a numerical model that helped the team derive accurate insights into the best in-package sensor position to limit interactions with the electromagnetic field.

Figure 2. Description of the polypropylene tray filled with a mashed potato sample. MW = microwave, PP = polypropylene, and PTFE = polytetrafluoroethylene.

Simulation-Based Insights

The team used the RF Module, an add-on to COMSOL Multiphysics^®, to predict the electromagnetic field in the microwave cavity relative to the microwave-absorbed power within the sample of mashed potato. For the solver settings, boundary mode analysis was used as a first step to compute the mode at port 2 (see Figure 1), which comes out of the waveguide. The second step was using a Frequency-Stationary study at 915 MHz.

Without a sensor inside the tray, the simulation revealed microwave-reflected power of 24% (S₁₁ = -6.16 dB), with the optimal sliding short circuit position at 790 mm. S₁₁, the input port reflection coefficient (i.e., the measure of the port’s mismatch), was then evaluated as a function of varying microwave shield positions following different z-axis rotations, from alpha = 0 to 180^°, with 30^° increments.

Figure 3. Local microwave-absorbed power within the mashed potato as a function of four sensor positions.

The results showed that the microwave shield position needed to be aligned to the microwave propagation direction (alpha = 0^°) to obtain a similar S₁₁ parameter to the food tray without the sensor (S₁₁ ≈ -6 dB, RF ≈ 25%). When rotated counterclockwise by 90^°, so the microwave shield position faced the same direction as the transverse electric plane, the S₁₁ parameter dropped to -7.8 dB. Furthermore, by changing the rotation angle of the microwave shield from 90 to 180^°, simulation results indicated S₁₁ variations from -7.8 to -4.7 dB. This rotation from 90 to 180° also resulted in the simulated microwave reflected power increasing by a factor of 2, with 17% microwave reflected power at 90° and 33% at the 180° rotation. These results suggest that the sensor positioning is crucial to take into account, as these reflected power variations can impact overall microwave power efficiency.

A portable vector network analyzer (VNA) was used in the experimental validation of the numerical model. The VNA was connected at port 1 to evaluate the microwave-reflected power from the load according to different operating conditions. The S₁₁ parameter was used to account for the microwave-reflected power at the antenna relative to the incident microwave power. Overall, the simulation results were in agreement with the experimental measurements performed with the VNA. Discrepancies between the experimental and simulation results can be attributed to slight variations in sensor positioning due to the difficulty of precisely replicating the uneven surface of the mashed potato.

Figure 4. Reflection factor following different rotation angles of the microwave shield.

An Important Discovery

This study led to an important discovery: Aligning the sensor position with the microwave propagation direction resulted in a similar microwave reflection coefficient as tests run without a sensor. However, a significant reduction of the microwave-reflected power is observed if the microwave shield is positioned within the transverse electric plane of the waveguide.

Testing showed that the model is also able to predict the high local microwave power densities at the interface between the sensor and the food sample, which could further elucidate the best position for a temperature and pressure (T-and-P) sensor, though this would require more experimental validation and accounting for potential temperature-dependent dielectric properties of the food, as well as their coupled impact on the electric field distribution.

Further Learning

For more information on this research, read the GEPEA-Oniris and CTCPA team’s full paper, which won a Best Paper award at the COMSOL Conference 2025 Amsterdam! The paper describes the team’s full modeling approach and results.

The vibration of structures is affected by the surrounding media. Air, one of these media, is light and compressible and thus can often be neglected when a structural engineer evaluates the vibration of a design. Simulations might be performed without a physical representation of the air — as if the structure were vibrating in a vacuum without any interference with a surrounding medium. Conversely, when a structure is submerged in water, with one thousand times the density of air and a higher speed of sound (stiffness), the water’s effect on the vibration characteristics is often noticeable.

Accounting for the Effects of the Surrounding Medium

Traditionally, this effect has been denoted “added mass”, derived from the following equation used to describe it.

In the equation, is the mass of the structure, and is the added mass. The value of the added mass has often been empirically derived by experiments where the eigenfrequencies of the structure have been measured in air and water, respectively. Eigenfrequencies, also known as natural or resonance frequencies, are the frequencies at which a structure naturally reinforces vibrations. The difference in eigenfrequencies is then used to calculate the added mass.

Unfortunately, this approach has led to the misunderstanding that the effect of the surrounding medium is an added mass that vibrates together with the structure, whereas actually, the surrounding media acts with a pressure on the surface of the moving object, and this pressure is highest when the surface speed is highest. Therefore, the pressure is not constant but changes in amplitude and direction depending on the direction of the moving surface.

The true influence of the surrounding fluid on the vibration response can be captured using COMSOL Multiphysics^® and a coupled acoustic–structure interaction (ASI) modeling approach. Similar multiphysics simulations have been described in the blog post “What is the Best Way to Analyze Fuel Tank Vibration” and also in our blog posts on muffler design, such as “Evaluating the Effect of Shell Thickness on Muffler Performance“.

Simulating the Surrounding Media

The focus of this blog post is highlighting the effect of the surrounding media and clearly describing an alternative to “added mass”. The physical example used in this blog is a Kaplan turbine blade in a water domain. The effect on its frequency response from a varied gap to the pipe wall is studied. The figure below shows the full Kaplan turbine geometry and the considered gap. Note that pressure equalization at the leading and trailing edge of the blade is always possible, irrespective of the varied gap size.

Kaplan turbine geometry with arrows showing the gap to the pipe wall.

Acoustic–Structure Interaction

The fundamental part of acoustic–structure interaction is the multiphysics coupling applied to the boundary between the fluid and the structural domain. The pressure degree of freedom (DOF) in the fluid is coupled to the three displacement degrees of freedom (DOFs) in the structure (in the x, y, and z directions), as described below.

The conditions for an external boundary are

in which is the structural acceleration (applied to the fluid), is the total acoustic pressure, is the density of the fluid (the subscript c indicates that this can be a complex-valued quantity used to model losses in the fluid ), is an optional acoustic dipole domain source term, is the surface normal, and is the (acoustic) load (force per unit area) on the structure.

The acoustic pressure is applied to the structural faces, and the acceleration of the surface is coupled to the acoustic domain. To enable the modeling of acoustic–structure interaction in COMSOL Multiphysics^®, both the Pressure Acoustics, Frequency Domain and Solid Mechanics interfaces need to be added to the model. Furthermore, the multiphysics coupling described above needs to be enabled and added (under the Multiphysics node) to all boundaries between structural and fluid volumes. Here is how it looks in the Model Builder tree:

Building a vibroacoustic simulation.

As the problem is solved in the frequency domain, the model is assumed linear, and only small structural deformations can be modeled. Losses in the fluid can be accounted for in several ways using a fluid model or even more advanced boundary conditions such as the Thermoviscous Boundary Layer Impedance feature. Damping due to vorticity in a moving fluid, and even the effect of a moving fluid, requires solving a more advanced acrostic formulation such as the linearized Navier–Stokes equations, also available in COMSOL^®.

Example Using a Kaplan Turbine Blade

In the following, the influence of the inner radius of a pipe has on the vibration response on the general Kaplan turbine, is investigated. With the functionality in COMSOL Multiphysics^®, it is possible to apply periodic boundary conditions and only model one blade while capturing the effect of overlapping blades. The water volume is created by offsetting the circular sections at the top and bottom so that their connecting sides pass between the blades. On those sides, the periodic boundary condition is applied (shown as yellow in the figure).

Water domain, with periodic boundary conditions in yellow.

The turbine blade was modeled as linear elastic with Young’s modulus 210 GPa, Poisson’s ratio 0.3, and density 7850 kg/m³. The water was assigned the properties of speed of sound at 1500 m/s and density at 1000 kg/m³.

For the model with a 10 mm gap to the pipe wall, 310,000 quadratic tetrahedral elements made up the mesh with refinements around the blade. The gap had at least 5 elements for all cases. The entire mesh was continuous, and no contact exists in the model.

All outer faces were set to the Sound Hard Boundary (Wall) boundary condition, which sets the normal component of the acceleration/velocity to zero. The condition is simply the natural Neumann condition for the total pressure, i.e., the normal derivative of the pressure is zero on the boundary

For the eigenfrequency analysis, the hub mount face of the blade was fully constrained

When determining the response to excitation in the frequency domain (the frequency response function, or FRF), an axial rotation of rad/s² was applied instead of the fully constrained boundary condition. That gives an acceleration of 1 rad/s², regardless of the frequency.

Results of the Acoustic–Structure Interaction Analysis

First, it was determined that the height of the water volume did not influence the results when the vertical distance between the blade and water boundary exceeded 10 cm above and under the blade. Therefore, this size of the water domain will be used throughout the study.

Second, the influence of the surrounding fluid (water or air) on vibration properties can be determined with the ASI approach. In the table below, the first six eigenfrequencies in water, with a gap large enough to not influence the turbine blade, are compared to the first six eigenfrequencies in free air. The table shows that water has a large influence on the eigenfrequencies, even when the gap has no effect.

Finally, the influence of a varied gap size on eigenfrequencies can be determined. As the gap increases, eigenfrequencies do as well, up to a gap of about 150 mm. Here, eigenfrequencies converge to a constant value, regardless of a continued increase in gap size. The first six real eigenfrequencies for each gap size were derived and are presented below.

The first six real eigenfrequencies as functions of the gap.

The decreasing eigenfrequencies for smaller gaps can be explained by the increase in resistance for pressure equalization around the blade. It is not to be confused with damping from viscous effects, nor should it be simplified to an additional mass. Instead, we see how the water, modeled explicitly, affects the eigenfrequency of the system consisting of the turbine blade and water domain coupled together.

The corresponding mode shape for the fifth eigenfrequency for a gap of 149 mm is shown below, left. In COMSOL Multiphysics^®, it is possible to plot the entire structure — not just the blade that is used in the calculation — giving a better overall picture of the state of the structure, as seen below, right.

Mode shape 5 with a 149 mm gap, showing the blade used in the calculation (left) and the entire structure (right).

3D illustration of a Kaplan turbine propeller.

When the structure is instead excited by adding a frequency-independent rotation on the hub mount face of the blade, it results in the frequency-domain response presented in the figure below. Note that the frequency domain is swept in steps of 1 Hz and therefore the amplitude of the eigenfrequencies might vary depending on how close to the exact value is reached in this simulation. Regardless, peak response amplitude is not of the highest interest since it is most often desired to avoid these frequencies with some margin of safety in product development.

Acceleration response in the frequency domain.

Using functionality in the COMSOL Multiphysics^® software, this example study has shown how modeling acoustic–structure interaction can be used to evaluate the frequency response for a structure submerged in fluid. This is an important step in the product development process for all products where the vibration response is a design criterion in the same way for underwater products as for products in air, where frequency analysis has been utilized for a long time.

In this example, it was made clear how the resistance of the flow of the fluid influences the frequency response of the structure in addition to the difference of submerging in water instead of in air. The eigenfrequency decreased by 57%, from 165 Hz to 71 Hz, by accounting for the surrounding water and minimizing the gap (to 3 mm).

Dynamic Loads on the Blade

A useful feature in COMSOL^® is the Rotating Frame feature, which simplifies rotor dynamics by keeping your mesh stationary. By directly applying centrifugal, Coriolis, and Euler forces as body loads, it accurately simulates high-speed spin effects such as stress stiffening. This feature was not used in this example due to the low rotational speeds for this general application, but it is a powerful tool if the application requires it.

About the Guest Author

Tobias Jonsson is a structural mechanics specialist and consultant at Lightness by Design in Stockholm, Sweden. With an MSc from KTH Royal Institute of Technology, Tobias uses numerical simulation to help companies with simulation-driven product development. As part of a COMSOL Certified Consultancy, he is dedicated to integrating advanced simulation into the core of the product development process.

The consultants at Lightness by Design have contributed to the field of acoustic–structure interaction by coauthoring several research papers published in scientific journals. Furthermore, they have used their expertise in acoustic–structure interaction approaches to simulation and vibration design to help multiple clients in the product development process.

The Resonator Essentials: Physics, Admittance, Damping, and Electromechanical Coupling

Before we dive into the circuit formalism, I will set the scene from a physics standpoint. A standard bulk acoustic wave (BAW) or surface acoustic wave (SAW) resonator uses the piezoeffect to trap acoustic energy at natural frequencies inside a microscale cavity. To analyze this behavior from the electrical point of view, we track terminal admittance, , where its:

• Real part, conductance, , represents energy transfer
• Imaginary part, susceptance, , tracks energy storage

This electromechanical interaction yields two primary frequencies. At the series resonance frequency , standing acoustic waves constructively interfere. Internal mechanical impedance drops to a minimum, causing a sharp peak in electrical conductance . The antiresonance or parallel resonance frequency is located slightly higher. Here, the piezoelectrically generated charge from mechanical strain cancels out the electrostatic plate charge on the electrodes. The device draws minimal net current, driving the total admittance magnitude to a minimum.

The gap between and defines the effective electromechanical coupling coefficient , which is a measure of the energy conversion efficiency: . A larger frequency split provides a higher , directly translating to a wider achievable filter bandwidth. The exact definition of the effective electromechanical coupling coefficient depends on the standard or approximation you follow. The formula may deviate a bit between BAW and SAW applications or if your conditions satisfy the small-coupling approximation.

In an ideal lossless resonator, the conductance peak at would approach infinity, and the admittance valley at would drop to absolute zero. Real-world piezoelectric devices, however, experience energy dissipation due to internal material damping, both mechanical and dielectric, as well as other mechanisms, such as anchor leakage into the substrate. This total dissipation dictates the device’s Q factor, . On the admittance curve, higher material damping directly suppresses and flattens the response. It rounds off the sharp resonance peaks and valleys, lowering the maximum conductance at and raising the minimum admittance floor at . For RF filter designers, minimizing this damping to achieve a high Q factor is essential to ensure low insertion loss and steep filter skirts.

A convenient admittance magnitude |Y| curve sketch with the series resonance and parallel antiresonance points mapped out.

Representing a Piezoelectric Resonator in a FEM Manner

If you were to create a representative FEM model of a piezoelectric resonator, it would be straightforward in COMSOL Multiphysics^® and reveal the powerful multiphysics capabilities of the software. The steps are as follows: Start with a Piezoelectricity interface that couples an Electrostatics interface with a Solid Mechanics interface. For Solid Mechanics, select the Piezoelectric Material material model and specify the relevant lossy behavior through its subnodes, such as Mechanical Damping and Dielectric Loss.

To run a Frequency Domain study, you might want to drive your device via a Terminal condition, added in the Electrostatics interface. Among other benefits, this feature automatically derives the terminal admittance, making it readily available for plotting and further results evaluation. The results evaluation includes extracting the mentioned resonance frequencies and Q factors. Alternatively, you may opt for an Eigenfrequency analysis, which also helps retrieve resonance frequencies and Q factors.

Bridging Finite Elements and Lumped Parameters

The beauty of a lumped-parameter circuit model lies in its simplicity, as solving a network of a few resistors, inductors, and capacitors takes fractions of a second but still provides you with the sought admittance curve. The electromechanical equivalent framework, which is central to this blog post, traces its lineage back to Stephen Butterworth’s filter designs and Walter Guyton Van Dyke’s early-20^th-century work mapping the behavior of quartz crystals. The “modified” modern version adds parasitic elements to fit the stringent demands of high-frequency MEMS resonators.

Modern RF MEMS resonators demand an engineering synergy between fast circuit models and modern numerical solvers. A circuit model is only as good as its input parameters. You cannot simply guess the equivalent capacitance of a complex, layered SAW interdigital finger topology or the motional resistance of a novel piezoelectric thin film. By using COMSOL Multiphysics^® to simulate the true geometric and material complexities of your device, you can extract the exact lumped parameters needed to populate your mBVD circuit. Altogether, the FEM provides the ground-truth physics, while the circuit representation allows for the system-level speed.

Note: The classic mBVD model configuration to be discussed here is chiefly designed for one-port resonators featuring a single signal terminal and a ground terminal. If your design boils down to a multiport configuration or the like, you will need to expand this framework into a more elaborate one.

Similarly, for the sake of clarity, we will focus on one particular resonance. You will need to add an extra parallel motional branch for every resonance mode you need to consider.

Breaking Down the mBVD Circuit

The conventional mBVD model captures a standard piezolelectric resonator behavior by splitting the electrical behavior after the initial series routing into two parallel paths: a static branch representing the physical dielectric layer stackup and a motional branch representing the electromechanical acoustic resonance. A circuit schematic is provided below.

A schematic of the one-port mBVD circuit topology showing the series resistance, , connected to the parallel combination of the motional and static branches. The motional branch includes , , and , which form a series contour, and the static branch includes and , which are connected in parallel.

The table below can be used to consider each element, how the elements impact the electrical admittance response, and how to cleanly extract their values from your FEM simulation.

Circuit Element	Physical Meaning	Admittance Impact	How to Extract It
Series resistance,	This element represents the purely ohmic resistance of the metal electrodes, busbars, and interconnect wires leading up to the active area of the resonator.	It limits the maximum peak of the admittance magnitude at series resonance .	can be calculated directly from the conductivity and geometry of the metal routing in your model, say via an auxiliary Electric Currents setup. You typically don’t need to include it explicitly in your main Piezoelectricity simulation.
Static capacitance,	represents the electrical plate capacitance, which is formed by the overlapping electrodes and the intervening piezoelectric material. is independent of mechanical movement.	It establishes the baseline “floor” of the imaginary admittance curve and determines the frequency separation between your series resonance and parallel antiresonance .	Run a frequency-domain simulation at an arbitrary point, , well below , where the device is mechanically “dead” and the motional branch acts as an open circuit. Evaluate the imaginary part of the admittance to isolate the static capacitance as
Static resistance,	represents the bulk dielectric material losses and electrical leakage paths through the piezoelectric layer.	It introduces a clean, flat baseline conductance to the real part of the admittance curve at low frequencies.	Assuming that your model goes with a nonzero dielectric loss, in the same low-frequency regime used to find , i.e., well below acoustic resonance, measure the real part of the admittance to directly isolate the flat baseline dielectric leakage as . For a case of small dielectric loss, is supposed to get a huge value.
Motional capacitance,	This element represents the mechanical elasticity or compliance of the resonator structure, which is intimately tied to the effective electromechanical coupling coefficient of the piezoelectric crystal.	It dictates the overall strength of the resonance. A higher coupling coefficient increases . Subsequently, it widens the window between and , which is a vital feature for wide-bandwidth RF filters.	can be pulled instantly using your extracted value and the resonance frequency ratio via .
Motional inductance,	This element represents the effective acoustic mass or inertia of the vibrating crystal lattice or thin-film membrane.	Together with , it determines the exact frequency location of the series resonance peak. This calculation follows the classical formula for resonance frequency .	Once you know and have identified from your simulation data, you can define the motional inductance as .
Motional resistance,	This resistance accounts for the total mechanical and acoustic energy dissipation in the device.	It dictates the height and sharpness of the admittance peak at resonance. A lower value signifies lower loss, leading to a much higher Q factor.	In the Frequency Domain sweep, the peak admittance at is dominated by the motional branch. The peak conductance simplifies to . Alternatively, you can pull the device’s mechanical Q factor, , from the Eigenfrequency study and use the following relation: .

Here we opted for the mBVD configuration featuring the static branch with and connected in parallel. The alternate convention with these elements connected in series also exists, and, if chosen, it will affect their extraction logic and interpretation. Two approaches are interchangeable in the sense that once the lumped parameters are obtained attentively, the outputs are to be identical.

Extracting the Parameters

Provided that we are now equipped with both theory and extraction logic, let’s extract the full set of lumped parameters and construct the mBVD circuit within COMSOL Multiphysics^®.

Note that this blog post acts as the primary resource for this model. If you’re interested in finding step-by-step instructions for leveraging similar models and plotting admittance curves, check out the example models we have available for the MEMS Module.

Here, we opt for a 3D layered one-port thin-film bulk acoustic resonator (FBAR), shown below. We also assume that is known beforehand and for further comparison, just add its contribution manually to the device’s admittance.

To accurately extract resonance frequencies and the peak conductance from the Frequency Domain study, you might want to refine the frequency step and use the Graph Markers for a Global admittance plot. Do not forget to assign the Dielectric Loss if you consider . Finally, it makes sense to run a separate study for a frequency several orders smaller than for the purpose of and extraction.

The UI in COMSOL Multiphysics^® showing a 3D FBAR model, with its admittance calculated via the Frequency Domain study.

Once all the parameters are calculated, you can bring these values directly into the built-in Electrical Circuit interface within the same model file to verify your work. Connect your newly defined lumped mBVD circuit to a voltage source and overlay its impedance response directly onto your original FEM data. Following the modern and handy AI trends, do not hesitate to ask the optional Chatbot window in COMSOL Multiphysics^® to provide you with COMSOL API code for adding the circuit with all the needed settings, as was done here.

The extracted mBVD parameters for the Electrical Circuit interface and the plot depicting good correlation between the full-fledged FEM model and the circuit representation.

The result is a great match around the operational bandwidth. From any practical standpoint, a minor discrepancy can be explained by the fact that the FEM model captures contributions from high-order modes present in the system. Once your mBVD circuit block matches the physics-based model, you are no longer constrained by the physical geometry. You can extend the circuit network right inside COMSOL Multiphysics^® by adding external tuning capacitors, inductors, or transmission lines, or by cascading multiple identical blocks to form complex ladder filters. Altogether, you end up with a verified, lightning-fast lumped model that seamlessly bridges microscale multiphysics with macroscopic system reality, enabling you to run massive sweeps in seconds.

Next Steps

Let’s wrap up by sharing the model featured in this blog post as well as related examples:

• Check out the model discussed in this blog post: Extracting mBVD Circuit Parameters from a 3D BAW Resonator
• To learn how to approach the same task from the Eigenfrequency side, see: Frequency-Domain Analysis of a SAW Unit Cell
• If you need a perfect match and/or a way to compensate for missing data (e.g., damping ratios), see this example, which makes use of the Parameter Estimation feature: Thin-Film BAW Resonator with Equivalent Circuit
• Check out tutorial models that leverage the MEMS Module: Application Gallery

The Agentic Loop

Until recently, most AI-assisted workflows focused on conversational interfaces: Users asked the AI tools modeling questions and received generated code or guidance while remaining in control of each individual step. The Chatbot window in COMSOL Multiphysics^® follows this paradigm, enabling engineers to ask modeling questions, generate and debug COMSOL API code, and streamline repetitive tasks, all within the software.

The next step in an AI-assisted workflow is moving from conversation to action. An AI agent can plan a sequence of tasks, execute them through software tools, inspect the results, and decide what to do next. In simple terms, an agent can be thought of as a wrapper or “harness” around a large language model (LLM). The harness applies the actions derived by the LLM. The LLM writes a piece of code and wraps it in a specific format to perform a tool call, which triggers the harness to, for example, open a command line and a command to run the piece of code.

Understanding Interactions Between Agents and COMSOL^®

COMSOL Multiphysics^® is a great fit for exploring the shift toward agentic workflows thanks to the COMSOL API, the software’s rich programming interface. Here, we will look at how agents can leverage the COMSOL API to assist engineers in their modeling workflows by acting as collaborators.

There are different ways to set up an agentic workflow in COMSOL Multiphysics^®. For interactive workflows, it is recommended to use a client–server setup to inspect an agent’s progress and the state of the COMSOL model, which lives on a COMSOL Multiphysics^® server (not to be mistaken with COMSOL Server™, an add-on product for deploying simulation apps). Both the COMSOL Multiphysics^® user and the agent can connect to the COMSOL Multiphysics^® server, see the model state, and make changes.

In a typical interaction between an agent and a live COMSOL model, the agent would:

1. Generate the source code for a Java program using the COMSOL API for use with Java
2. Compile and run the code using the Java commands bundled with the COMSOL Multiphysics^® installation. Through the Java program, the agent would then:
   1. Connect to the COMSOL Multiphysics^® server
   2. Query model properties, run studies, and/or make controlled changes
   3. Disconnect from the server
3. Reason over the output, plan, and perform next steps

How to Connect COMSOL^® and Codex

Now that we have covered the basics of how an agent interacts with a COMSOL^® model, we’ll go over the steps for setting up an agentic workflow in the software.

Since we have an agent available, we won’t have to complete the setup manually. Instead, we can ask the agent for help based on the architecture description. To do this, we start with a COMSOL Multiphysics^® server instance, connect the COMSOL^® user interface (UI), and load a model. (Also see: Running COMSOL Multiphysics^® in Client–Server Mode.) Then, we install and set up an agent. For this example, we will use Codex from OpenAI. Codex can create a skill that enables us to manipulate the model object located on the server. A skill bundles instructions and scripts in a package that you can reuse and share.

Here is the prompt we used to create a simple version of the skill:

Create a lightweight, local Codex skill named comsol-api. Install it as a repo-scoped skill. The skill lets Codex inspect and modify a COMSOL Multiphysics model located on a COMSOL Multiphysics server through the COMSOL API *for* use with Java.
 
The skill should work in a client-server configuration. The model is located on an already running COMSOL Multiphysics Server. The main part of the skill is a small wrapper script that compiles Java code using the COMSOL’s bundled Java compiler and API plugin classpath, and runs this program. The Java code should connect to the server, manipulate the model object and disconnect again. Before performing any changes, connect to the server and read the model tag.
 
Save the .java file of each agent interaction with the server in a workspace folder alongside a log of the output of the Java programs and the COMSOL progress log (showing computation time, memory requirement, degrees of freedom) for later reference.
 
Research here to learn more about the COMSOL API.
 
https://www.comsol.com/support/learning-center/article/overview-of-the-comsol-api-107912
 
Based on that, create an implementation plan and implement it. The skill should be intentionally light and only focus on running Java code against the already running server. No need for tests. Don’t install additional packages. Keep it package free.
 
Test the skill by adding a parameter and reading back the same parameter.

Luna — the most affordable current-gen OpenAI model, with a reasoning setting of medium — was able to implement a skill in under ten minutes.

The Codex agent implementing and validating a skill.

To test the skill, the agent added a global parameter and then read the same parameter back from the live model on the COMSOL Multiphysics^® server. During the test, the parameter appeared in the COMSOL Multiphysics^® UI.

Possibilities and Challenges of Agents

Now that we have established a connection between COMSOL^® and Codex, what can we do with it? Every click in the UI has a corresponding sequence of COMSOL API commands. Due to the comprehensive scope of the COMSOL API (which makes it possible to perform any part of the modeling workflow for any physics), Codex can perform any action.

The API can expose a broad variety of information for an agent to inspect, such as:

• Loaded model tags and component structure
• Parameters, expressions, descriptions, and units
• Feature types, labels, and tags
• Geometry and mesh sequences
• Physics interface and boundary condition settings
• Study and solver configuration
• Numerical results and evaluated expressions
• Warnings, errors, and solver logs

This information can be extracted, reasoned over, and manipulated, enabling agents to plan, run, and adjust experiments as necessary, as well as summarize the results in a report. This way, the agent can, for example, define parameters and settings, add physics interfaces, and create plots or numerical evaluations. Additionally, this agentic workflow enables you to perform sweeps that are usually not possible with the Parametric Sweep study. With the right setup, you can keep such a system running for an extended period of time.

Demonstration: Heat Sink Model

To demonstrate this capability, we will use the agent to experiment with nonlinear solver parameters in order to reduce the number of nonlinear iterations and speed up model computation, using the heat sink tutorial model as an example. The task is to analyze the available settings, select the three most promising ones, and run them with three values each. While the agent can traverse the model object and identify the corresponding node or search the COMSOL^® documentation, we extracted the relevant code using the Copy as Code to Clipboard submenu to help the agent. This functionality is especially useful here, as the connection we implemented is quite bare bones. The prompt we used is as follows:

I am looking to speed up the computation, specifically reduce the number of nonlinear iterations necessary to achieve convergence. Analyze the parameters of the segregated solver, identify which might be good candidates to tune. Pick the most promising three and run three values each. Analyze the computation time, number of iterations in a report. Make only changes that are meaningful. I extracted the values that can be set in the segregated solver node.
 
model.sol("sol1").feature("s1").feature("se1").set("segterm", "tol");
model.sol("sol1").feature("s1").feature("se1").set("maxsegiter", 200);
model.sol("sol1").feature("s1").feature("se1").set("ntolfact", 1);
model.sol("sol1").feature("s1").feature("se1").set("segtermonres", "auto");
model.sol("sol1").feature("s1").feature("se1").set("segreserrfact", 1000);
model.sol("sol1").feature("s1").feature("se1").set("segstabacc", "segcflcmp");
model.sol("sol1").feature("s1").feature("se1").set("segcfltech", "interp");
model.sol("sol1").feature("s1").feature("se1").set("subinitcfl", 5);
model.sol("sol1").feature("s1").feature("se1").set("submincfl", 10000);
model.sol("sol1").feature("s1").feature("se1").set("subforcecfl", true);
model.sol("sol1").feature("s1").feature("se1").set("subkppid", 0.65);
model.sol("sol1").feature("s1").feature("se1").set("subkipid", 0.15);
model.sol("sol1").feature("s1").feature("se1").set("subkdpid", 0.15);
model.sol("sol1").feature("s1").feature("se1").set("subcfltol", 0.2);
model.sol("sol1").feature("s1").feature("se1").set("subadaptcfltol", true);
model.sol("sol1").feature("s1").feature("se1").set("segcflaa", "on");
model.sol("sol1").feature("s1").feature("se1").set("segcflaacfl", 9000);
model.sol("sol1").feature("s1").feature("se1").set("segcflaaset", "pid");
model.sol("sol1").feature("s1").feature("se1").set("segcflaadim", 10);
model.sol("sol1").feature("s1").feature("se1").set("segcflaamix", 1);
model.sol("sol1").feature("s1").feature("se1").set("segcflaadelay", 0);
model.sol("sol1").feature("s1").feature("se1").set("segcflaafact", 1);
model.sol("sol1").feature("s1").feature("se1").set("segcfljtech", "off");
model.sol("sol1").feature("s1").feature("se1").set("consistencycheck", "strictinitially");
model.sol("sol1").feature("s1").feature("se1").set("plot", "off");
model.sol("sol1").feature("s1").feature("se1").set("probesel", "none");

Your prompt can, of course, include additional resources, including Learning Center articles, blog posts, or Knowledge Base articles from the COMSOL^® website.

The results of the agent performing solver tuning studies.

Running nine experiments, the agent identified a relevant setting and reduced the number of nonlinear iterations by 25%, leading to a 30% faster computation time!

Another test we can run is to try out the difference between the three different compressibility options in the Laminar Flow interface. Although we cannot use the Parametric Sweep study in the UI to switch between the options in the dropdown menu directly, we can ask Codex to test each option, print the temperature in the chip, and list the number of nonlinear iterations. To instruct Codex to perform these tasks, here is the prompt we used:

There are three different compressibility options in the Laminar Flow interface. Run a test and compare all of them. I added an Evaluation Group that evaluated the temperature in the chip that you can run with the command:
 
model.result().evaluationGroup("eg2").run();
 
Create a report with the average temperature, computation time and number of iterations for the three different compressibility options.

There is no substantial change between the options, and the computation time is almost identical. So, we will stick with the default option. Please note that this was only an example on how to perform sweeps, and any settings in the interfaces should be taken based on physical considerations, not purely on computational speed.

Considerations When Using Agents

Although our experiments were successful, there are still challenges associated with agentic workflows. One of the challenges is spatial perception and spatial reasoning. For a human user, it is easy to select a boundary for a boundary condition by looking at the geometry. Under the hood, the number of the boundary is assigned to the selection of the respective boundary condition. Identifying this number purely by using the API can be difficult. Possible ways to circumvent this challenge is the rigorous use of named selections. Alternatively, you can activate the entity number in the view and instruct the agent to automatically export images from different view angles. This method is helpful as long as the geometry is not too complex.

It is also important to check the output to see if the changes are really meaningful. Considering the example of nonlinear solver tuning, the agent could simply change the termination technique to Iterations instead of Tolerance and set the maximum number of iterations to one. This method would reduce the number of iterations but would not give a numerically sound result. Monitoring key quantities, like the temperature of the chip, and ensuring that they do not change substantially is one possible way to filter out such bad candidates.

Since agents like Codex can do more than simply output text, there are risks associated with running them. They can make incorrect assumptions, delete files, expose sensitive data, or read malicious prompts on the websites they ingest. Because agents are inherently stochastic systems, they can make mistakes. To use these tools more safely, the agents should be operated in a sandbox, only have access to the files they really need, and require approval for high-risk actions.

Final Thoughts on Connecting Codex to COMSOL^®

What we showed here is a minimal working example of a connection of Codex to COMSOL Multiphysics^® for agentic workflows. This skill and wrapper script can be further extended to enable a more robust interface handling, which would be especially useful for identifying overly long simulation or mesh creation times for robust recovery.

These advancements can change how engineers interact with simulation software. The general scope of the COMSOL API makes COMSOL Multiphysics^® well suited for updating simulation workflows to meet the era of AI.

OpenAI is a trademark of OpenAI, Inc. Oracle and Java are registered trademarks of Oracle and/or its affiliates.