Simulating Smart Temperature Control in Electric Kettles with COMSOL: A Multiphysics Modeling Guide

Publish on: 2025/12/10 Classify at: RESEARCH/Simulation

Words: 3162 Read:≈ 15min

Summary

Develop a numerical model for the temperature maintenance process in electric kettles using COMSOL Multiphysics. This article explores how free convection, thermostat control, heat transfer, and species transport interact in a common household appliance—from theory to MATLAB implementation.

Introduction

Picture this: you’ve just boiled water for tea, but get distracted by a phone call. Ten minutes later, your electric kettle has kept the water perfectly warm—not too hot, not too cold. This seemingly simple “keep warm” function actually involves a sophisticated dance of physics.

Natural convection, often overlooked in everyday life, plays a starring role here. The same phenomenon that drives ocean currents and weather patterns also governs the gentle circulation of water in your kettle. From chemical reactors to fermentation vessels, understanding these convective flows is key to optimizing thermal systems.

What makes this problem particularly interesting is the coupling of multiple physical processes:

Heat Conduction: Energy flowing through the glass walls and steel heating plate
Natural Convection: Gravity-driven water circulation caused by temperature differences
Heat Loss: Energy escaping from the kettle surface to the surrounding air
Species Transport: How dissolved substances get mixed by the swirling water
Thermostat Control: The on-off switching logic that maintains your target temperature

In this article, we build a 2D axisymmetric COMSOL model to simulate these phenomena, drawing inspiration from the classic “Free Convection in a Water Glass” example in the COMSOL Application Gallery.

Physical Problem Description

Geometric Model

To keep the simulation tractable, we simplify the electric kettle as an axisymmetric frustum-shaped container. Thanks to rotational symmetry, we can capture all the essential physics in a 2D model—dramatically reducing computational cost without sacrificing accuracy.

Geometric Parameters:

Parameter	Value	Description
R_top	7 cm	Top radius of kettle
R_bottom	7.5 cm	Bottom radius of kettle
H_kettle	16 cm	Kettle height
T_wall	1.3 mm	Glass wall thickness
T_bottom	3 mm	Steel bottom thickness

Physical Processes

Heat Conduction

In the solid regions—the glass walls and steel heating plate—heat travels by conduction, following Fourier’s law:

$q = -k \nabla T$

Here, $k$ is the thermal conductivity and $T$ is temperature. For the glass wall, we use the thermal properties of silica glass.

Natural Convection (Boussinesq Approximation)

The magic happens in the water. When the bottom heats up, warmer water becomes less dense and rises, while cooler water near the walls sinks down to take its place. This creates a continuous circulation pattern—classic Rayleigh-Bénard convection.

We model this density variation using the Boussinesq approximation:

$\rho = \rho_{ref} [1 - \alpha_p (T - T_{ref})]$

where $\alpha_p$ is the thermal expansion coefficient. The resulting buoyancy-driven flow forms characteristic recirculation cells that efficiently mix and heat the water.

Species Transport

Beyond heat, we also track how dissolved species (think minerals or tea extract) get distributed. The concentration evolves through:

Convective transport: Species carried along by the flowing water
Diffusion: Molecules spreading from high to low concentration

This coupling beautifully demonstrates how natural convection accelerates mixing far beyond what diffusion alone could achieve.

Thermostat Control Logic

The thermostat employs a hysteresis control strategy to prevent rapid on-off cycling:

1
2
When T_avg > 50.5°C → Turn off heater (HeaterState = 0)
When T_avg < 49.5°C → Turn on heater (HeaterState = 1)

The heating plate temperature switches according to state:

Heating state: T_heater = 80°C
Off state: T_heater = 25°C (room temperature)

COMSOL Model Setup

Physics Interfaces

Translating our physical understanding into COMSOL requires four key physics interfaces:

Interface	Identifier	Purpose
Laminar Flow	spf	Simulate buoyancy-driven water flow
Heat Transfer in Fluids	ht	Simulate heat conduction and convection
Transport of Diluted Species	tds	Track species concentration distribution
Events	ev	Implement thermostat switching logic

Multiphysics Coupling

Non-Isothermal Flow: This is where the magic happens. By coupling Laminar Flow with Heat Transfer and enabling the Boussinesq approximation, temperature differences directly drive fluid motion.

Species Transport Coupling: The diluted species interface “rides along” with the velocity field, letting us track how substances get stirred by the convective currents.

Material Properties

Water (Temperature-Dependent Properties)

Water’s physical properties vary significantly with temperature—a critical factor for accurate natural convection modeling:

Temperature	Density ρ (kg/m³)	Viscosity η (mPa·s)	Thermal Conductivity k (W/m·K)	Heat Capacity Cp (J/kg·K)
20°C	998.2	1.002	0.598	4182
50°C	988.1	0.547	0.644	4181
80°C	971.8	0.355	0.670	4198

Note: The model uses piecewise polynomial fits of these properties. See the Appendix for the complete MATLAB implementation.

Solid Materials

Material	Thermal Conductivity	Density	Heat Capacity
Silicate Glass	1.38 W/(m·K)	2203 kg/m³	703 J/(kg·K)
Structural Steel	44.5 W/(m·K)	7850 kg/m³	475 J/(kg·K)

Boundary Conditions

The boundary conditions are carefully chosen to represent realistic physical scenarios:

Boundary	Condition Type	Setting
Heating plate bottom	Temperature boundary	T = T_heater_off + (T_heater_on - T_heater_off) × HeaterState
Kettle sidewall exterior	Natural convection	External natural convection with heat transfer coefficient, T_amb = 25°C
Water surface	Convective heat flux	h = 2 W/(m²·K), driven by temperature difference: $q = h(T_{ext} - T)$
Water surface	Slip wall	No tangential stress (free surface approximation)
Glass-water interface	No-slip wall	Zero velocity at solid boundaries
Axis of symmetry	Axial symmetry	Symmetry condition for all fields

Note: Because of the relatively low values of heat transfer film coefficients at the top and side surfaces, a significant portion of heat may be conducted through the bottom boundary when the heater is off. The Heat Transfer Module provides a library of heat transfer coefficient functions that can be utilized for more accurate modeling.

Event Control (Thermostat)

The COMSOL Events interface implements the temperature control logic:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
% Discrete state variable
HeaterState = 1  % Initially on

% Indicator states (for detecting temperature crossings)
Upwards = T_avg - 50.5[degC]    % Upward crossing
Downwards = T_avg - 49.5[degC]  % Downward crossing

% Implicit events
When Upwards > 0: HeaterState = 0 (turn off heating)
When Downwards < 0: HeaterState = 1 (turn on heating)

Mesh Configuration

Structured mesh with refinement in boundary layer regions to capture the thin thermal and velocity boundary layers near walls:

Water domain: Mapped mesh, maximum element size 1 mm
Boundary layers: 6 boundary layer mesh layers
Solid regions: Free triangular mesh

Solver Settings

Time Stepping

Time range: 0 to 3 minutes
Time step: 3 seconds
Relative tolerance: 1×10⁻³
Absolute tolerance: 2.5×10⁻⁵

Nonlinear Solver

Fully Coupled solver with direct solver (PARDISO). This approach handles the strong coupling between momentum, energy, and species transport equations effectively.

Results and Discussion

Let’s see what our simulation reveals!

Temperature Field Evolution

Temperature field animation showing the evolution

The simulation starts with water at a uniform 49°C—just below the thermostat’s lower threshold. As soon as the heater kicks in, we see the physics spring to life.

The temperature distribution shows characteristic features of Rayleigh-Bénard convection:

Hot thermal plumes rising from the heating plate
Cold water descending along the cooler side walls
Gradual homogenization due to convective mixing

Natural Convection Flow Field

Velocity field streamline plot showing recirculation zones

The streamlines tell a beautiful story of organized chaos. Buoyancy drives clear recirculation patterns—convection cells that act like tiny heat pumps, efficiently moving energy from the hot bottom to the cooler top.

Notice how:

Central upward flow of heated water
Peripheral downward flow of cooled water
Multiple recirculation cells may form depending on the Rayleigh number

Species Concentration Distribution

Here’s where things get practical. Imagine you’ve dropped a tea bag or dissolved some minerals. The species transport results show how convection dramatically speeds up mixing:

Initial concentration hotspots quickly spread throughout the kettle
Convection outpaces diffusion by orders of magnitude
Concentration contours follow the flow streamlines like tracer particles

Average Temperature Curve

The thermostat does its job beautifully! The average water temperature oscillates between 49.5°C and 50.5°C—exactly as designed. The sawtooth pattern reveals:

Rapid initial heating when heater is on
Gradual cooling when heater is off due to heat loss to the environment
Hysteresis prevents rapid on-off switching

Heat Flux Analysis

Boundary	Heat Flux Direction	Description
Heating plate (when on)	Inflow (+Q_in)	Primary heat source
Sidewall	Outflow (-Q_wall)	Convective loss to ambient
Water surface	Outflow (-Q_top)	Convective and radiative loss

At quasi-steady state, energy in equals energy out: $Q_{in} \approx Q_{wall} + Q_{top}$

Summary: The Complete Workflow

Build Geometry: Create 2D axisymmetric model of kettle with water domain and solid regions
Define Materials: Assign temperature-dependent water properties and solid material properties
Set Boundary Conditions: Configure heating plate, sidewall convection, and water surface conditions
Configure Events: Implement thermostat hysteresis control logic
Generate Mesh: Create structured mesh with boundary layer refinement
Run Simulation: Execute time-dependent solver with appropriate tolerances
Analyze Results: Examine temperature fields, flow patterns, and species distribution

Conclusions and Future Work

What We Learned

This simulation journey revealed several interesting insights:

Convection dominates: The high Rayleigh number confirms that natural convection—not conduction—is the primary heat transfer mechanism. The beautiful recirculation patterns in our streamline plots prove it.
Smart thermostat design: The hysteresis strategy works exactly as intended, keeping water within ±0.5°C of the target while avoiding the rapid on-off cycling that would plague a simple threshold controller.
Where the heat goes: Most heat escapes through the sidewalls rather than the water surface—a consequence of the larger contact area with the glass.
Mixing for free: Natural convection dramatically accelerates species transport, dispersing dissolved substances far faster than diffusion alone.

Model Limitations

No model is perfect. Ours makes several simplifying assumptions:

2D axisymmetric geometry (can’t capture 3D asymmetric effects)
No evaporation at the water surface (important at higher temperatures)
Spatially uniform heating plate temperature
Perfect thermal contact at material interfaces

Where to Go From Here

Several exciting extensions await:

Phase change: Add boiling and evaporation models for high-temperature operation
Realistic heating: Model the actual Joule heating distribution in the heating element
Full 3D: Capture asymmetric effects from handles and spouts
Smarter control: Implement PID control for improved energy efficiency
Surface effects: Include Marangoni convection at the water-air interface

References

COMSOL Multiphysics Documentation - Heat Transfer Module
COMSOL Application Gallery - Free Convection in a Water Glass
Incropera, F.P., DeWitt, D.P., Bergman, T.L., & Lavine, A.S. Fundamentals of Heat and Mass Transfer (7th ed.). Wiley.
COMSOL Application Gallery - Nonisothermal Flow Examples

Appendix

Click to expand complete MATLAB script

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
%% COMSOL Model: Electric Kettle Water Heating with Thermostat Control
% This script models the heating and keep-warm process of an electric kettle.
% The thermostat maintains water temperature between 49.5 deg C and 50.5 deg C
% by switching the heater ON/OFF automatically.
% 
% Physics Interfaces:
%   - Events (ev): Thermostat state control
%   - Laminar Flow (spf): Fluid dynamics
%   - Diluted Species (tds): Concentration transport
%   - Heat Transfer in Fluids (ht): Thermal analysis
% 
% Multiphysics Couplings:
%   - Non-Isothermal Flow (nitf1)
%   - Reacting Flow, Diluted Species (rfd1)

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

model = ModelUtil.create('Model');
model.title('Electric Kettle Water Heating with Keep-Warm Function');
model.description('This model simulates the heating and temperature maintenance process of an electric kettle using thermostat control. The heater switches ON when water temperature drops below 49.5 deg C and OFF when it exceeds 50.5 deg C.');

%% ==================== Parameters ====================
% Kettle geometry parameters
model.param.set('R_top', '7[cm]', 'Kettle top radius');
model.param.set('R_bottom', '7.5[cm]', 'Kettle bottom radius');
model.param.set('H_kettle', '16[cm]', 'Kettle height');
model.param.set('T_wall', '0.13[cm]', 'Kettle wall thickness');
model.param.set('T_bottom', '0.3[cm]', 'Kettle bottom thickness');
model.param.set('L_wall', 'sqrt((R_top-R_bottom)^2+H_kettle^2)', 'Outer wall slant length');
model.param.set('rho_ref', '1000[kg/m^3]', 'Reference density');

%% ==================== Component & Geometry ====================
model.component.create('comp1', true);
model.component('comp1').geom.create('geom1', 2);
model.component('comp1').geom('geom1').axisymmetric(true);

% Water domain
model.component('comp1').geom('geom1').create('pol1', 'Polygon');
model.component('comp1').geom('geom1').feature('pol1').set('source', 'table');
model.component('comp1').geom('geom1').feature('pol1').set('table', { ...
    '0' 'T_bottom'; ...
    'R_bottom-T_wall' 'T_bottom'; ...
    'R_top-T_wall' 'H_kettle'; ...
    '0' 'H_kettle'});

% Kettle bottom (heating plate)
model.component('comp1').geom('geom1').create('pol4', 'Polygon');
model.component('comp1').geom('geom1').feature('pol4').set('source', 'table');
model.component('comp1').geom('geom1').feature('pol4').set('table', { ...
    '0' 'T_bottom'; ...
    '0' '0'; ...
    'R_bottom-T_wall+T_bottom/H_kettle*sqrt(L_wall^2-H_kettle^2)' '0'; ...
    'R_bottom-T_wall' 'T_bottom'});

% Kettle wall
model.component('comp1').geom('geom1').create('pol2', 'Polygon');
model.component('comp1').geom('geom1').feature('pol2').set('source', 'table');
model.component('comp1').geom('geom1').feature('pol2').set('table', { ...
    '0' 'H_kettle'; ...
    '0' '0'; ...
    'R_bottom' '0'; ...
    'R_top' 'H_kettle'});

model.component('comp1').geom('geom1').run;
model.component('comp1').geom('geom1').run('fin');

%% ==================== Selections ====================
model.component('comp1').selection.create('sel1', 'Explicit');
model.component('comp1').selection('sel1').set([2]);
model.component('comp1').selection('sel1').label('Water');

%% ==================== Variables ====================
model.component('comp1').variable.create('var1');
model.component('comp1').variable('var1').set('T_heater_on', '80[degC]', 'Heater ON temperature');
model.component('comp1').variable('var1').set('T_heater_off', '25[degC]', 'Heater OFF temperature');
model.component('comp1').variable('var1').set('T_target', '50[degC]', 'Target temperature');
model.component('comp1').variable('var1').set('T_avg', 'aveop1(T)', 'Average water temperature');

%% ==================== Views ====================
model.view.create('view2', 3);
model.component('comp1').view('view1').set('showmaterial', true);
model.component('comp1').view('view1').axis.set('xmin', -0.158);
model.component('comp1').view('view1').axis.set('xmax', 0.233);
model.component('comp1').view('view1').axis.set('ymin', -0.008);
model.component('comp1').view('view1').axis.set('ymax', 0.168);
model.view('view2').set('showgrid', false);

%% ==================== Materials ====================
% --- Material 1: Silica Glass (Kettle Wall) ---
% Only thermal properties needed for heat transfer
model.component('comp1').material.create('mat1', 'Common');
model.component('comp1').material('mat1').selection.set([3]);
model.component('comp1').material('mat1').label('Silica glass');
model.component('comp1').material('mat1').set('family', 'custom');
model.component('comp1').material('mat1').set('transparency', 0.5);
model.component('comp1').material('mat1').propertyGroup('def').set('heatcapacity', '703[J/(kg*K)]');
model.component('comp1').material('mat1').propertyGroup('def').set('density', '2203[kg/m^3]');
model.component('comp1').material('mat1').propertyGroup('def').set('thermalconductivity', '1.38[W/(m*K)]');

% --- Material 2: Water ---
% Temperature-dependent properties for heat transfer and fluid flow
model.component('comp1').material.create('mat2', 'Common');
model.component('comp1').material('mat2').selection.named('sel1');
model.component('comp1').material('mat2').label('Water, liquid');
model.component('comp1').material('mat2').set('family', 'water');

% Define temperature-dependent functions
model.component('comp1').material('mat2').propertyGroup('def').func.create('eta', 'Piecewise');
model.component('comp1').material('mat2').propertyGroup('def').func.create('Cp', 'Piecewise');
model.component('comp1').material('mat2').propertyGroup('def').func.create('rho', 'Piecewise');
model.component('comp1').material('mat2').propertyGroup('def').func.create('k', 'Piecewise');
model.component('comp1').material('mat2').propertyGroup('def').func.create('an1', 'Analytic');

% Dynamic viscosity eta(T)
model.component('comp1').material('mat2').propertyGroup('def').func('eta').set('arg', 'T');
model.component('comp1').material('mat2').propertyGroup('def').func('eta').set('pieces', { ...
    '273.15' '413.15' '1.3799566804-0.021224019151*T^1+1.3604562827E-4*T^2-4.6454090319E-7*T^3+8.9042735735E-10*T^4-9.0790692686E-13*T^5+3.8457331488E-16*T^6'; ...
    '413.15' '553.75' '0.00401235783-2.10746715E-5*T^1+3.85772275E-8*T^2-2.39730284E-11*T^3'});
model.component('comp1').material('mat2').propertyGroup('def').func('eta').set('argunit', 'K');
model.component('comp1').material('mat2').propertyGroup('def').func('eta').set('fununit', 'Pa*s');

% Heat capacity Cp(T)
model.component('comp1').material('mat2').propertyGroup('def').func('Cp').set('arg', 'T');
model.component('comp1').material('mat2').propertyGroup('def').func('Cp').set('pieces', { ...
    '273.15' '553.75' '12010.1471-80.4072879*T^1+0.309866854*T^2-5.38186884E-4*T^3+3.62536437E-7*T^4'});
model.component('comp1').material('mat2').propertyGroup('def').func('Cp').set('argunit', 'K');
model.component('comp1').material('mat2').propertyGroup('def').func('Cp').set('fununit', 'J/(kg*K)');

% Density rho(T)
model.component('comp1').material('mat2').propertyGroup('def').func('rho').set('arg', 'T');
model.component('comp1').material('mat2').propertyGroup('def').func('rho').set('smooth', 'contd1');
model.component('comp1').material('mat2').propertyGroup('def').func('rho').set('pieces', { ...
    '273.15' '293.15' '0.000063092789034*T^3-0.060367639882855*T^2+18.9229382407066*T-950.704055329848'; ...
    '293.15' '373.15' '0.000010335053319*T^3-0.013395065634452*T^2+4.969288832655160*T+432.257114008512'});
model.component('comp1').material('mat2').propertyGroup('def').func('rho').set('argunit', 'K');
model.component('comp1').material('mat2').propertyGroup('def').func('rho').set('fununit', 'kg/m^3');

% Thermal conductivity k(T)
model.component('comp1').material('mat2').propertyGroup('def').func('k').set('arg', 'T');
model.component('comp1').material('mat2').propertyGroup('def').func('k').set('pieces', { ...
    '273.15' '553.75' '-0.869083936+0.00894880345*T^1-1.58366345E-5*T^2+7.97543259E-9*T^3'});
model.component('comp1').material('mat2').propertyGroup('def').func('k').set('argunit', 'K');
model.component('comp1').material('mat2').propertyGroup('def').func('k').set('fununit', 'W/(m*K)');

% Thermal expansion coefficient alpha_p(T) - needed for Boussinesq approximation
model.component('comp1').material('mat2').propertyGroup('def').func('an1').set('funcname', 'alpha_p');
model.component('comp1').material('mat2').propertyGroup('def').func('an1').set('expr', '-1/rho(T)*d(rho(T),T)');
model.component('comp1').material('mat2').propertyGroup('def').func('an1').set('args', {'T'});
model.component('comp1').material('mat2').propertyGroup('def').func('an1').set('fununit', '1/K');
model.component('comp1').material('mat2').propertyGroup('def').func('an1').set('argunit', {'K'});
model.component('comp1').material('mat2').propertyGroup('def').func('an1').set('plotargs', {'T' '273.15' '373.15'});

% Set material properties
model.component('comp1').material('mat2').propertyGroup('def').set('thermalexpansioncoefficient', {'alpha_p(T)' '0' '0' '0' 'alpha_p(T)' '0' '0' '0' 'alpha_p(T)'});
model.component('comp1').material('mat2').propertyGroup('def').set('dynamicviscosity', 'eta(T)');
model.component('comp1').material('mat2').propertyGroup('def').set('heatcapacity', 'Cp(T)');
model.component('comp1').material('mat2').propertyGroup('def').set('density', 'rho(T)');
model.component('comp1').material('mat2').propertyGroup('def').set('thermalconductivity', {'k(T)' '0' '0' '0' 'k(T)' '0' '0' '0' 'k(T)'});
model.component('comp1').material('mat2').propertyGroup('def').addInput('temperature');

% --- Material 3: Structural Steel (Heating Plate) ---
% Only thermal properties needed
model.component('comp1').material.create('mat3', 'Common');
model.component('comp1').material('mat3').selection.set([1]);
model.component('comp1').material('mat3').label('Structural steel');
model.component('comp1').material('mat3').set('family', 'steel');
model.component('comp1').material('mat3').propertyGroup('def').set('heatcapacity', '475[J/(kg*K)]');
model.component('comp1').material('mat3').propertyGroup('def').set('thermalconductivity', '44.5[W/(m*K)]');
model.component('comp1').material('mat3').propertyGroup('def').set('density', '7850[kg/m^3]');

%% ==================== Coupling Operators ====================
model.component('comp1').cpl.create('aveop1', 'Average');
model.component('comp1').cpl('aveop1').selection.set([2]);
model.component('comp1').cpl('aveop1').set('axisym', true);

%% ==================== Ambient Properties ====================
model.component('comp1').common.create('ampr1', 'AmbientProperties');
model.component('comp1').common('ampr1').set('T_amb', '25[degC]');

%% ==================== Physics: Events (Thermostat Control) ====================
model.component('comp1').physics.create('ev', 'Events', 'geom1');

% Discrete state: Heater ON/OFF
model.component('comp1').physics('ev').create('ds1', 'DiscreteStates', -1);
model.component('comp1').physics('ev').feature('ds1').set('dim', 'HeaterState');
model.component('comp1').physics('ev').feature('ds1').set('dimInit', 1);
model.component('comp1').physics('ev').feature('ds1').set('dimDescr', 'State of the Heater');

% Indicator states for temperature crossing detection
model.component('comp1').physics('ev').create('is1', 'IndicatorStates', -1);
model.component('comp1').physics('ev').feature('is1').set('indDim', {'Upwards'; 'Downwards'});
model.component('comp1').physics('ev').feature('is1').set('g', {'T_avg-50.5[degC]'; 'T_avg-49.5[degC]'});
model.component('comp1').physics('ev').feature('is1').set('dimInit', [0; 0]);
model.component('comp1').physics('ev').feature('is1').set('dimDescr', {'Upwards crossing'; 'Downwards crossing'});

% Implicit event: Switch OFF heater when T > 50.5°C
model.component('comp1').physics('ev').create('impl1', 'ImplicitEvent', -1);
model.component('comp1').physics('ev').feature('impl1').set('condition', 'Upwards>0');
model.component('comp1').physics('ev').feature('impl1').set('reInitName', {'HeaterState'; ''});
model.component('comp1').physics('ev').feature('impl1').set('reInitValue', [0; 0]);
model.component('comp1').physics('ev').feature('impl1').label('Switch off heater');

% Implicit event: Switch ON heater when T < 49.5°C
model.component('comp1').physics('ev').create('impl2', 'ImplicitEvent', -1);
model.component('comp1').physics('ev').feature('impl2').set('condition', 'Downwards<0');
model.component('comp1').physics('ev').feature('impl2').set('reInitName', 'HeaterState');
model.component('comp1').physics('ev').feature('impl2').set('reInitValue', 1);
model.component('comp1').physics('ev').feature('impl2').label('Switch on heater');

%% ==================== Physics: Laminar Flow (spf) ====================
model.component('comp1').physics.create('spf', 'LaminarFlow', 'geom1');
model.component('comp1').physics('spf').selection.named('sel1');

% Pressure point constraint
model.component('comp1').physics('spf').create('prpc1', 'PressurePointConstraint', 0);
model.component('comp1').physics('spf').feature('prpc1').selection.set([4]);

% Slip wall boundary condition (water surface)
model.component('comp1').physics('spf').create('wallbc2', 'WallBC', 1);
model.component('comp1').physics('spf').feature('wallbc2').selection.set([5]);
model.component('comp1').physics('spf').feature('wallbc2').set('BoundaryCondition', 'Slip');

% Physics properties
model.component('comp1').physics('spf').prop('ShapeProperty').set('order_fluid', 2);
model.component('comp1').physics('spf').prop('PhysicalModelProperty').set('IncludeGravity', true);
model.component('comp1').physics('spf').prop('PhysicalModelProperty').set('rref', {'R_top-T_wall'; '0'; 'H_kettle'});
model.component('comp1').physics('spf').prop('AdvancedSettingProperty').set('UsePseudoTime', true);

%% ==================== Physics: Diluted Species (tds) ====================
model.component('comp1').physics.create('tds', 'DilutedSpecies', 'geom1');
model.component('comp1').physics('tds').selection.set([2]);
model.component('comp1').physics('tds').feature('init1').set('initc', '1*(z>15[cm])');

%% ==================== Physics: Heat Transfer (ht) ====================
model.component('comp1').physics.create('ht', 'HeatTransferInFluids', 'geom1');

% Solid heat transfer for glass and heating plate
model.component('comp1').physics('ht').create('solid1', 'SolidHeatTransferModel', 2);
model.component('comp1').physics('ht').feature('solid1').selection.set([1 3]);

% Initial temperature
model.component('comp1').physics('ht').feature('init1').set('Tinit', '49[degC]');

% Temperature boundary (heating plate - thermostat controlled)
model.component('comp1').physics('ht').create('temp1', 'TemperatureBoundary', 1);
model.component('comp1').physics('ht').feature('temp1').selection.set([2]);
model.component('comp1').physics('ht').feature('temp1').set('T0', 'T_heater_off + (T_heater_on - T_heater_off) * HeaterState');

% Heat flux boundary 1: Side wall (natural convection)
model.component('comp1').physics('ht').create('hf1', 'HeatFluxBoundary', 1);
model.component('comp1').physics('ht').feature('hf1').selection.set([8]);
model.component('comp1').physics('ht').feature('hf1').set('HeatFluxType', 'ConvectiveHeatFlux');
model.component('comp1').physics('ht').feature('hf1').set('HeatTransferCoefficientType', 'ExtNaturalConvection');
model.component('comp1').physics('ht').feature('hf1').set('ExtNaturalConvectionType', 'InclinedWall');
model.component('comp1').physics('ht').feature('hf1').set('Lwall', 'L_wall');
model.component('comp1').physics('ht').feature('hf1').set('phi_tilt', 'acos(H_kettle/L_wall)');
model.component('comp1').physics('ht').feature('hf1').set('pA_src', 'root.comp1.ampr1.p_amb');
model.component('comp1').physics('ht').feature('hf1').set('Text_src', 'root.comp1.ampr1.T_amb');

% Heat flux boundary 2: Top surface (convection)
model.component('comp1').physics('ht').create('hf2', 'HeatFluxBoundary', 1);
model.component('comp1').physics('ht').feature('hf2').selection.set([5 6]);
model.component('comp1').physics('ht').feature('hf2').set('HeatFluxType', 'ConvectiveHeatFlux');
model.component('comp1').physics('ht').feature('hf2').set('h', '2[W/(m^2*K)]');
model.component('comp1').physics('ht').feature('hf2').set('Text_src', 'root.comp1.ampr1.T_amb');

%% ==================== Multiphysics Couplings ====================
model.component('comp1').multiphysics.create('nitf1', 'NonIsothermalFlow', 2);
model.component('comp1').multiphysics('nitf1').set('BoussinesqApproximation', true);

model.component('comp1').multiphysics.create('rfd1', 'ReactingFlowDS', 2);

%% ==================== Mesh ====================
model.component('comp1').mesh.create('mesh1');

% Mapped mesh for water domain
model.component('comp1').mesh('mesh1').create('map1', 'Map');
model.component('comp1').mesh('mesh1').feature('map1').selection.geom('geom1', 2);
model.component('comp1').mesh('mesh1').feature('map1').selection.set([2]);
model.component('comp1').mesh('mesh1').feature('map1').set('smoothcontrol', true);
model.component('comp1').mesh('mesh1').feature('map1').create('size1', 'Size');
model.component('comp1').mesh('mesh1').feature('map1').feature('size1').set('custom', 'on');
model.component('comp1').mesh('mesh1').feature('map1').feature('size1').set('hmax', 0.001);
model.component('comp1').mesh('mesh1').feature('map1').feature('size1').set('hmaxactive', true);

% Boundary layer mesh
model.component('comp1').mesh('mesh1').create('bl1', 'BndLayer');
model.component('comp1').mesh('mesh1').feature('bl1').create('blp', 'BndLayerProp');
model.component('comp1').mesh('mesh1').feature('bl1').feature('blp').selection.set([4 7]);
model.component('comp1').mesh('mesh1').feature('bl1').feature('blp').set('blnlayers', 6);

% Free triangular mesh for remaining domains
model.component('comp1').mesh('mesh1').create('ftri1', 'FreeTri');
model.component('comp1').mesh('mesh1').feature('ftri1').set('smoothcontrol', true);

model.component('comp1').mesh('mesh1').run;

%% ==================== Study ====================
model.study.create('std1');
model.study('std1').create('time', 'Transient');
model.study('std1').feature('time').set('tlist', 'range(0,3[s],3[min])');
model.study('std1').feature('time').set('usertol', true);
model.study('std1').feature('time').set('rtol', '1e-3');
model.study('std1').feature('time').set('plot', true);
model.study('std1').feature('time').set('probesel', 'none');

%% ==================== Solver Configuration ====================
model.sol.create('sol1');
model.sol('sol1').attach('std1');

model.sol('sol1').create('st1', 'StudyStep');
model.sol('sol1').feature('st1').label('Compile Equations: Time Dependent');

model.sol('sol1').create('v1', 'Variables');
model.sol('sol1').feature('v1').label('Dependent Variables 1.1');
model.sol('sol1').feature('v1').set('clist', {'{range(0[s], 3[s], 3[min])}' '0.18[s]'});

model.sol('sol1').create('t1', 'Time');
model.sol('sol1').feature('t1').label('Time-Dependent Solver 1.1');
model.sol('sol1').feature('t1').set('tlist', 'range(0,3[s],3[min])');
model.sol('sol1').feature('t1').set('rtol', '1e-3');
model.sol('sol1').feature('t1').set('atolglobalvaluemethod', 'manual');
model.sol('sol1').feature('t1').set('atolglobal', '2.5e-5');
model.sol('sol1').feature('t1').set('maxorder', 2);
model.sol('sol1').feature('t1').set('stabcntrl', true);
model.sol('sol1').feature('t1').set('bwinitstepfrac', 0.01);
model.sol('sol1').feature('t1').set('estrat', 'exclude');
model.sol('sol1').feature('t1').set('plot', true);


% Direct solver (PARDISO) - must be created before fc1 references it
model.sol('sol1').feature('t1').create('d1', 'Direct');
model.sol('sol1').feature('t1').feature('d1').label('Direct, nonisothermal flow (nitf1) (Merged)');
model.sol('sol1').feature('t1').feature('d1').set('linsolver', 'pardiso');
model.sol('sol1').feature('t1').feature('d1').set('pivotperturb', 1.0E-13);

% Iterative solver with multigrid
model.sol('sol1').feature('t1').create('i1', 'Iterative');
model.sol('sol1').feature('t1').feature('i1').label('AMG, nonisothermal flow (nitf1)');
model.sol('sol1').feature('t1').feature('i1').set('maxlinit', 100);
model.sol('sol1').feature('t1').feature('i1').set('rhob', 20);

model.sol('sol1').feature('t1').feature('i1').create('mg1', 'Multigrid');
model.sol('sol1').feature('t1').feature('i1').feature('mg1').label('Multigrid 1.1');
model.sol('sol1').feature('t1').feature('i1').feature('mg1').set('prefun', 'saamg');
model.sol('sol1').feature('t1').feature('i1').feature('mg1').set('maxcoarsedof', 80000);
model.sol('sol1').feature('t1').feature('i1').feature('mg1').set('saamgcompwise', true);

model.sol('sol1').feature('t1').feature('i1').feature('mg1').feature('pr').create('sc1', 'SCGS');
model.sol('sol1').feature('t1').feature('i1').feature('mg1').feature('pr').feature('sc1').set('linesweeptype', 'ssor');
model.sol('sol1').feature('t1').feature('i1').feature('mg1').feature('pr').feature('sc1').set('iter', 0);

model.sol('sol1').feature('t1').feature('i1').feature('mg1').feature('po').create('sc1', 'SCGS');
model.sol('sol1').feature('t1').feature('i1').feature('mg1').feature('po').feature('sc1').set('linesweeptype', 'ssor');
model.sol('sol1').feature('t1').feature('i1').feature('mg1').feature('po').feature('sc1').set('iter', 1);

model.sol('sol1').feature('t1').feature('i1').feature('mg1').feature('cs').create('d1', 'Direct');
model.sol('sol1').feature('t1').feature('i1').feature('mg1').feature('cs').feature('d1').set('linsolver', 'pardiso');
model.sol('sol1').feature('t1').feature('i1').feature('mg1').feature('cs').feature('d1').set('pivotperturb', 1.0E-13);

% Fully coupled solver - now d1 exists
model.sol('sol1').feature('t1').create('fc1', 'FullyCoupled');
model.sol('sol1').feature('t1').feature('fc1').label('Fully Coupled 1.1');
model.sol('sol1').feature('t1').feature('fc1').set('linsolver', 'd1');
model.sol('sol1').feature('t1').feature('fc1').set('maxiter', 8);
model.sol('sol1').feature('t1').feature('fc1').set('ntolfact', 0.5);
model.sol('sol1').feature('t1').feature('fc1').set('jtech', 'once');
model.sol('sol1').feature('t1').feature('fc1').set('stabacc', 'aacc');
model.sol('sol1').feature('t1').feature('fc1').set('aaccdim', 5);
model.sol('sol1').feature('t1').feature('fc1').set('aaccmix', 0.9);
model.sol('sol1').feature('t1').feature('fc1').set('aaccdelay', 1);
model.sol('sol1').feature('t1').feature('fc1').set('damp', '0.9');

model.sol('sol1').feature('t1').feature('aDef').set('cachepattern', true);
model.sol('sol1').feature('t1').feature.remove('fcDef');

%% ==================== Run Solver ====================
% model.study('std1').runNoGen;

%% ==================== Save Model ====================
mphsave(model, 'electric_kettle_thermostat.mph');