Crystal Growth Simulation

The problem

Czochralski growth pulls a single silicon crystal from a rotating melt. The flow inside that melt — buoyancy-driven convection coupled to rotation-driven swirl and heat transfer — determines the quality of the resulting crystal. Resolving it with computational fluid dynamics is accurate but slow, which makes exploring the process parameter space expensive.

A surrogate model that reproduces the converged flow field in milliseconds, across operating conditions it was never trained on, turns that search from an overnight batch into an interactive one.

Reference data

We validate against certified CFD reference data for a 2D axisymmetric silicon melt (radius 0.30 m, height 0.15 m), solved in ANSYS Fluent as steady-state cases on a consistent 8,181-node mesh. The simulation geometry and boundary conditions follow Huang et al. (2025), AIP Advances 15(10), 105202.

The benchmark spans three independent parameter sweeps around a shared anchor condition:

• Crystal rotation — 4 to 20 rpm
• Hot wall temperature — 1730 to 1785 K
• Crucible rotation — −1 to −10 rpm

Every case passed a four-pillar certification before inclusion: residual convergence, mesh independence, mass and energy conservation, and programmatic checks of the physical boundary conditions. CFD data generation and certification by Bertwin Kurisinkal Shine, TU Chemnitz.

How we measure

The test set is held out by operating point, not by mesh node. Because every case shares the same mesh, splitting individual nodes at random would place neighbouring points of the same solution on both sides of the split — that measures interpolation within a case already seen, not generalisation to a new one. Whole cases are therefore reserved, and the model is scored on operating points it was not trained on.

We are explicit about what that does and does not demonstrate. The reference data consists of three one-dimensional sweeps — crystal rotation, temperature and crucible rotation, each varied on its own from a shared anchor condition. Holding out whole cases therefore measures interpolation along a parameter that was swept, including at the ends of each range, where accuracy holds up. It does not establish that the model generalises to a parameter direction absent from training: when we withhold an entire sweep and ask the model to predict its effect, accuracy collapses. Extending the benchmark with cases that vary several parameters together is the next step, and until those exist we report the current figures as interpolation across the swept operating range rather than as general predictive capability.

Accuracy is reported as the coefficient of determination in physical units, separately for each of the five predicted fields: radial, axial and swirl velocity, pressure, and temperature. A single error figure averaged over normalised fields can fall while physical accuracy worsens, because it rewards ignoring whichever field happens to be smallest. Per-field reporting makes that failure mode visible: a score below zero means a prediction is worse than simply guessing that field's mean.

Results at unseen operating points

Trained on 21 operating points and evaluated on held-out conditions, the surrogate reproduces the converged fields with the following coefficients of determination:

Every field is above 0.997, averaged over five independent training runs on 21 operating points and scored on six held out. Agreement is tight across runs: the weakest single result across all seeds and fields is 0.9968.

The meridional components — radial and axial velocity — used to be the weak point at around 0.8. They are one to two orders of magnitude smaller than the swirl in this regime, which makes them sensitive both to how the fields are scaled before training and to which spatial frequencies the network can represent. Fixing those two things closed the gap without enlarging the model: the network above has roughly 127,000 parameters and trains in about seven minutes on a single CPU node.

Where the physics term earns its place

Our models combine a data loss against reference solutions with a residual penalty from the governing equations. On a densely sampled mesh the data alone already determines the field, and the physics residual adds little. Its value appears when reference data is scarce: below roughly one percent of the mesh, including the physics residual produces a statistically significant improvement over the identical model trained on data alone.

That is the regime that matters commercially. High-fidelity CFD is expensive precisely because each case costs compute, so a surrogate that needs fewer reference cases to reach a usable accuracy is worth more than one that needs many.

The hybrid architecture

The network combines classical residual blocks with a parameterised quantum submodule whose expectation values are concatenated into the prediction head. Physical boundary conditions are enforced structurally rather than learned — radial and swirl velocities are constructed to vanish on the axis of symmetry.

We hold ourselves to a like-for-like comparison: every hybrid configuration is run against a width-matched classical control, so any difference is attributable to the quantum submodule and not to parameter count. On this benchmark the clearest measured effect of the quantum submodule is improved robustness — under deliberately unfavourable loss scaling, where the classical control degrades sharply, the hybrid retains substantially more accuracy. We report that as we find it, and continue to investigate the regimes where hybrid models offer an advantage.

Benchmark figures are from internal evaluation on the reference dataset described above, cover interpolation across the swept operating range, and are updated as the work progresses. For methodology details or collaboration enquiries, please get in touch.