To easily generate random-looking geometric surfaces, the COMSOL Multiphysics® software provides a powerful set of built-in functions and operators, such as functions for uniform and Gaussian random distributions and a very useful *sum* operator. In this blog post, we show you how to generate a randomized surface with what amounts to a “one liner” expression with detailed control of the constituent spatial frequency components that determine the nature of the surface’s roughness.

### Characterizing Surface Roughness

There are many ways to characterize a rough surface. One way is to use its approximate fractal dimension, which is a value between 2 and 3 for a surface. A surface of fractal dimension 2 is an ordinary, almost everywhere smooth surface; the value 2.5 represents a fairly rugged surface; and values close to 3 represent something that is close to “3D space filling”. Correspondingly, a curve of fractal dimension 1 is smooth almost everywhere, the value 1.5 represents a fairly rugged line, and values close to 2 represent something that is close to “2D space filling”.

*The range of fractal dimension values for curves going from 1 (left) to about 1.2 (center) and to 1.6 (right).*

Using a fractal dimension measure can be a useful approximation, but we need to remember that real surfaces aren’t fractal in nature over more than a few orders of magnitude of scale. Real surfaces have a spatial frequency “cutoff” due to their finite size and due to the fact that when “zooming in”, you will eventually hit some new type of microstructure behavior.

Another way of characterizing surface roughness is with respect to its spatial frequency content. This can be turned into a constructive method of synthesizing surface data by using a sum of trigonometric functions similar to a Fourier series expansion. Each term in such a sum represents a certain frequency of oscillation through space. This is the method that we will use here. Let’s quickly review the concepts of spatial frequencies and elementary wave shapes before moving on to trigonometric series.

### Spatial Frequencies

In physics, the frequency of oscillations over time occurs in mathematical expressions like

cos(2\pi f t)

where the unit of the frequency *f* is 1/s, also known as hertz or Hz.

Oscillations through space have a corresponding spatial frequency, as in the following expression, where we simply have replaced the time variable *t* by a spatial variable *x* and the time frequency *f* with the spatial frequency *v*.

cos(2\pi \nu x)

where the SI unit of the spatial frequency is 1/*m*.

Spatial frequencies are commonly represented by a wave number *k* = 2*πv*.

A related quantity is the wavelength \lambda=\frac{1}{\nu}, which is related to the frequency and wave number as follows:

k=2\pi \nu=\frac{2\pi}{\lambda}

There may be more than one dimension of space and, accordingly, there may be multiple spatial frequencies. In 2D, using Cartesian coordinates, we have:

cos(2\pi (\nu_x x + \nu_y y))=cos(\bf{k} \cdot \bf{x})

where \bf{k}=(k_x,k_y)=(2\pi \nu_x,2\pi\nu_y) is the wave vector and \bf{x}=(x,y).

The wave vector \bf{k} represents the direction of the wave.

### Elementary Waves

A rough surface *f*(*x*,*y*) can be seen as composed of many elementary waves of the form

cos(\bf{k} \cdot \bf{x}+\phi)

where *φ* is a phase angle.

The phase angle also makes it possible to express sine functions due to the relationship sin(\theta)=cos(\pi/2-\theta).

For a completely random surface, it should hold that the phase angle *φ* can take any value in, say, the interval 0 to *π* or –*π*/2 to *π*/2. When synthesizing elementary waves for a random surface, we can pick *φ* from a uniform random distribution in such an interval of length *π*, since we then allow for the expression cos(\phi) to span all possible values between -1 and +1. Note that there may be end-point or wrap-around effects if we choose an interval with a size bigger than *π*. This is due to the cosine function being its own mirror image in steps of *π*, according to cos(\pi-\theta)=-cos(\theta).

In order to get an efficient representation that can be used for simulations, we will only allow for a discrete set of spatial frequencies:

*ν _{x}* =

*m, ν*=

_{y}*n*

where *m* and *n* are integers.

Let’s consider a surface that is composed of elementary waves of the following form:

cos(\bf{k}_{mn} \cdot \bf{x}+\phi)= cos(2 \pi (mx+ny)+\phi) , \bf{k}_{mn}=2\pi(m,n)

By letting *m* and *n* take both positive and negative values with equal probabilities, we should be able to get a method of synthesizing a surface with no preferred direction of oscillations.

Note that, in this way, each wave direction is represented twice. For example, the direction (-2,-3) is the same as (2,3); (2,-1) is the same as (-2,1); and so on.

If we allow the spatial frequencies *m* and *n* to take values up to maximum integers *M* and *N*, respectively, then this corresponds to a high-frequency cutoff at:

\nu_{xmax}=M, \nu_{ymax}=N

Since we also allow for negative values, there are negative cutoffs at:

\nu_{xmin}=-M, \nu_{ymin}=-N

Having a spatial frequency cutoff at \nu_{xmax}=M in the *x* direction means that the shortest wavelength we can represent is \lambda_{xmin}=\frac{1}{M}, and similarly for the *y* direction, \lambda_{ymin}=\frac{1}{N}.

### Associated Amplitudes for Elementary Waves

Each elementary wave will have an associated amplitude so that each constituent wave component has the following form:

A_{mn}cos(\bf{k}_{mn} \cdot \bf{x}+\phi)

The final surface will be a sum over such wave components:

f(\bf{x})=\sum_{m,n}A_{mn}cos(\bf{k}_{mn} \cdot \bf{x}+\phi)

The simplest choice of amplitude would be to choose the coefficients *A _{mn}* from a uniform or perhaps Gaussian distribution. However, it turns out that this will not generate a particularly natural-looking surface. In nature, different processes, such as wearing and erosion, make it more likely that slow oscillations have a larger amplitude than fast ones. In the discrete case, this corresponds to the amplitudes tapering off according to some distribution:

A_{mn} =a(m,n) \sim h(m,n)=\frac{1}{\vert m^2+n^2\vert^{\beta}}=\frac{1}{(m^2+n^2)^{\frac{\beta}{2}}}

where the spectral exponent *β* indicates how quickly higher frequencies are attenuated. Following *The Science of Fractal Images* (Ref.1), the spectral exponent can be related to the fractal dimension of a surface, but only for an infinite series of waves covering arbitrarily high frequencies and only for certain ranges of the exponent. In practice, the amplitudes *a*(*m*,*n*) of our synthesized surface will be generated using a limited number of frequencies, multiplied with a random function *g*(*m*,*n*) having a Gaussian distribution:

*a*(*m*,*n*) = *g*(*m*,*n*)*h*(*m*,*n*)

A Gaussian, or normal, distribution is chosen to get a smooth but random variation in amplitudes with no limit on the magnitude.

The phase angles *φ* will be sampled from a function *u* with a uniform random distribution between –*π*/2 and *π*/2:

*φ*(*m*,*n*) = *u*(*m*,*n*)

#### Summing It Up

To represent our rough surface, we want to use the following double sum:

f(x,y)=\sum_{m=-M}^{M} \sum_{n=-N}^{N} a(m,n) cos(2 \pi(mx+ny)+\phi(m,n))

where *x* and *y* are spatial coordinates; *m* and *n* are spatial frequencies; *a*(*m*,*n*) are amplitudes; and *φ*(*m*,*n*) are phase angles. This expression is similar to a truncated Fourier series. Although the series is expressed in terms of cosine functions, the phase angles make it so this sum can express a quite general trigonometric series due to the angle sum rule:

cos(\alpha+\beta)=cos(\alpha)cos(\beta)-sin(\alpha)sin(\beta)

### Determining Periodicity

Due to its definition, the function *f*(*x*,*y*) will be periodic. In order to get a natural-looking surface, we should “cut out” a suitably small portion by letting *x* and *y* vary between some limited values; otherwise, the periodicity of the synthesized data will be apparent. What should these values be?

The overall periodicity will be determined by the slowest oscillations, which correspond to the spatial frequencies *m* = 1 or *n* = 1 in the *x* direction and *y* direction, respectively. This gives a period length of 1 in each direction.

We could generate the surface over a rectangle [*a*, *a* + 1] × [*b*, *b* + 1] or smaller in order to “avoid” the periodicity.

### Defining Parameters and Random Functions in COMSOL Multiphysics®

For the COMSOL Multiphysics implementation, start by defining a couple of parameters for the spatial frequency resolution and spectral exponent according to the following figure:

The amplitude generation will require a random function with a Gaussian distribution in two variables. This functionality is available under the *Global Definitions* node:

Here, the *Label* and *Function name* have been changed to *Gaussian Random* and *g1*, respectively. In addition, the *Number of arguments* is set to *2* instead of the default *1* and the *Distribution type* is set to *Normal*, which corresponds to a normal or Gaussian distribution.

In a similar way, for the phase angle, we need a uniform random function in the interval between –*π*/2 and *π*/2:

The *Label* is changed to *Uniform Random*, the *Function name* to *u1*, the *Number of arguments* to *2*, and the *Range* to *pi*.

You can optionally use random seeds to get the same surface each time you use the same input parameters.

### Defining the Parametric Surface

The next step is to add a *Parametric Surface* node under *Geometry* using a fairly lengthy *z*-coordinate expression, as follows:

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

where *x* = *s1* and *y* = *s2* vary between 0 and 1.

The factor 0.01 is used to scale the data in the *z* direction. Alternatively, this scaling factor can be absorbed into the amplitude coefficients.

*A parametric surface geometry feature is used to generate a synthesized random surface.*

Note that whenever you update any of the parameters or expressions for the Parametric Surface, you need to click the *Rebuild with Updated Functions* button in the *Advanced Settings* section of the Settings window.

This expression is a double-sum over the integer parameters *m* and *n* each running from –*N* to *N*. If we compare this to the mathematical discussion earlier, we can see that we have set *M* = *N*, resulting in a square surface patch. The term where *m* and *n* are simultaneously zero corresponds to an unwanted “DC” term and is eliminated from the sum by the *if* statement.

The syntax for the *sum*() operator is as follows:

*sum(expr,index,lower,upper)*

which evaluates a sum of a general expression *expr* for all indices *index* from *lower* to *upper*.

The syntax for the *if*() operator is as follows:

*if(cond,expr1,expr2)*

for which the conditional expression *cond* is evaluated to *expr1* or *expr2* depending on the value of the condition.

In this example, the resolution of the parametric surface has been increased by setting the *Maximum number of knots* to 100 (the default is 20). In addition, the *Relative tolerance* is relaxed to 1e-3 (the default is 1e-6). The underlying representation of the parametric surface is based on nonuniform rational B-splines (NURBS). More knots correspond to a finer resolution of the NURBS representation. The tolerance is increased, since we are not overly concerned about the approximation accuracy of the generated surface for this example.

By generating a mesh, we can get a useful visualization of the surface, as seen in the figure below.

*A meshed random surface.*

Note that N = 20 means that the fastest oscillations are 1/20 = 0.05 m, assuming SI units. The periodicity in the *x* and *y* directions can be seen by following the curves parallel to the *y*– and *x*-axes at *x* = 0, *x* = 1 and *y* = 0, *y* = 1; respectively.

To see the periodicity even more clearly, we can plot the surface on the square [0,2] × [0,2]:

*The periodicity of the surface on the square [0,2] × [0,2]. The surface height is represented by color.*

*Surfaces generated on the square [0,1] × [0,1] by superimposing 20 frequency components with amplitude spectral exponents β = 0.5, β = 1.0, β = 1.5, and β = 1.8, clockwise from the top-left image. The surface height is represented by color.*

### Using the Surface Data in Analyses

This type of randomly generated surface can, in COMSOL Multiphysics, be used in any kind of physics simulation context, including for electromagnetics, structural mechanics, acoustics, fluid, heat, or chemical analysis. The expression for the double sum is not limited for use in geometry modeling, but can also be used for material data, equation coefficients, boundary conditions, and more. Using methods, a large number of surface realizations can be used in a loop to gather statistics of the results.

By generalizing the double-sum to a triple-sum, you can synthesize 3D inhomogeneous material data. However, you have to be prepared for long and memory-intensive computations when performing triple-sums for 3D simulations.

*A fracture flow simulation based on synthetically generated fracture aperture data. The Rock Fracture Flow tutorial model is part of the COMSOL Multiphysics Application Library.*

*A generic thermal expansion analysis of two 1-centimeter-sized metal blocks with a material interface based on the parametric surface described in this blog post. The bottom material slab is aluminum and the top material slab is steel. The visualization shows the von Mises stress at the material interface and on the surface of the aluminum slab.*

### Relationship to Discrete Cosine and Fourier Transforms

The sum

f(x,y)=\sum_{m=-M}^{M} \sum_{n=-N}^{N} a(m,n) cos(2 \pi(mx+ny)+\phi(m,n))

is similar to a discrete cosine transform or to the real part of a discrete Fourier transform:

f_c(x,y)=\sum_{m=-M}^{M} \sum_{n=-N}^{N} F_c(m,n)e^{i(2 \pi(mx+ny))}

where the subscript *c* is used to indicate complex quantities and *x* and *y* now take discrete values. Here, the phase angle information is encoded in the complex Fourier coefficients.

Due to the definition of the discrete Fourier transform, we are allowed to perform a shift in index in order to generate the following more familiar form:

f_c(x,y)=\sum_{m=0}^{2M} \sum_{n=0}^{2N} F_c(m,n)e^{i(2 \pi(mx+ny))}

or by using discrete values:

f_c(k,l)=\sum_{m=0}^{2M} \sum_{n=0}^{2N} F_c(m,n)e^{i(2 \pi(m \frac{k}{2M+1}+n \frac{l}{2N+1}))}

More commonly, the discrete Fourier transform is indexed like this:

f_c(k,l)=\sum_{m=0}^{\mathfrak{M}-1} \sum_{n=0}^{\mathfrak{N}-1} F_c(m,n)e^{i(2 \pi(m \frac{k}{\mathfrak{M}}+n \frac{l}{\mathfrak{N}}))}

where

\mathfrak{M}=2M+1, \mathfrak{N}=2N+1.

Note that in order to generate real-valued data, the Fourier coefficients need to fulfill conjugate symmetry relationships in order to eliminate the imaginary-valued contributions from sine functions. Using a sum of cosine functions (i.e., a cosine transform) avoids this problem.

A fast way of generating a large number of Fourier coefficients is to use a fast cosine transform (FCT) or fast Fourier transform (FFT). This could be done in another program and then imported to the COMSOL Desktop® user interface as an interpolation table. The trigonometric interpolation method described above is slower, but has the advantage that it can be used directly on an unstructured mesh and is automatically refined by simply refining the mesh in the user interface.

For a description of using FFT for synthesizing surfaces, see Ref.1.

### 1D and Cylindrical Cases

Let’s conclude with a few interesting, special cases of random surface generation in COMSOL Multiphysics, including curves and cylinders.

#### Random Curve

In a 2D simulation, a random curve can be generated using the following expression:

0.01*sum(if((m!=0),((m^2)^(-b/2))*g1(m)*cos(2*pi*m*s+u1(m)),0),m,-N,N)

where g1 and u1 are 1D random functions.

Note that when generating a curve, the spectral exponent will have a lower value as compared to that of a surface for the “same level of randomness”.

*A randomized curve with spectral exponent 0.8.*

#### Random Polar Curve

A randomized curve in polar coordinates representing random deviations from a circle can be generated:

x=cos(2*pi*s)*(1+0.1*sum(if((m!=0),((m^2)^(-b/2))*g1(m)*cos(2*pi*m*s+u1(m)),0),m,-N,N))

y=sin(2*pi*s)*(1+0.1*sum(if((m!=0),((m^2)^(-b/2))*g1(m)*cos(2*pi*m*s+u1(m)),0),m,-N,N))

This corresponds to a parametric curve in 2D polar coordinates:

x=r(\phi) cos(\phi)

y=r(\phi) sin(\phi)

*A randomized polar curve with spectral exponent 0.8.*

#### Random Cylinder

A randomized cylinder in 3D can be generated using a parametric surface with parameters as follows:

x=cos(2*pi*s1)*(1+0.1*sum(sum(if((m!=0)||(n!=0),((m^2+n^2)^(-b/2))*g1(m,n)*cos(2*pi*(m*s1+n*s2)+u1(m,n)),0),m,-N,N),n,-N,N))

y=sin(2*pi*s1)*(1+0.1*sum(sum(if((m!=0)||(n!=0),((m^2+n^2)^(-b/2))*g1(m,n)*cos(2*pi*(m*s1+n*s2)+u1(m,n)),0),m,-N,N),n,-N,N))

z=s2*2*pi

where the parameters *s1* and *s2* vary between 0 and 1.

This corresponds to a parametric surface in cylindrical coordinates:

x=r(\phi,z) cos(\phi)

y=r(\phi,z) sin(\phi)

z=z

Such a single-piece random cylinder represents a type of self-intersecting surface that is not allowed in COMSOL Multiphysics. You can easily get around this by, for example, creating four surface patches corresponding to the parameter *s1* varying from 0 to 0.25, 0.25 to 0.5, 0.5 to 0.75, and 0.75 to 1.0. One such patch corresponds to a polar angle span of size \frac{\pi}{2}.

*A randomized tubular surface using polar coordinates.*

### New Part in COMSOL Multiphysics® Version 5.5

As of version 5.5, the Part Library for COMSOL Multiphysics has been extended with several new parts, including a part for creating a random flat surface.

*The* Random Flat Surface *part in the Part Library in COMSOL Multiphysics.*

Get the Model Files

### Reference

*The Science of Fractal Images*, Editors: Peitgen, Heinz-Otto, Saupe, Dietmar. Eds.

## FAQs

### How do I create a selection in COMSOL? ›

To add selection nodes, **right-click a Definitions node and choose from the Selections options as listed in Table 6-6**. You can also right-click the Geometry node and choose from Selections options similar to those in Table 6-6 for defining selections based on the geometry objects in the geometry sequence.

**How do I add material properties to COMSOL? ›**

The materials tab or in the model builder window by right-clicking the materials node once the new

**What is a free surface on COMSOL? ›**

The free liquid surface **corresponds to the phase boundary between the liquid and the gas and is represented on a fixed mesh**. The figure below shows the surface of two droplets in a channel, taken from the droplet breakup model in the Application Library of the add-on Microfluidics Module.

**What is Workplane in Comsol? ›**

A work plane is **a 2D plane oriented anywhere in the 3D space**.

**How do I Draw a polygon in Comsol? ›**

To create a polygon, **right-click a 2D Geometry node and select Polygon ( )**. For a 3D model, in the Geometry toolbar, from the More Primitives ( ) menu, select Polygon. You can also right-click the Geometry node to add this node from the context menu.

**What is moving mesh in Comsol? ›**

Moving Mesh Features. **Features added under Moving Mesh control the spatial frame**. They can be used to study both stationary states and time-dependent deformations where the geometry changes its shape due to motion of solid boundaries and deformation of solid domains.

**How do I write expressions in Comsol? ›**

**Click the Insert Expression button ( ) at the bottom of the section or press Ctrl+Space to insert a predefined expression or variable into the selected cell in the Value column**. A list of available Mathematical Functions such as trigonometric functions.

**What is explicit Comsol? ›**

Definitions / Explicit Selections

Explicit selections in the COMSOL Multiphysics^{®} software **can be used to create selections of individual geometric entities**. This helps simplify your model and streamline your workflow.

**How is geometry defined in Comsol? ›**

(03) Defining Geometry - COMSOL 4.2 Tutorial - YouTube

**How do you copy geometry in Comsol? ›**

To copy geometry features, you can also **right-click the geometry feature in the model tree (for example, Rectangle or Sphere) and select Copy ( )**. Then right-click the Geometry node and select Paste (for example, Paste Rectangle or Paste Sphere) ( ).

### How do you make a hexagon in Comsol? ›

How to create Hexagon shape easily in COMSOL Multiphysics - YouTube

**How do you Draw an ellipse in Comsol? ›**

Ellipse. **On the 2D Geometry toolbar Draw group, from the Circle ( ) menu, select Ellipse ( ) or Ellipse (Corner)( )**. Then draw the ellipse in the Graphics window.