COMSOL 6.4 - Heating Circuit

Small heating circuits are used in many applications. For example, in manufacturing processes they heat up reactive fluids. Figure 1 illustrates a typical heating device for this model. The device consists of an electrically resistive layer deposited on a glass plate, which causes Joule heating when a voltage is applied to the circuit. The layer’s properties determine the amount of heat produced.

Figure 1: Geometry of a heating device.

In this particular model, you must observe three important design considerations:

•

Noninvasive heating

•

Minimal deflection of the heating device

•

Avoidance of overheating the process fluid

The heater must also work without failure. You achieve the first and second requirements by inserting a glass plate between the heating circuit and the fluid; it acts as a conducting separator. Glass is an ideal material for both of these purposes because it is nonreactive and has a low coefficient of thermal expansion.

You must also avoid overheating due to the risk of self-ignition of the reactive fluid stream. Ignition is also the main reason for separating the electrical circuit from direct contact with the fluid. The heating device is tailored for each application, making virtual prototyping very important for manufacturers.

For heating circuits in general, detachment of the resistive layer often determines the failure rate. This is caused by excessive thermally induced interfacial stresses. Once the layer is detached it gets locally overheated, which accelerates the detachment. Finally, in the worst case, the circuit might overheat and burn. From this perspective, it is also important to study the interfacial tension due to the differences in temperature as well as the different thermal-expansion coefficients of the resistive layer and the substrate. The geometric shape of the layer is a key parameter to design circuits for proper functioning. You can investigate all of the abovementioned aspects by modeling the circuit.

This multiphysics example simulates the electrical heat generation, the heat transfer, and the mechanical stresses and deformations of a heating circuit device. The model uses the Heat Transfer in Solids interface of the Heat Transfer Module in combination with the Electric Currents, Layered Shell interface from the AC/DC Module and the Solid Mechanics interface from the Structural Mechanics Module.

This model requires the AC/DC Module, Heat Transfer Module, and Structural Mechanics Module.

Model Definition

Figure 2 shows a drawing of the modeled heating circuit.

Figure 2: Drawing of the heating circuit deposited on a glass plate.

The device consists of a serpentine-shaped Nichrome resistive layer, 10 μm thick and 5 mm wide, deposited on a glass plate. At each end, it has a silver contact pad measuring 10 mm-by-10 mm-by-10 μm. When the circuit is in use, the deposited side of the glass plate is in contact with surrounding air, and the back side is in contact with the heated fluid. Assume that the edges and sides of the glass plate are thermally insulated.

Table 1 gives the resistor’s dimensions.

Table 1: Dimensions.
object	Length	Width	Thickness
Glass Plate	130 mm	80 mm	2 mm
Pads and Circuit	-	-	10 µm

During operation the resistive layer produces heat. Model the electrically generated heat using the Electric Currents, Layered Shell interface from the AC/DC Module. An electric potential of 12 V is applied to the pads. In the model, you achieve this effect by setting the potential at one edge of the first pad to 12 V and that of one edge of the other pad to 0 V.

To model the heat transfer in the thin conducting layer, use the Thin Layer feature from the Heat Transfer in Solids interface. The heat rate per unit area (measured in W/m2) produced inside the thin layer is given by

(1)

where QDC = J · E = σ|∇tV|2 (W/m3) is the power density. The generated heat appears as an inward heat flux at the surface of the glass plate.

At steady state, the resistive layer dissipates the heat it generates in two ways: on its up side to the surrounding air (at 293 K), and on its down side to the glass plate. The glass plate is similarly cooled in two ways: on its circuit side by air, and on its back side by a process fluid (353 K). You model the heat fluxes to the surroundings using heat transfer coefficients, h. For the heat transfer to air, h = 5 W/(m2·K), representing natural convection. On the glass plate’s back side, h = 20 W/(m2·K), representing convective heat transfer to the fluid. The sides of the glass plate are insulated.

The model simulates thermal expansion using static structural-mechanics analyses. It uses the Solid Mechanics interface for both the glass plate and the circuit layer. The stresses in the circuit layer are modeled using a thin layer approximation by adding a Thin Layer feature. The equations of these two interfaces are described in the Structural Mechanics Module User’s Guide. The stresses are set to zero at 293 K. The displacements in the Solid Mechanics interface are constrained in such a way that the plate can freely deform while its rigid body motions are suppressed through appropriate constraints.

Table 2 summarizes the material properties used in the model.

Table 2: Material properties.
Material	E [GPa]	ν	α [1/K]	k [W/(m·K)]	ρ [kg/m3]	Cp [J/(kg·K)]
Silver	83	0.37	1.89e-5	420	10500	230
Nichrome	213	0.33	1e-5	15	9000	20
Glass	73.1	0.17	5.5e-7	1.38	2203	703

Results and Discussion

Figure 3 shows the heat that the resistive layer generates.

Figure 3: Stationary heat generation in the resistive layer when 12 V is applied.

The highest heating power occurs at the inner corners of the curves due to the higher current density at these spots. The total generated heat, as calculated by integration, is approximately 13.8 W.

Figure 4 shows the temperature of the resistive layer and the glass plate at steady state.

Figure 4: Temperature distribution in the heating device at steady state.

The highest temperature is approximately 428 K, and it appears in the central section of the circuit layer. It is interesting to see that the differences in temperature between the fluid side and the circuit side of the glass plate are quite small because the plate is very thin. Using boundary integration, the integral heat flux on the fluid side evaluates to approximately 8.5 W. This means that the device transfers the majority of the heat it generates — 8.5 W out of 13.8 W — to the fluid, which is good from a design perspective, although the thermal resistance of the glass plate results in some losses.

The temperature rise also induces thermal stresses due the materials’ different coefficients of thermal expansion. As a result, mechanical stresses and deformations arise in the layer and in the glass plate. Figure 5 shows the equivalent stress distribution in the device and the resulting deformations. During operation, the glass plate bends toward the air side.

Figure 5: The thermally induced von Mises stress plotted with the deformation.

The highest equivalent stress, approximately 13 MPa, occurs at the inner corners of the curves of the Nichrome circuit. The yield stress for high quality glass is roughly 250 MPa, and for Nichrome it is 360 MPa. This means that the individual objects remain structurally intact for the simulated heating power loads.

You must also consider stresses in the interface between the resistive layer and the glass plate. Assume that the yield stress of the surface adhesion in the interface is in the region of 50 MPa — a value significantly lower than the yield stresses of the other materials in the device. If the equivalent stress increases above this value, the resistive layer locally detaches from the glass. Once it has detached, heat transfer is locally impeded, which can lead to overheating of the resistive layer and eventually cause the device to fail.

Figure 6 displays the effective forces acting on the adhesive layer during heater operation. As the figure shows, the device experiences a maximum interfacial stress that is an order of magnitude smaller than the yield stress. This means that the device is well designed in terms of adhesive stress.

Figure 6: The effective forces at the interface between the resistive layer and the glass plate.

Finally, study the warping of the device, that is, its deviation from a plane surface, shown in Figure 7.

Figure 7: Warping displacement on the fluid side of the glass plate.

The maximum deviation from being a planar surface, is approximately 50 μm. For high-precision applications, such as semiconductor processing, this might be a significant value that limits the device’s operating temperature.

ACDC_Module/Electromagnetic_Heating/heating_circuit

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Geometry 1

The interface includes and . In the volume, these two interfaces solve for temperature and displacement, respectively. In the shell representing the circuit, the temperature, the electric potential and displacement are solved by , , and interfaces, respectively.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
V_in	12[V]	12 V	Input voltage
d_layer	10[um]	1E-5 m	Layer thickness
sigma_silver	6.3e7[S/m]	6.3E7 S/m	Electric conductivity of silver
sigma_nichrome	9.3e5[S/m]	9.3E5 S/m	Electric conductivity of nichrome
T_air	20[degC]	293.15 K	Air temperature
h_air	5[W/(m^2*K)]	5 W/(m²·K)	Heat transfer film coefficient, air
T_fluid	353[K]	353 K	Fluid temperature
h_fluid	20[W/(m^2*K)]	20 W/(m²·K)	Heat transfer film coefficient, fluid

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Block 1 (blk1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 80.

In the text field, type 130.

In the text field, type 2.

Click

Work Plane 1 (wp1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2.

Click

Work Plane 1 (wp1) > Plane Geometry

Click the

button in the toolbar.

Work Plane 1 (wp1) > Square 1 (sq1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 10.

Locate the section. In the text field, type 7.

In the text field, type 10.

Click

Work Plane 1 (wp1) > Square 2 (sq2)

Right-click > > > > and choose .

In the window for , locate the section.

In the text field, type 30.

In the text field, type 8.

Click

Work Plane 1 (wp1) > Polygon 1 (pol1)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Click

Browse to the model’s Application Libraries folder and double-click the file heating_circuit_polygon.txt.

Click

Work Plane 1 (wp1) > Fillet 1 (fil1)

In the toolbar, click

On the object , select Points 2–8, 23–29, 34, 36, 37, 41, and 42 only.

It might be easier to select the points by using the window. To open this window, in the toolbar click and choose . (If you are running the cross-platform desktop, you find in the main menu.)

In the window for , locate the section.

In the text field, type 10.

Click

Work Plane 1 (wp1) > Fillet 2 (fil2)

In the toolbar, click

On the object , select Points 6–12, 26–31, 37, 40, 43, 46, 49, and 50 only.

In the window for , locate the section.

In the text field, type 5.

In the toolbar, click

Form Union (fin)

In the toolbar, click

The geometry should look like the figure below.

In the window, under > click .

Definitions

Add a selection that you can use later when applying boundary conditions and shell physics settings.

Circuit

In the toolbar, click

In the window for , type Circuit in the text field.

Locate the section. From the list, choose .

Select Boundaries 6–8 only.

Before creating the materials for use in this model, it is a good idea to specify which boundaries are to be modeled as conducting shells. Using this information, COMSOL Multiphysics can detect which material properties are needed.

Heat Transfer in Solids (ht)

Thin Layer 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose . In the Lth text field, type d_layer.

Locate the section. From the list, choose .

Electric Currents in Shells (ecis)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Lth text field, type d_layer.

Conductive Shell 1

In the window, under > click .

In the window for , locate the section.

From the εr list, choose .

Use in the interface.

Solid Mechanics (solid)

In the window, under click .

Thin Layer 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Lth text field, type d_layer.

To model the thermal effects in the thin layers, add a multiphysics coupling.

Multiphysics

Thermal Expansion, Thin Layer 1 (tetl1)

In the toolbar, click

and choose > .

Now set up the materials.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Materials

Silver Layer

In the window, under right-click and choose > .

In the window for , type Silver Layer in the text field.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Heat capacity at constant pressure	Cp	230[J/(kg*K)]	J/(kg·K)	Basic
Young’s modulus	E	83e9[Pa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.37	1	Young’s modulus and Poisson’s ratio
Density	rho	10500[kg/m^3]	kg/m³	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	420[W/(m*K)]	W/(m·K)	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	sigma_silver	S/m	Basic
Coefficient of thermal expansion	alpha_iso ; alphaii = alpha_iso, alphaij = 0	18.9e-6[1/K]	1/K	Basic
Thickness	lth	d_layer	m	Shell

Nichrome Layer

Right-click and choose > .

Select Boundary 7 only.

In the window for , type Nichrome Layer in the text field.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Heat capacity at constant pressure	Cp	20[J/(kg*K)]	J/(kg·K)	Basic
Young’s modulus	E	213e9[Pa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.33	1	Young’s modulus and Poisson’s ratio
Density	rho	9000[kg/m^3]	kg/m³	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	15[W/(m*K)]	W/(m·K)	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	sigma_nichrome	S/m	Basic
Coefficient of thermal expansion	alpha_iso ; alphaii = alpha_iso, alphaij = 0	10e-6[1/K]	1/K	Basic
Thickness	lth	d_layer	m	Shell

Electric Currents in Shells (ecis)

Ground 1

In the toolbar, click

and choose .

Select Edge 43 only.

Electric Potential 1

In the toolbar, click

and choose .

Select Edge 10 only.

In the window for , locate the section.

In the V0 text field, type V_in.

With the materials defined, set up the remaining physics of the model. In the next section, the resistive loss within the circuit is defined as a heat source for the thermal stress physics. The resistive loss is calculated automatically within the interface. Add the coupling feature to take the resistive loss into account.

Multiphysics

Electromagnetic Heating 1 (emh1)

In the toolbar, click

and choose > .

Heat Transfer in Solids (ht)

Heat Flux 1

In the toolbar, click

and choose .

Select Boundaries 4 and 6–8 only.

In the window for , locate the section.

From the list, choose .

In the h text field, type h_air.

In the Text text field, type T_air.

Heat Flux 2

In the toolbar, click

and choose .

Select Boundary 3 only.

In the window for , locate the section.

From the list, choose .

In the h text field, type h_fluid.

In the Text text field, type T_fluid.

Solid Mechanics (solid)

In order for the problem to be well posed, the glass plate must be constrained so that it does not have any possible rigid body translations or rotations. The constraints must be such that no stresses are induced by inhibited thermal expansion.

Rigid Motion Suppression 1

In the toolbar, click

and choose .

The absolute value of the maximum or minimum displacement is not important in itself, since it is just a function of how the rigid body constraints are applied. Instead, you want to see how much the boundary deviates from being planar. Use the feature to display that.

Warpage 1

In the toolbar, click

and choose .

Select Boundary 3 only.

In the window for , locate the section.

From the list, choose .

Locate the section. Click to select the

toggle button.

Select Point 1 only.

Locate the section. Click to select the

toggle button.

Select Point 3 only.

Locate the section. Click to select the

toggle button.

Select Point 63 only.

Mesh 1

Free Triangular 1

In the toolbar, click

and choose .

Select Boundaries 4 and 6–8 only.

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 2.

Swept 1

In the toolbar, click

Distribution 1

Right-click and choose .

In the window for , locate the section.

In the text field, type 3.

Click

Results

Preferred Units 1

In the toolbar, click

and choose .

In the window for , locate the section.

Select the checkbox.

Click

In the dialog, select > in the tree.

Click .

In the window for , locate the section.

In the table, enter the following settings:


Quantity	Unit	Preferred unit
Displacement	m	µm

Click

In the dialog, select > in the tree.

Click .

In the window for , locate the section.

In the table, enter the following settings:


Quantity	Unit	Preferred unit
Stress tensor	N/m^2	MPa

Study 1

In order to improve the solver’s performance, scale the degrees of freedom in to 1e-3.

Solution 1 (sol1)

In the toolbar, click

In the window, expand the node.

In the window, expand the > > > node, then click .

In the window for , locate the section.

From the list, choose .

In the text field, type 1e-3.

In the toolbar, click

Results

The first default plot shows the von Mises stress and the deformation of the device (Figure 5). The third default plot shows the temperature distribution in the volume of the full 3D geometry (Figure 4). Two plots are also generated to visualize the electric potential and the von Mises stress on the circuit layer.

Temperature (ht)

To generate Figure 3 follow the steps below.

Surface Losses

In the toolbar, click