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:

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:

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

and

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:

0.01*sum(sum(if((k!=0)||(l!=0),((k^2+l^2)^(-b/2))*g1(k,l)*cos(2*pi*(k*s1+l*s2)+u1(k,l)),0),k,-N,N),l,-N,N)

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

0.01*sum(sum(sum(if((k!=0)||(l!=0)||(m!=0),((k^2+l^2+m^2)^(-b/2))*g1(k,l,m)*cos(2*pi*(k*x+l*y+m*z)+u1(k,l,m)),0),k,-N,N),l,-N,N),m,-N,N)

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.

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.

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.

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

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*.

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

0.01*sum(sum(sum(if((k!=0)||(l!=0)||(m!=0),((k^2+l^2+m^2)^(-b/2))*g1(k,l,m)*cos(2*pi*(k*x+l*y+m*z)+u1(k,l,m)),0),k,-N,N),l,-N,N),m,-N,N)

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.

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.

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

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.

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.

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.

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/m^{3} 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.

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.

Now, let’s generate a default mesh.

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.

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.

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.

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*.

Then, change the material data accordingly.

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.

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.

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.

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 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.

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.

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.

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.

##### Rehema Ndeda

June 24, 2017Thank you for the post. It has been quite useful.

##### Yu Zhang

April 2, 2018Hi, 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.

Yu

##### Bjorn Sjodin

April 3, 2018Hi 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:

https://www.comsol.com/blogs/how-to-generate-random-surfaces-in-comsol-multiphysics/

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,

Bjorn

##### Yu Zhang

April 4, 2018Hi, 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.

Yu

##### Bjorn Sjodin

April 4, 2018Hi Yu,

1)

The relationship between the continuous variables x, y, and z and the discrete k,l, and m are:

x=k/(2K+1)

y=l/(2L+1)

z=m/(2M+1)

To compare with the usual way of writing discrete transforms you may need to look at the previous blog post again

https://www.comsol.com/blogs/how-to-generate-random-surfaces-in-comsol-multiphysics/

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

##### Bjorn Sjodin

April 4, 2018Hi 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.

Bjorn

##### Yu Zhang

April 4, 2018Hi, 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,

Yu

## Comments (7)