Using MATLAB® Functions in Your COMSOL Multiphysics® Models
Fanny Littmarck | April 18, 2014
Did you know that you can use MATLAB® functions in your COMSOL Multiphysics® models? Well, you can, and in this video tutorial we will show you how, using LiveLink™ for MATLAB®.
Video: LiveLink™ for MATLAB® Demonstration
- Shown in the video: MEMS Heating Circuit
- Another example for using MATLAB functions: Temperature Distribution in a Vacuum Flask
In this tutorial, we will learn how to model the joule heating and thermal expansion in a MEMS heating circuit. The device consists of a thin, nichrome layer pattern on a thicker glass layer. Electric current conduction through the nichrome layer produces joule heating. Most of the material properties are assumed to be constant, except the thermal conductivity of glass.
A MATLAB® function is created that reads in temperature-dependent thermal conductivity data and interpolates to provide the conductivity value as a function of temperature at any point within the glass layer. The effect of spatial inhomogeneity in the material property is incorporated by adding a random component to the conductivity, which varies as a function of the Y-coordinate of the glass layer.
In order to set up the model in COMSOL Desktop®, we will create the geometry, set up the physics, assign material properties, mesh, and run the model.
Once we have run the basic model in the COMSOL Desktop, we will start the LiveLink™ for MATLAB® and load the COMSOL® model. We will then perform data fitting to obtain the parameters for temperature-dependent electrical conductivity and specify the limits of thickness of glass and nichrome layers. These parameters will be updated in the model. The model will be solved within a nested for-loop that will allow us to investigate the design space.
So, let’s create the model in COMSOL Desktop.
We will start by selecting a 3D space dimension and then, from the Add Physics menu, we’re going to select our three different physics interfaces of interest.
We will go to “Structural Mechanics”, select and add “Thermal Stress”. Then, we’ll go to “AC/DC” and select and add “Electric Currents, Shell”.
Then, we’ll go back to “Structural Mechanics” once again and we will add the Shell or the Structural Shell interface.
Let us go ahead and get into the geometry building steps to save some time. I’m going to get it from a file where these steps are already available, and let’s go ahead and build the geometry.
Now, we can go back to the global definitions and create a MATLAB function.
The name of the function will be thermal_conductivity and the arguments that we are going to send to this MATLAB function will be the temperature and the Y-coordinate.
The simulation will actually start with the Structural Shell interface. This will be assigned only to the nichrome layer, which is what we are modeling as a shell.
Now let’s go to the Electric Currents, Shell interface. Very similarly, we will assign the nichrome layer as the shell or the boundary of interest.
We will go ahead and then add a Boundary Heat Source boundary condition to the same surface and here we will basically connect the Thermal Stress interface to the surface losses produced by the Electric Currents, Shell interface.
So this is how we are telling the model that there’s some heat being generated on this boundary.
Let’s go to the Materials branch.
We’ll first go to the Material Browser, go to the built-in library and select “Silica glass”. This gets assigned to the glass layer.
Here, among all the material properties listed for silica glass, all these numbers are constant, and we’ll go ahead and change the value of “thermal conductivity”.
So, instead of that being a number, it will now be the MATLAB function thermal_conductivity and its value will be evaluated at temperature, T, and at any y-coordinate.
Once we’re done with setting up the material properties, let’s customize the mesh.
Once we are done with setting up the mesh, it is time to select “Compute”.
Once we have solved the model, we can create several plots to look at some of the quantities of interest. Let’s also look at the thermal conductivity of glass, which we obtain from the MATLAB function, thermal_conductivity. I mention that there is a randomness that we added to incorporate the effect of spatial inhomogeneity. We can see that randomness when we compare the conductivity distribution over space with the color map, here.
So, now that we have solved the basic model in COMSOL Desktop, let’s switch gears and take a look at the MATLAB graphical user interface that you get to see when you start LiveLink™ for MATLAB®.
We have a little script here and the first step that we perform is to load our model file that we created in COMSOL Desktop into the MATLAB workspace using the mphload command.
Once you do that, you get to see the model structure in the workspace.
Finally, let’s take a look at the main component in this script, where we use a nested for-loop to update the values of the glass layer thickness and the nichrome layer thickness, and solve the model repeatedly within this nested for-loop. We also create custom plots, which will only plot the temperature for the extreme design points.
We would also use a few specific functions, such as mphglobal and mphmean to obtain information such as the values of global quantities, like the “Maximum Interface Stress” and efficiency, and the average value of temperature on the bottom surface of the glass layer.
Once we obtain these values, we can assign them and store them in MATLAB variables.
We can also combine native MATLAB® functions with LiveLink™ functions. For example, mphinterp is a LiveLink™ for MATLAB® function that you can combine with the std, or standard deviation, function in MATLAB to obtain the standard deviation of temperature on the bottom surface of the glass layer, evaluated at a set of arbitrary coordinate points.
Once we have solved the problem for all these design cases, we can also create plots of interest. Such as the average temperature on the lower surface, the standard deviation of temperature on the lower surface, and surface plots of the Maximum Interface Stress and Heat Transfer Efficiency. All of these as functions of the thickness of the glass layer as well as the nichrome layer.
So, let’s go ahead and take a look at these plots. Here is our first plot, bunch of subplots actually, which shows the temperature at the different extremes of our design space. Minimum and maximum value of the glass thickness, and also minimum and maximum value of the nichrome layer thickness.
Let’s take a look at the plot of the average temperature on the lower surface of glass as a function of both glass layer thickness and the nichrome layer thickness. We see that the effect of changing the glass thickness is not so much on the average temperature, but there is some significant effect when we change the thickness of the nichrome layer. The thicker it is, the higher the temperature.
In the next figure, we are looking at the standard deviation of temperature as a function of, again, the glass thickness and the nichrome thickness — and we definitely see some nonlinear behavior here.
Let’s take a look at the surface plot where we can see the Maximum Interface Stress as a function of the two different layered thicknesses — and once again we see some nonlinearities. We see that the Maximum Interface Stress is obtained when we have a thicker glass layer and a much thinner nichrome layer.
MATLAB is a registered trademark of The MathWorks, Inc.