How to Generate Randomized Inhomogeneous Material Data

June 20, 2017

You can generate and visualize randomized material data with specified statistical properties determined by a spectral density distribution by using the tools available under the Results node in the COMSOL Multiphysics® software. In this blog post, we show examples that are quite general and have potential uses in many application areas, including heat transfer, structural mechanics, subsurface flow, and more.

Synthesizing Randomized Inhomogeneous Material Data

We have already discussed how to generate random-looking geometric surfaces by using sum and if operators in combination with uniform and Gaussian random distribution functions. The idea is that by summing up a set of spatially varying waves with careful choices of amplitudes and phase angles, we can mimic the type of randomness frequently found in natural materials and many natural phenomena in general.

To generate synthetic material data in 2D, we can use the same double-sum expression that we used to create randomized CAD surface data in the previous blog post:

f(x,y)=\sum_{k=-K}^{K} \sum_{l=-L}^{L} a(k,l) cos(2 \pi(kx+ly)+\phi(k,l))

Note that this sum could also be used to generate random data for use in boundary conditions on surfaces in a 3D model.

For the 3D volume case, we will need triple sums:

f(x,y,z)=\sum_{k=-K}^{K} \sum_{l=-L}^{L} \sum_{m=-M}^{M} a(k,l,m) cos(2 \pi(kx+ly+mz)+\phi(k,l,m))

The frequency-dependent amplitudes will take their values from a random distribution according to

a(k,l) = g(k,l) h(k,l)=\frac{g(k,l)}{\vert k^2+l^2\vert^{\beta}}=\frac{g(k,l)}{(k^2+l^2)^{\frac{\beta}{2}}}


a(k,l,m) = g(k,l,m) h(k,l,m)=\frac{g(k,l,m)}{\vert k^2+l^2+m^2 \vert^{\beta}}=\frac{g(k,l,m)}{(k^2+l^2+m^2)^{\frac{\beta}{2}}}

for the 2D and 3D cases, respectively.

The functions g(k,l) and g(k,l,m) have a random Gaussian, or normal, distribution and h(k,l) and h(k,l,m) are frequency-dependent amplitude functions with values that taper off for higher frequencies in accordance with the spectral exponent β. The higher the value of the spectral exponent, the smoother the generated data will be. For a variety of reasons, many natural processes have this property, which is characterized by slower variations that are more dominant than fast.

The phase angles φ are sampled from a function u. The function has a uniform random distribution that is between -\frac{\pi}{2} and \frac{\pi}{2}:
φ(k,l) = u(k,l) and φ(k,l,m) = u(k,l,m)

The sums run over a set of discrete frequencies. More nonzero terms in a sum imply a larger number of higher-frequency contributions, resulting in data that contains finer details. The maximum frequencies in each direction are determined by the integers K, L, and M, respectively. These types of sums are closely linked to discrete inverse cosine transforms. They essentially correspond to an inverse cosine transform of the amplitude coefficients g(k,l) and g(k,l,m), with some additional manipulations of the phase angles. For details, see the previous blog post on how to generate random surfaces in COMSOL Multiphysics.

Expressions for Double and Triple Summation

In COMSOL Multiphysics, the following double-sum expression can be entered in various edit fields, such as 2D material properties or 3D boundary conditions:


We can use a similar expression for the triple sum, which can be used for 3D material data, loads, sources, and sinks:


where we have set K = L = M = N.

For more details about the underlying theory and syntax used here, see the blog post mentioned above.

Working with models that contain triple sums is computationally quite expensive. It is more efficient to first generate the data and export it to file and then import it again as an interpolation function, perhaps in a separate model. This interpolation function can then be used in a variety of ways, as we will explain later. Alternatively, you can use external software to generate the data by means of inverse FFT.

Let’s now take a look at how to generate 3D material data.

Generating a Randomized Interpolation Table in COMSOL Multiphysics®

Creating a 3D volume matrix of random data is surprisingly easy. It amounts to creating a couple of random functions, some parameters, a Grid data set, and an Export node.

Start by creating a random function for the amplitudes with 3 input arguments based on a normal, or Gaussian, distribution. This corresponds to the function g(k,l,m) in the mathematical description above. In this case, we arbitrarily use the default settings for a random function with the mean value set to 0 and the standard deviation set to 1.

A screenshot of the settings for the Gaussian Random function in COMSOL Multiphysics.

Next, we create a random function for the phase angles with 3 input arguments based on a uniform distribution between -\frac{\pi}{2} and \frac{\pi}{2} corresponding to the function u(k,l,m) above.

A screenshot of the settings for the Uniform Random function in COMSOL Multiphysics.

Now, create a data set of the type Grid 3D that references the random functions as a source. We need this data set to give an evaluation context to the triple-sum expression that we will define later in the Export node.

A screenshot of the settings for Grid 3D in COMSOL Multiphysics.

We will use two results parameters, N and b, for the spatial frequency resolution and spectral exponent, respectively.

A screenshot of the Parameters settings in COMSOL Multiphysics.

To make it easier to work with the large data sets that are generated, you can turn off the Automatic update of plots option. This setting is available in the Settings window of the Results node. Turning it off avoids recomputing the expressions each time you click on a plot group under Results.

A cropped screenshot of the COMSOL Multiphysics GUI showing the Automatic update of plots option.

To visualize the data before exporting to file, add a Slice plot and type (or paste) the expression:


A screenshot of the COMSOL Multiphysics GUI featuring a Slice plot in the Graphics window.

To export the data, add a Data node under Export and type in the same expression as for the Slice plot above. In the Settings window of the Data node, make sure to set the data set to Grid 3D and to specify a file name that the data will be written to. Here, we can let the points be evaluated in a way that is independent of the Grid 3D data set.

For the setting Points to evaluate in, select Regular grid. For Data format, select Grid. You can freely choose the number of x, y, and z points to evaluate in. In the figure below, these points have each been set to 50. Note that data generation corresponding to numbers higher than 50 may take a very long time to generate. For a 50x50x50 grid, we already get 125,000 data points.

A screenshot showing the Data settings in COMSOL Multiphysics.

The text file that is generated and exported can now be imported to a new file for the purpose of setting up a physics analysis where we use the generated data in material properties. This can be done for any type of physics, including electromagnetics, structural mechanics, acoustics, CFD, heat transfer, and chemical analysis. By using the COMSOL® API in methods or applications, for example, such export and import operations can be automated and set in the context of a for-loop in order to generate statistics over a larger sample set. In this example, we only generate one set of data.

Using the Generated Data in a Heat Transfer Analysis

To illustrate how this type of data can be used, let’s create a test model of the simplest possible kind, based on a heat transfer analysis.

Start by creating a new 3D model using a Heat Transfer in Solids interface.

A screenshot showing the geometry settings for a 3D example model in COMSOL Multiphysics.

Now, import the data from file as an Interpolation function. This function will be available under the Global Definitions node.

A screenshot showing the Interpolation function settings in COMSOL Multiphysics.

The Interpolation function is given the name cloud and can later be accessed using expressions like cloud(x,y,z).

To make unit handling easy, when using this interpolation function, we will set the input argument units to m and the function unit to 1. The Function unit corresponds to the unit of f(x,y,z)=cloud(x,y,z) and setting it to 1 makes it dimensionless.

A screenshot showing interpolation, extrapolation, and unit settings in COMSOL Multiphysics.

To keep things simple, let’s use a Block geometry object that matches the imported data exactly, with the corner at the origin and sides at 1. This corresponds to the size and position of the Grid 3D data set used earlier for generating the data.

A COMSOL screenshot showing the Block geometry object settings.

For a “real” case, you can instead import or create a CAD geometry, which can be used to truncate the interpolation function in a suitable way. This truncation of data is automatic in COMSOL Multiphysics. The figure below shows such an interpolation of randomized data over a CAD model of a wrench. When evaluating over an arbitrary geometry, it can be useful to scale the coordinate values in the triple-sum. In the wrench example, instead of k*x+l*y+m*z, as in the expressions above, the scaled expression k*(x/0.05)+l*(y/0.05)+m*(z/0.05) is used.

A COMSOL model of a wrench using a scaled expression.

This type of irregular material data may have uses in statistical modeling of materials such as those found in additive manufacturing, where perfect material homogeneity of a 3D-printed component may not always be possible to achieve. The data can be used for any type of material property, such as conductivity, permeability, and elasticity properties, to name a few.

Getting back to our unit cube example, we now add a Blank Material node. We will, somewhat arbitrarily, set the Density to 2000 kg/m3 and the Heat capacity to 1 J/kg/K. Since we are performing a stationary analysis, the Heat capacity is irrelevant. The Thermal conductivity is set to the expression 1+2[W/m/K]*cloud(x,y,z). We can see from the earlier Slice plot visualization that the values for the interpolation table are roughly between -0.2 and 0.2. This means that this expression will generate an interesting spatial distribution of thermal conductivity values between about 0.6 and 1.4.

A screenshot of the COMSOL software with the thermal conductivity settings highlighted.

The coefficient 2[W/m/K] is used to assign a consistent unit to the expression. The constant 1 will be automatically converted to the correct unit: [W/m/K].

Let’s define some simple boundary conditions. Set the temperature at the top surface to 393.15[K] and the bottom surface to 293.15[K], corresponding to a 100-K temperature difference.

A demonstration of defining simple boundary conditions in a COMSOL model.

Now, let’s generate a default mesh.

An illustration of generating a default mesh in a COMSOL model.

COMSOL Multiphysics will automatically interpolate values such as those from the material properties to this unstructured mesh from the imported interpolation function. Alternatively, we could generate a swept mesh with hexahedral elements of the same size as the original data, 50x50x50. Such a representation would be more “true” to the original data.

You can experiment with different element orders, such as linear and quadratic types. Unless you use a very fine mesh that “oversamples” the data, the results will depend somewhat on the element order.

Running the Study will produce a couple of temperature plots, the second of which is an Isosurface plot.

An example of running a Study in COMSOL Multiphysics to generate temperature plots.

Notice how the Isosurface plot looks a little bit jagged, which is due to the underlying irregularity of the material data. We can create another Slice plot to yet again visualize the data. This time, we do so under the guise of thermal conductivity by using the variable ht.kmean, which equals the expression 1+2[W/m/K]*cloud(x,y,z) defined earlier.

An example of creating a Slice plot to visualize COMSOL model data.

Here, the data is sampled at a lower sampling density than the original interpolation function, since we used the default mesh with the Element size set to the default Normal setting. Successively refining this unstructured mesh will ultimately sample the data at more or less the same level of detail as the original synthesized data.

As mentioned earlier, the approach used here for heat transfer is applicable to virtually any other type of simulation. For example, in a porous media flow simulation, the randomized quantity would be permeability rather than thermal conductivity. In the case of porous media flow, a more advanced random distribution may be needed, but let’s save that discussion for a future blog post.

Binary Data and Percolation Effects

We can also use the synthesized data in a different way: by using Boolean expressions to convert it to binary data. This method can be used for simulating two or more materials where the material interface is randomized and the material properties change abruptly from one point to another. COMSOL Multiphysics will automatically handle the sharp interpolations needed for this case.

The following picture shows a visualization of the Boolean expression cloud(x,y,z)>-0.03, which evaluates to 1 at points where the inequality is true and 0 at the other points.

A visualization of a certain Boolean expression in COMSOL Multiphysics.

To get a nicer plot, you can set the resolution of the Slice plot to Extra fine. This setting is available in the Quality section of the Settings window for the Slice plot.

We would now like to use this type of binary information in a simulation. It can be interesting, for example, to use it in a heat transfer simulation to see the so-called percolation effects. For certain threshold values, you get a large connected component in the material so that the entire slab of material starts conducting much more efficiently.

To try this, change the expression for the thermal conductivity to 1-0.9[W/m/K]*(cloud(x,y,z)>thold), where thold is a global parameter. Start by defining thold in Parameters under Global Definitions.

A cropped screenshot of the Model Builder with the Parameters node highlighted.

Then, change the material data accordingly.

A cropped screenshot of the Material Contents section in the COMSOL GUI.

For each point in space, the Thermal conductivity will, in a binary fashion, evaluate to 1 or 0.1, depending on the value of the inequality.

Now, let’s see how different values of this Boolean threshold will affect the simulation. For this purpose, run a parametric sweep over the parameter thold from -0.2 to 0.2.

A screenshot showing the Parametric Sweep settings in COMSOL Multiphysics.

Add a Surface Integration node under Derived Values to integrate the total heat flux that goes through one of the surfaces. This will be given by the surface integral of -ht.ntflux or +ht.ntflux, depending on if you are integrating over the top or bottom surface. In the figure below, we used the top surface.

A screenshot showing the Surface Integration settings in COMSOL Multiphysics.

The resulting Table plot shows the amount of heat power transferred (in watts). We can see that for threshold values around 0, the conductivity rises quickly from a low value to a high value. This is due to the sudden appearance of one or more large connected components where the expression 1-0.9[W/m/K]*(cloud(x,y,z)>thold) evaluates to 1.

A Table plot of transferred heat power for an example model in COMSOL Multiphysics.

The figures below show a Volume plot with a Filter attribute for three threshold values around 0. The filter shows the parts of the domain for which cloud(x,y,z)<thold corresponds to the locations of higher conductivity.

We can see from these figures how the highly conductive parts start connecting for the higher threshold values.

The corresponding Filter settings are shown in the figure below.

A screenshot showing the Filter settings in COMSOL Multiphysics.

A similar type of percolation effect, seen here for binary data, is also happening for the continuous data case shown earlier. However, when using binary data, the effects are easier to see.

Point Cloud Visualization

Finally, let’s look at an alternative way of visualizing this type of random data. We will visualize the data set using a large number of randomized points (or rather, small spheres) and let the radius and color of the points vary according to the interpolation function cloud(x,y,z). In addition, we will only allow the points to be visualized for positive values of cloud(x,y,z). This technique will allow us to “see inside the data” in a way that is difficult to achieve using other methods. Note that this visualization technique works for any type of data, including real measured data.

Start by generating three random variables with uniform distribution, with the Range set to 1 and the Mean value set to 0.5.

A screenshot showing the settings for one random variable in COMSOL Multiphysics.

To generate this type of plot, we use a Scatter Volume plot type. This is available by right-clicking a 3D Plot Group and selecting More Plots > Scatter Volume.

In the Settings window of the Scatter Volume plot, set the expression for the X-, Y-, and Z-components: rn1(x), rn2(x), and rn3(x), respectively. Here, we are using the x-coordinate in an unusual way, in that we are using it merely as a long vector of arbitrary values.

Next, in the Evaluation Points section, set the Number of points for the X grid points to 100,000; 1,000,000; or more, depending on how many points your computer can handle. Set each of the Y grid point and Z grid point values to 1. This is a trick for getting a long vector of values that we can feed into the random functions in order to generate a lot of random points within the unit cube.

A COMSOL plot of random material data in grayscale.

To make the plot appear as in the above figure, go to the Radius section and set Expression to cloud(rn1(x),rn2(x),rn3(x)) and the Radius scale factor to 0.3. In addition, in the Color section, set the Expression to cloud(rn1(x),rn2(x),rn3(x)) and the Color table to GrayScale.

One additional noteworthy fact about this plot: negative values will be ignored. This helps our visualization, since roughly half of the generated data is negative and we can more easily see through the data and get an intuitive feel for the variations. This method will only work for a rectangular block. To instead generate this type of plot over an arbitrary CAD geometry, you can use the Particle Tracing Module, which allows you to generate random points inside any type of CAD model.

A point cloud.

By the way, a similar-looking plot can be achieved in a 2D model by simply creating a 2D Surface plot using a double-sum expression, as shown in the figure above.

Comments (20)

Leave a Comment
Log In | Registration
Rehema Ndeda
Rehema Ndeda
June 24, 2017

Thank you for the post. It has been quite useful.

Yu Zhang
Yu Zhang
April 2, 2018

Hi, Bjorn:

Thank you for sharing us this method to generated random field. Now I am going to generate a Gaussian field (mean 0, and variance 1) which is spatial correlated, and found that the fast Fourier Transformation is commonly used to generate such a random field. I am not sure whether I can use your method to create a spatial correlated Gaussian field which has mean 0 and variance 1. In you method, is the generated f(x,y) just the Gaussian field I wanted to create?

Thank you for your time and help.


Bjorn Sjodin
Bjorn Sjodin
April 3, 2018

Hi Yu,

Thanks for your comment. The method described here is more or less equivalent to an FFT-based method. However, this method will of course not scale as well as an FFT-based method. You will find a comparison between the two methods in a previous blog posting:
in the section called
Relationship to Discrete Cosine and Fourier Transforms
There you will also find a book reference for the FFT-based method.

Some benefits with the method used in this blog post are: a) ease-of-use b) mesh independence, you can refine the mesh and get consistent results without redefining the random process.

Best regards,

Yu Zhang
Yu Zhang
April 4, 2018

Hi, Bjorn:

Thank you so much for your detailed reply, which helped me get more insights about the this method. However, I am still have some confusions about it. Could you please give me a further explanation? Any help would be highly appreciated.

(1) How does the autocorrelation function (like Gaussian type rho(h)=exp(-(h/b)^2), where b is the correlation length) relate to the parameters (like k, l, m) in your method?

(2) After the transform, does the generated random field still have standard Gaussian distribution (mean 0, and variance 1)?

Thank you so much for your time and help.


Bjorn Sjodin
Bjorn Sjodin
April 4, 2018

Hi Yu,
The relationship between the continuous variables x, y, and z and the discrete k,l, and m are:
To compare with the usual way of writing discrete transforms you may need to look at the previous blog post again
where we also describe the index shift that you can perform in the discrete Fourier transform.

When all that is sorted out you need to make sure you have a fine enough mesh to resolve the individual cos (or sin) functions for the highest frequencies.
We need to remember here that we are performing these calculations on an unstructured triangulated/tetrahedrized mesh instead of a structured grid as in the discrete case.
With finite elements (such as in COMSOL) those cos functions are approximated with a piecewise polynomial on each element. By refining the mesh you should in theory get as good autocorrelation as your computer can handle a fine mesh. The finer the mesh the closer to the “exact” autocorrelation you will get.
2) The mean will be 0 but the variance calculation will not be 1. For example, the arbitrary factor 0.01 in front of the series is controlling the variance and the more terms you include the greater the variance. I haven’t looked for a literature reference for what the variance is in terms of this type of truncated Fourier series, or tried to derive it, but it would certainly be interesting to see such an expression. Perhaps some reader would know the answer to this question before any of us have looked it up. If I can find a nice expression for this I will make sure to post it here.


Bjorn Sjodin
Bjorn Sjodin
May 8, 2020 COMSOL Employee

There is an update and answer to the question 2 above about how to get the variance from this type of function.

The files that are available for download as part of the blog post now contains a file random_flat_surface_roughness.mph (made with version 5.5) . This is a geometry part file and you can, for example, copy it to the part library at:
C:\Program Files\COMSOL\COMSOL55\Multiphysics\parts\COMSOL_Multiphysics\Random_Surfaces
This part will allow you to specify the surface roughness by an amplitude scaling factor (the “old” way) or the surface RMS height (Sq). The method used to give the surface RMS height demonstrates how to control the standard deviation or, more or less equivalently, variance of this type of data and is also applicable to this posting on 3D random data. You just need a triple sum instead of a double sum. The file shows how to scale the amplitude to get a function that is proportional to the standard deviation.


Bjorn Sjodin
Bjorn Sjodin
April 4, 2018

Hi again Yu,

I think for understanding the correlation for Fourier synthesis methods, like this one, case you can read
The Science of Fractal Images, Editors: Peitgen, Heinz-Otto, Saupe, Dietmar. Eds.
which is the book referenced in the first blog posting. Describing this here would be too much to fit in this comment field and more of text book chapter character.

Then you also have the discretization error described in the above comment coming from the finite element method, so only when the mesh element size goes to zero you would get perfect agreement with the “exact” theory of text books.


Yu Zhang
Yu Zhang
April 4, 2018

Hi, Bjorn,

Thanks a lot for your generous help and detailed explanations. I will study the book systematically. If I have some new findings, I will post it here. Thank you.

Best regards,


March 23, 2020

Hi Bjorn,
Thank you for the post ! It’s was very useful ! However, I’m wondering if the method is still valid for the case when I already have a distribution of the properties ? (it means actually it’s no longer be something randomly distributed in the spatial domain). In this case, how could we randomly distribute the properties (with knowing distribution) over a surface or volume ?
Could we replace the gaussian distribution g1 by a distribution of the parameter we want to use ?


Bjorn Sjodin
Bjorn Sjodin
March 23, 2020 COMSOL Employee

Hi Tu,

I’m glad you liked the post. Do you mean that you have a known distribution function for your random properties? In that case, do you have a mathematical description of this distribution?


March 23, 2020

Thank Bjorn for your quick reply.
Indeed, as you said, I know the distribution of my properties (as a look-up table, but I can try to fit with an explicit function if needed).
To be more specific, for ex when I used the li-ion battery physics, I’d like to load the particle size distribution instead of a mean radius, and so the same for the conductivity, the diffusivity, etc…

Bjorn Sjodin
Bjorn Sjodin
March 24, 2020 COMSOL Employee

Hi Tu,

Without knowing the details of your application, I think you can use techniques similar to those in this blog post to approximate distribution functions based on Normal (or Gaussian) distribution functions. In this posting an approximate probability distribution is generated by a sum of Normal distribution functions. In a similar way you can approximate a wide range of probability distributions based on COMSOL’s built-in random function. I suggest you search the web for references on such techniques.


March 24, 2020

Thank you for the suggestion. I’ve found this thread , that can actually solve my problem.
For the coefficient 0.01, it’s arbitrary taken I suppose ? And if I don’t want to have a negative values In my results, what should I change in this case ?
Thank you,


Bjorn Sjodin
Bjorn Sjodin
March 27, 2020 COMSOL Employee

Hi Tu,

For this specific problem and for further discussions I suggest you contact our technical support team.

Best regards,

Eduardo Moreno
Eduardo Moreno
June 19, 2020

Do you have and mph example for the case of Randomized Inhomogeneous Material Data for 2D axial symmetric?
best regards
Eduardo Moreno

Bjorn Sjodin
Bjorn Sjodin
June 22, 2020 COMSOL Employee

Hi Eduardo,

We don’t have any examples of axial symmetry yet. Could you tell me a bit more about your application? This may be something we could cover by a future blog post.

Best regards,
Bjorn Sjodin

Diletta Giuntini
Diletta Giuntini
August 11, 2020

Hello Bjorn,
Thank you for this post, very useful.
I implemented this approach to generate an initial condition that is randomly distributed in space (3D). However, when I then use the cloud function into my problem’s geometry (an elongated bar), the random distribution does not occupy my entire domain, but only a small portion – even if I scale the expression as you suggested in the wrench expression above. Any insight, or should I perhaps share my model with the technical support?

Bjorn Sjodin
Bjorn Sjodin
August 11, 2020 COMSOL Employee

Hi Diletta,

I’m glad you found it useful. Probably some translation of the random data is also needed in your case. Usually you can accomplish this by replacing expressions in x, y, and z with corresponding expressions in x-x0, y-y0, and z-z0 where x0, y0, z0 is some kind of representative center of the geometry. In addition, you may need scaling and even potentially rotation on top of that. If you need further help don’t hesitate contacting technical support.

Best regards,

Rukshan Azoor
Rukshan Azoor
August 25, 2020

Hi Bjorn,

This and your previous article have been very useful. I just wanted to share another method that I eventually adopted to generate spatially correlated random field realisations to be input into COMSOL (feel free to remove this comment if not relevant).
I used a Python library to generate a spatially correlated Gaussian random field with mean 0 and variance 1 (similar to what Yu required above), so I can transform it to the required material property. I have described the method and shared the code here:
The code can be used in a Jupyter notebook and the spatial correlation lengths, size of domain and the resolution can be changed. The output is a point cloud text-file that can be imported directly into COMSOL as an interpolation function as you have done.

Also, is there an easy method that you could suggest to automate the process of importing data, running and saving results for multiple realisations?


Bjorn Sjodin
Bjorn Sjodin
August 26, 2020 COMSOL Employee

Hi Rukshan,

I’m glad you liked the articles. Automating this is possible by using the COMSOL API for use with Java, for example using the Method Editor in the Application Builder, or with MATLAB, using the LiveLink for MATLAB product.

I hope this helps,