COMSOL

COMSOL Multiphysics solves coupled PDEs with the finite element method in real geometry. Learn the modelling workflow (geometry, materials, physics, mesh, study, results), coupling flow + heat + species + reaction, boundary conditions, mesh convergence, and verifying against analytical solutions — with a scripted example.

Some problems cannot be solved with lumped balances: temperature hot spots in a catalyst pellet, flow maldistribution in a microreactor, concentration gradients in a membrane. COMSOL resolves the spatial detail by solving the governing PDEs on a mesh.

Its strength is multiphysics — coupling fluid flow, heat transfer, chemical species transport and reactions in one model.

An analogy

COMSOL is like a microscope for equipment: it zooms into the inside of a reactor or exchanger and shows how temperature and concentration vary point by point.

The modelling workflow

Geometry (draw or import CAD) → Materials → Physics interfaces (Laminar Flow, Heat Transfer in Solids/Fluids, Transport of Diluted Species, Chemistry) → Boundary conditions → Mesh → Study (stationary, time-dependent, parametric) → Results (plots, probes, derived values).

Boundary conditions

Dirichlet (fixed value, e.g. wall temperature, inlet concentration), Neumann (fixed flux, e.g. insulation or symmetry = zero flux), Robin (convection, h(T − T∞)), inlet/outlet velocity or pressure, no-slip walls. Wrong BCs are the most common source of wrong results.

Mesh and convergence

Finer meshes near walls and reaction fronts (boundary-layer meshes). Refine until key outputs change by < 1–2 %. Check solver logs for convergence and conservation (mass/energy balances integrated over boundaries).

Chemical engineering multiphysics

Packed-bed reactors (porous media flow + reaction + heat), microreactors, electrochemical cells (batteries, electrolysers), membrane separations, CVD reactors, and mixing tanks with rotating machinery. Use nondimensional numbers (Re, Pe, Da, Bi) to check regimes; couple to Aspen or MATLAB for system-level studies.

Formulas

Worked example: Verify a heat-generating slab

Given: q̇ = 10⁶ W/m³, k = 0.5 W/m·K, half-thickness 10 mm, T_s = 300 K

Find: Maximum temperature.

Answer: 400 K — the COMSOL result must match before adding complexity

Worked example: Is convection important? (Péclet)

Given: Microchannel: u = 1 mm/s, L = 1 mm, D = 10⁻⁹ m²/s

Find: Pe.

Answer: Pe ≫ 1 → convection dominates; need fine mesh or stabilisation

COMSOL model via LiveLink for MATLAB

import com.comsol.model.*
import com.comsol.model.util.*

model = ModelUtil.create('Model');
comp  = model.component.create('comp1', true);
geom  = comp.geom.create('geom1', 1);          % 1D slab
geom.create('i1', 'Interval');
geom.feature('i1').set('coord', [0 0.01]);     % 0 .. L = 10 mm (half-thickness)
geom.run;

ht = comp.physics.create('ht', 'HeatTransfer', 'geom1');
ht.feature('solid1').set('k', {'0.5'});        % W/m/K
hs = ht.create('hs1', 'HeatSource', 1);
hs.selection.all; hs.set('Q0', '1e6');         % W/m^3
ht.create('sym1', 'Symmetry', 0);  ht.feature('sym1').selection.set(1);
ht.create('temp1', 'TemperatureBoundary', 0);
ht.feature('temp1').selection.set(2);
ht.feature('temp1').set('T0', '300[K]');

comp.mesh.create('mesh1').autoMeshSize(3);     % refine and compare
model.study.create('std1').create('stat', 'Stationary');
model.study('std1').run;

Tmax = mphmax(model, 'T', 'line');             % expect 300 + 100 = 400 K

The same steps as the GUI — geometry, physics, boundary conditions, mesh, study — scripted for automation and parametric sweeps. The analytical answer (ΔT = q̇L²/2k = 100 K) verifies the model before adding complexity.

What the solver does — finite differences in Python

import numpy as np
# -k T'' = q on 0<x<L, T'(0)=0 (symmetry), T(L)=Ts
k, q, L, Ts = 0.5, 1e6, 0.01, 300.0
for n in (5, 10, 20, 40):                 # mesh refinement study
    x = np.linspace(0, L, n + 1); h = x[1] - x[0]
    A = np.zeros((n + 1, n + 1)); b = np.full(n + 1, q * h**2 / k)
    for i in range(1, n):
        A[i, i-1:i+2] = [-1, 2, -1]
    A[0, 0:2] = [2, -2]                   # ghost-node symmetry
    A[n, n] = 1; b[n] = Ts
    T = np.linalg.solve(A, b)
    print(n, round(T[0], 3))              # -> 400.0 at every mesh (exact for quadratic)

FEM and FD both turn the PDE into a linear system on a mesh. For real problems the answer changes with mesh size — a mesh-convergence study (refine until the result stops changing) is mandatory before trusting results.

Common mistakes

Applications

Key takeaways

Going further