How to Generate Randomized Inhomogeneous Material Data

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