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 Multiphysics 6.4 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.

Read the previous blog post here.

Why Accessibility Matters

Color vision deficiency (CVD) can make it difficult to distinguish between certain colors. The most common form is red–green color blindness, which affects the ability to differentiate between shades of red and green. People with CVD may find it difficult to analyze or may misinterpret results presented in a standard rainbow color table or grayscale. Because simulation results are often shared with a wide audience, choosing an accessible color table can help make the results easier for more people to interpret as intended.

The electric potential of an atmospheric busbar shown using the Cividis color table option available in COMSOL Multiphysics^®.

Built-In Accessible Color Tables

COMSOL Multiphysics^® includes several built-in color tables that are both perceptually uniform and well suited for accessible scientific visualization including Cividis, Viridis, Inferno, Magma, and Plasma.

The tables being “perceptually uniform” means equal changes in your data appear as equal changes in the colors. Accessible color tables provide smooth, consistent changes in both color and brightness, making it easier to identify gradients. These color tables provide a range of options for presenting simulation results clearly while remaining accessible to a broad audience. Whether you’re visualizing temperature distributions, stress concentrations, fluid flow, or electric fields, your choice of color table plays an important role in communicating the model’s underlying physics.

Example mesh plots, propagating waves, wake regions, singularities, and smooth gradients for each color table.

While all five color tables can be used across a wide range of applications, some may be better suited for different types of simulation results. Viridis and Cividis are good choices for stress, strain, pressure, velocity, and concentration, where smooth gradients help reveal variations. Inferno and Magma can help emphasize temperature distributions and thermal gradients, providing strong contrast. Plasma can be useful for quantities such as electric and magnetic field magnitude and current density, where its distinctive color progress can help highlight variations in the field.

In addition to these perceptually uniform color tables, COMSOL Multiphysics^® includes the GrayScale linear color table. This option can be particularly useful when hue should not influence the interpretation of a result, or when you want to emphasize the intensity of a quantity without introducing additional color information. However, note that humans can only distinguish between about 30 different shades of gray, so it may be more challenging to notice small changes in a grayscale format.

A plot showing the interpolation function (left) of a model resembling the Alpine mountain the Matterhorn. The grayscale image (right) shows the height of the mountain. Note that the color bar is normalized to go from 0 to 1.

Transforming Color Tables

Color table transformations available in the settings of results plots, including Reverse, Nonlinear, and Nonlinear symmetric options, provide additional flexibility in how colors can be used to represent simulation results. They can be used to adjust the appearance of a color table and change how variations in the data are displayed. A transformation can be used to adjust the appearance of a color table and change how variations in the data are displayed. It can be used to emphasize or de-emphasize different parts of a result and provides another way to tailor a color table to a particular application or type of data.

The available color table transformations.

The von Mises stress in a bracket without a color table transformation applied (left). The von Mises stress in a bracket with the Nonlinear color table transformation applied. (right)

Custom Color Tables

In addition to the built-in color tables, version 6.1 of COMSOL Multiphysics^® introduced support for custom color tables. You can import color table files, edit the RGB components of a color table directly, and refer to color tables stored on disk or in a database.

This flexibility makes it possible to develop custom visualization styles for specific applications while maintaining control over how simulation results are represented.

To add a Color Tables node to the model tree, open the Show More Options dialog and select the Color Tables checkbox.

The Color Tables node used to define a custom table. When values are defined for a custom color table, a preview will appear in the Settings window.

Accessible Visualizations

Choosing an appropriate color table is an important part of communicating simulation results. By using perceptually uniform and accessible color tables, you can create visualizations that accurately represent your data while making your results easier for a wider audience to interpret. With the built-in Cividis, Viridis, Inferno, Magma, and Plasma color tables, as well as GrayScale and support for custom color tables, COMSOL Multiphysics^® provides a range of options for creating clear and accessible simulation visualizations.

Read Part 1 of this two-part blog series here.

The Importance of Acoustics in a Meeting Room

Have you ever been in a meeting where you can barely hear what the other person is saying? This is often due to excessive reverberation in the room. It is well known that a bad sound environment in the workplace can lead to poor communication, focus, and productivity. The best solution to this problem is to install acoustic absorbers, such as a suspended porous ceiling or an absorbing curtain, in a suitable place in the room. These absorbers have an air layer behind them to absorb low-frequency sound, making them highly extended reactive. It is interesting to see the difference that using different sound-absorbing boundary models makes when simulating a room with such extendedly reacting sound-absorbing functions.

Room Model and Absorption Setting

The simple room (3.5 m 6.8 2.5 m) below was analyzed with octave bands ranging from 63 to 500 Hz (approximately 40 to 710 Hz) using the Pressure Acoustics, Time Explicit interface.

A diagram showing the geometry of an analyzed meeting room.

The Gaussian pulse was applied to the source point (0.7 m, 4.15 m, 1.2 m), and at the receiver point (2.8 m, 2.35 m, 1.2 m), the pressure response was evaluated up to 0.8 s with a sampling frequency of 10000 Hz.

To solve the acoustic wave equation in the time domain, the Pressure Acoustics, Time Explicit interface uses the discontinuous Galerkin finite element method (dG-FEM), which is a memory-efficient and scalable method. The room model has typical sound absorption functions, i.e., the suspended porous ceiling and the curtain. In this blog post, these absorbers are modeled using both locally and extendedly reacting models, and the resulting differences are examined. The image above is for the extended reaction model, which includes the thickness of the porous material and the air space behind the absorbing materials. In the case of the local reaction model, the thickness and air space are not modeled explicitly, while the surface impedance at normal incident condition is assigned to the surface boundaries of the materials. The surface impedances of those sound absorbers were evaluated with submodels in the frequency domain, assuming a two-dimensional periodic structure.

The suspended ceiling is made of a 50 mm porous material. The air space behind the porous material is 180 mm and was made to absorb low-frequency sound. Since the thickness of the porous layer is large, an equivalent fluid model is used to model the ceiling absorber, including the extended reaction. The Johnson–Champoux–Allard (JCA) model is used to model the fluid property in the porous material. With this model, the complex effective density and bulk modulus of the porous material are expressed as follows:

Therein, , , and are respectively the absolute pressure, the density, the ratio of specific heats, and the dynamic viscosity of the fluid in the porous material. , , , and represent the porosity, the flow resistivity, the viscous characteristic length, the thermal characteristic length, and the tortuosity factor of the porous material. is the angular frequency. The table below summarizes the parameters of the porous material used in the example model.

Parameter (Unit)	(1)	(Pa·s/m2)	(mm)	(mm)	(1)
Value	0.99	17000	0.14	0.15	1.01

Porous matrix properties used in the example model. These values are from Ref. 1.

The Poroacoustics feature, available in the Pressure Acoustics, Time Explicit interface, can be used to model the porous layer based on the equivalent fluid model. The Porous Absorber with Local and Extended Reacting Approximations for Time-Domain Modeling document is a good reference for our implementation of the time-domain equivalent fluid model. To use this feature, we should prepare the frequency-dependent data of the complex effective density and compressibility (a reciprocal of the bulk modulus) of the porous material. In the example model, these properties were calculated using the same submodel used to evaluate the surface impedances for local reaction modeling. The following screenshot shows the Poroacoustics settings for the submodel.

Poroacoustics feature settings for submodeling of the porous material.

The following screenshots show more calculations of the material properties from the result of the submodeling.

Calculation of the complex effective density (left) and compressibility (right) of the porous material from the result of the 2D frequency-domain analysis.

You could also use the measured data of these material properties. Then, the frequency-dependent properties of the porous materials are imported and approximated to the rational function forms with the Partial Fraction Fit function for the analytic inverse Fourier transformation.

Use of the Partial Fraction Fit function to approximate the complex effective density (left) and compressibility (right) to the rational function forms.

The fitted results are further imported in the Poroacoustics feature in the Pressure Acoustics, Time Explicit interface.

Poroacoustics feature in the Pressure Acoustics, Time Explicit interface.

On the other hand, the extended reacting absorbing curtain was modeled using a transfer impedance, on the assumption that it was very thin. We assumed that the flow resistance (transfer impedance) of the curtain was 416 Pa·s/m. Note that the mass effect was neglected here for simplicity. The curtain was placed on a window frame with an air space of 200 mm. The Interior Impedance feature was used for modeling.

Note that to compare the different boundary models, the boundary surfaces of the spaces behind the extended reactive materials were assumed to be rigid. Other boundaries were assumed to be plaster board and flooring. They were modeled using the frequency-dependent impedance boundary. The random incidence absorption coefficients for all materials are shown below:

Random incidence absorption coefficients for all materials used in the room model.

The ceiling and curtain show the different absorption coefficients for the different boundary models due to the high extended reactivity of these absorbers. The local reaction models show larger absorption coefficient values than the extended models over a wide frequency range, except above 630 Hz. You can see the angle-dependent absorption coefficients for the absorbing materials from the model file.

Normalizing and Filtering of the Impulse Response

The Gaussian pulse is a typical source signal in wave-based, time-domain modeling, but its frequency spectrum is not flat, as will be shown later. Thus, to use dG-FEM results for room acoustics evaluation, the source spectrum needs to be normalized. Here, we introduce how to perform the normalization in the frequency domain.

First, we introduce the frequency characteristic of the Gaussian pulse. The Gaussian pulse was excited to sound fields with the following initial sound pressure distribution:

Here, represents the sound pressure. is the distance from the source, and is the parameter characterizing the source spectrum. Substituting this initial sound distribution function to the wave equation for a spherical wave yields the following relation in a free field (Ref. 2).

Therein, is the time variable, is the coordinate vector, and is the sound speed. With the Green’s function for the three-dimensional free field, the volume acceleration of the point source is expressed as follows (Ref. 3):

Here, is the density of the medium. By Fourier transformation, the spectrum of the source signal using the Gaussian pulse is expressed as follows:

Here, is the frequency variable.

The spectrum of the source signal used in the example model.

The source was designed so that the gain of -3 dB from the peak was obtained at the upper-limit frequency in the 500 Hz octave band.

Thus, we can obtain the normalized impulse response with a source strength of 1 m³/s² volume acceleration (around 0.14 mW source power in the air) by dividing the results by the analytic spectrum above. However, this normalization also amplifies the unwanted frequency components, such as a higher frequency component that is not resolved in the dG calculation. To remove such components, we also perform filtering. The following equation expresses the normalization procedure:

Here, and represent the normalized and dG-FEM-calculated sound pressure. is the Fourier transformation operator. and are the high-pass and the low-pass filters in the frequency domain, respectively.

Since the filtering process may result in a noncausal signal, zero padding is recommended before converting the dG results to the frequency domain via discrete Fourier transformation. The example model exported the pressure waveform calculated by the dG-FEM to a WAV file and imported it as an interpolation function. Then, zero padding was performed over the Grid 1D dataset as follows:

Definition of an interpolation function using a WAV file exported from the calculation result. To perform the zero padding, the extrapolated value is set to 0.

Zero-padded Grid 1D for discrete Fourier transformation. The negative part of the interval is the zero-padded region. The resolution is set so that the sampling frequency corresponds to 10,000 Hz.

In the example model, the normalization and filtering were done with the following Global ODEs and DAEs interface.

Here, pE_normalized represents the normalized sound pressure, and rpE(freq) and ipE(freq) express, respectively, the real and imaginary part of the Fourier coefficient of the dG-FEM-calculated sound pressure for the extended reaction model. pL_normalized, rpL(freq), and ipL(freq) represent those same values, but for the local reaction model. source(freq) denotes the spectrum of the source signal defined with an analytic function. HP(freq) and LP(freq) are the high-pass and low-pass filters defined using step functions. You can see more detail of the setting from the example model. The frequency domain data of pE_normalized and pL_normalized are calculated using the Frequency Domain study step. Then, the normalized impulse responses are reconstructed via the Frequency to Time FFT study step.

Another benefit of the normalization is that the results can be combined with the ray acoustics model results to obtain the broadband impulse response using a hybrid approach. This is done by setting the ray release condition to be omnidirectional with the following total source power :

Here, is the unit volume acceleration of 1 m³/s².

Results

The normalized band-limited room impulse responses and their spectra (room transfer function) for two sound-absorbing boundary models are shown below.

Normalized room impulse responses for the extended and the local reaction models.

Normalized room transfer function for the extended and the local reaction models.

The differences between the two models using different sound absorbing models are clear. The local reaction model overestimates the absorption by the sound-absorbing materials over a wide frequency range, which results in a faster decay of the impulse response. According to the comparison of the spectra, the local reaction model has more absorption below 160 Hz and around 500 Hz, while it underestimates the absorption above 600 Hz. These tendencies correspond to the differences in the random incidence absorption coefficients between two sound-absorbing boundary models.

The following comparisons of the room acoustics metrics also indicate the differences between the two models.

A comparison of the reverberation parameters (EDT and T20) between the extended and local reaction models.

A comparison of the clarity (C50) between the extended and local reaction models.

As for the reverberation parameters (EDT and T20), the differences between the two models are larger than the two times of just noticeable difference (JND) of 5% in a broadband. The clarity parameter (C50) also showed a larger difference than the JND (1.1 dB) at most frequencies, with the maximum difference reaching 11 dB. These results indicate that the difference between the local and the extended models is audible. You can hear the normalized and unnormalized room impulse responses for the different absorbing boundary models in the following:

Playback audio waveforms of the normalized impulse responses for the extended (left) and the local (right) reaction models.

Playback audio waveforms of the unnormalized impulse responses for the extended (left) and the local (right) reaction models.

You can feel the longer reverberation for the extended model. Also, the unnormalized results will be heard as emphasizing high-frequency sounds due to the frequency characteristic of the Gaussian pulse.

These results are at one receiver location, and it is important to set up more receiver locations for the practical evaluation of room acoustics. However, the clear differences at the possible receiver location are sufficient to demonstrate the effect of the type of the sound-absorbing boundaries on the meeting room.

The wave-based methods are considered to be inherently very accurate because they include all wave phenomena. However, the results of the example model showed that depending on the characteristics of the sound absorber, a simple use of the local reaction model may lead to erroneous results. In order to perform a more reliable simulation, we should carefully confirm the configuration and the absorption characteristics of the absorbers.

Conclusion

As a continuation of the previous blog post, this blog post examines the effects using different sound-absorbing boundary types on a small meeting room. Typical sound absorbing elements, like the suspended ceiling and the curtain, are highly extended reactive. As a result, applying the conventional local reaction model (the impedance boundary) to these absorbers produced a less accurate result. The demonstration revealed that the capability of the extended sound-absorbing boundary model is very crucial to perform a proper simulation. Of course, the impedance boundary is still convenient because of its lower computational cost than the extended model. For efficient modeling, the use of different absorbing models, depending on the design stage or the characteristics of absorbers, is very important. This blog post also introduced the procedure for obtaining a band-limited impulse response from the dG-FEM analysis with the Gaussian pulse excitation. Note that when exciting a sound field with the initial pressure distribution condition, the source location must be apart from any boundary surfaces reflecting sounds.

The Pressure Acoustics, Time Explicit interface in the Acoustics Module, an add-on to the COMSOL Multiphysics^® software, is well suited for wave-based room acoustics modeling because of its high-memory efficiency and accuracy, thanks to the ability to model both locally and extendedly reacting boundary models, including the frequency dependency. In addition, you can significantly reduce the computation time by using a GPU-accelerated formulation (currently available for only local reaction model). Even with an entry-level GPU (NVIDIA^® T400), the local reaction model (1,307,650 DOFs and 40,000 timesteps) can be solved in about an hour. If you have a higher-grade graphics card, the reduction of computation time is expected to be much greater. You can find more details about the performance gain by GPU in the following document: Acoustics of an Open-Plan Office Space.

Let’s model room acoustics with the exceptionally accurate and efficient simulation capabilities available in COMSOL Multiphysics^®!

Next Steps

Interested in trying out the model for yourself? Download the related MPH file below:

Reference

1. H. Wang and M. Hornikx, “Extended reacting boundary modeling of porous materials with thin coverings for time-domain room acoustic simulations,” J. Sound Vib., vol. 548, 117550, 2023; https://doi.org/10.1016/j.jsv.2022.117550.
2. S. Sakamoto, “Phase-error analysis of high-order finite difference time domain scheme and its influence on calculation results of impulse response in closed sound field,” Acoust. Sci. Technol., vol. 28, 295-309, 2007; https://doi.org/10.1250/ast.28.295.
3. T. Okuzono, T. Otsuru, R. Tomiku, and N. Okamoto, “Application of modified integration rule to time-domain finite-element acoustic simulation of rooms,” J. Acoust. Soc. Am., 1vol. 32, 804–813, 2012; https://doi.org/10.1121/1.4730920.