  # Part 2: Modeling the Harmonic Excitations of Nonlinear Systems

August 11, 2016

Extending our discussion on modeling the harmonic excitations of linear systems, we will now shift focus to nonlinear systems. We will look at problems where the loading on the system has some sinusoidal components as well as cases where the material properties or loads and constraints depend directly on the solution. As you will see, COMSOL Multiphysics can address these apparently nonlinear cases with some very efficient solution algorithms. Let’s find out how.

### Loads with Time-Harmonic Components and Nonlinear Systems

In a previous blog post on modeling the harmonic excitations of linear systems, we identified two conditions under which it is possible to model the transient response of a system in the frequency domain. These two conditions are:

1. All time-varying loads and constraints on the system must vary sinusoidally at the same fixed frequency.
2. All loads, constraints, and material properties must be independent of the solution.

Today, we will look at relaxing both of these assumptions a bit. Let’s begin our discussion with a very common case of a system where the loads do not vary purely sinusoidally in time: the vibration of a tensioned guitar string. As we all know, when you increase the tension on a guitar string, the pitch (the fundamental frequency of vibration) increases. The strings of a guitar vibrate at different frequencies based on the string tension.

Let’s consider such a case in terms of the generic partial differential equation from our earlier blog post:

(1)

\begin{align}
M u_{tt} + C u_t + \nabla \cdot (-K \nabla u) = F_0+\tilde F sin(\omega t) &\text { on } \Omega \\
\mathbf{n} \cdot (K \nabla u) + Au = f _0 + \tilde f sin(\omega t ) &\text{ on } \Gamma_1 \\
u = g_0 + \tilde g sin ( \omega t) &\text{ on } \Gamma_2
\end{align}

where the body loads, F_0+\tilde F \sin(\omega t), boundary loads, f_0+\tilde f sin(\omega t), and constraints, g_0+\tilde g \sin(\omega t), are decomposed into a constant component and a sinusoidally time-varying component. As a structural engineer, you might refer to this as a decomposition into static and dynamic components. As an electrical engineer, you would describe these elements as the DC and AC components, respectively.

We can further assume that the dynamic components are much smaller in magnitude than the static components, that is: \tilde F \ll F_0, \tilde f \ll f_0, and \tilde g \ll g_0. If the dynamic loads are relatively small, then it is reasonable to assume that the dynamic loads will lead to small oscillations for the statically loaded case. That is, we are assuming that the solution will take the form: u = u_0 + \tilde u. This is also known as a harmonic perturbation, or a frequency-domain perturbation. In COMSOL Multiphysics, we refer to this as a Prestressed Analysis study for the structural mechanics physics interfaces; a Small-Signal Analysis study for the electromagnetics physics interfaces; and an AC Impedance study for the electrochemistry physics interfaces.

These study types all have the same underlying solution approach. They start out by ignoring the sinusoidally time-varying component and first solve the stationary problem:

(2)

\begin{align}
\nabla \cdot (-K \nabla u_0) = F_0 &\text { on } \Omega \\
\mathbf{n} \cdot (K \nabla u_0) + Au_0 = f _0 &\text{ on } \Gamma_1 \\
u_0 = g_0 &\text{ on } \Gamma_2
\end{align}

where, in fact, both the material properties and loads can be directly dependent on the solution. That is, K = K(u_0), F_0 = F_0(u), f_0 = f_0(u_0), and g_0 = g_0(u_0). This means that the stationary problem can actually be nonlinear, which doesn’t really complicate the analysis as long as the combination of the material nonlinearities, static loads, and constraints represent a well-posed problem. For a better understanding of how to address such nonlinear stationary problems, take a look at our Solver series.

After solving for the stationary response of the system, we can use the static solution, u_0, to evaluate the nonlinear material properties used by the frequency-domain perturbation form of the governing equation:

(3)

\begin{align}
-\omega^2 M(u_0) \tilde u + j \omega C(u_0) \tilde u + \nabla \cdot (-K(u_0) \nabla \tilde u) = \tilde F &\text { on } \Omega \\
\mathbf{n} \cdot (K(u_0) \nabla \tilde u) + A(u_0) \tilde u = \tilde f &\text{ on } \Gamma_1 \\
\tilde u = \tilde g &\text{ on } \Gamma_2
\end{align}

For structural problems, the static solution will also represent the deformed state of the system. The frequency-domain perturbation equations are formed on this deformed state, introducing a so-called geometric nonlinearity, in addition to the material nonlinearity.

The above governing equation can also be posed as an eigenfrequency problem by homogenizing all loads and constraints, that is, setting each of them to zero magnitude.

Our Application Gallery features a number of tutorials that demonstrate this solution approach. Here are some examples:

### What About Nonlinearities Due to Frequency-Domain Excitation?

Up to this point, the examples presented have all assumed that the response of the system due to the frequency-domain loading is linear. That is, we previously assumed that if any of the frequency-domain loads or constraints (\tilde F, \tilde f, \tilde g) are increased in magnitude by a scalar value, then the frequency-domain solution will increase in magnitude by the same scalar value. Of course, for such an assumption to hold, the loads, constraints, and material properties must all be independent of, \tilde u, the frequency-domain solution. There are, however, many cases where this is not true. Let’s find out how to address some of these situations.

We can start with a case where material properties are dependent on the cycle-averaged solution. Since the field, \tilde u, is complex valued, the cycle-averaged magnitude is given by: |\tilde u| = \frac{1}{2} \Re (\tilde u \cdot \tilde u^*). This allows us to solve governing equations of the form:

(4)

\begin{align}
-\omega^2 M(|\tilde u|) \tilde u + j \omega C(|\tilde u|) \tilde u + \nabla \cdot (-K(|\tilde u|) \nabla \tilde u) = \tilde F(|\tilde u|) &\text { on } \Omega \\
\mathbf{n} \cdot (K(|\tilde u|) \nabla \tilde u) + A(|\tilde u|) \tilde u = \tilde f(|\tilde u|) &\text{ on } \Gamma_1 \\
\tilde u = \tilde g(|\tilde u|) &\text{ on } \Gamma_2
\end{align}

This nonlinear partial differential equation can actually be solved with exactly the same algorithms described in our Solver series. So, in fact, there is almost nothing “extra” that we need to learn in order to address such problems, other than to enter the appropriate expression for the material nonlinearity. For an example of such a problem, see our tutorial model of self focusing in a BK-7 optical glass, where the refractive index, n=n_0+\gamma I, is directly dependent on I, the electromagnetic field intensity.

The material property does not need to be such a simple function of the cycle-averaged magnitude. The AC/DC Module, for instance, includes an effective H-B curve modeling approach that uses the nonlinear material relationships described in this previous blog post on modeling magnetic materials in the frequency domain.

So far, we have only considered a material response that is nonlinear with respect to the cycle-averaged field magnitude. It is, however, also possible for the sinusoidal excitation on the system to couple into a higher harmonic excitation of the system. To understand why this is so, keep in mind that our excitations vary in time as \sin(\omega t ) and that we have assumed that our solution will also vary sinusoidally in time and at the same frequency: u(t) = \tilde u \sin(\omega t ) . But, if there is any material response that depends directly on the instantaneous (not cycle-averaged) field magnitude, then we can use the trigonometric identity \sin^2(\omega t ) = (1- \cos ( 2 \omega t ))/2 to deduce that the response is actually of the form: u(t) = \tilde u_1 \sin(\omega t ) +\tilde u_2 \sin(2 \omega t ) + \cdots.

This type of response, referred to as frequency-doubling or higher harmonic generation, is fairly common in electromagnetics, especially optical systems. Although there may, in fact, be many higher harmonics, in practice, only one or two higher harmonics are likely of engineering interest. In such cases, we can write a new set of general governing equations:

(5)

\begin{align}
-\omega^2 M \tilde u_1 -j \omega C \tilde u_1 + \nabla \cdot (-K \nabla \tilde u_1) = \tilde F_1 -Q &\text { on } \Omega \\
\mathbf{n} \cdot (K \nabla \tilde u_1) + A \tilde u_1 = \tilde f_1 &\text{ on } \Gamma_1 \\
\tilde u = \tilde g_1 &\text{ on } \Gamma_2 \\
-4 \omega^2 M \tilde u_2 -j 2 \omega C \tilde u_2 + \nabla \cdot (-K \nabla \tilde u_2) = + Q &\text { on } \Omega \\
\mathbf{n} \cdot (K \nabla \tilde u_2) + A \tilde u_2 = 0 &\text{ on } \Gamma_1 \\
\tilde u_2 = 0 &\text{ on } \Gamma_2
\end{align}

where there is a general domain coupling term, Q, between the two sets of equations. For simplicity, the static components are omitted. Also keep in mind that all material properties can be dependent on the frequency and thus vary for the different harmonics. Now, what we have here is a set of governing partial differential equations with a nonlinear coupling term. Solving these equations simply requires looking to the same solution approaches that we introduced here.

Download this tutorial model to see an example that demonstrates the modeling of second-harmonic generation in the frequency domain.

There does come a point past which you can no longer exploit the assumptions that we’ve discussed here to simplify the problem, even if you have sinusoidal excitations. In such situations, you will want to shift focus to time-domain modeling. Although modeling in the time domain will take longer than it would in the frequency domain, you can capture the full temporal evolution of the solution and incorporate any kind of nonlinearity, even one that leads to a nonsinusoidal response.

Within the Application Gallery, you can find several general examples of time-domain modeling. They include:

If you have any questions about the modeling approaches discussed here and think they would be useful for your multiphysics modeling needs, please don’t hesitate to contact us.

#### Categories ##### Syed Tauseef Hossain
March 7, 2019

Very interesting. If a system is to be driven at different amplitudes, is it necessary to solve the problem in time domain considering only geometric nonlinearity? For example, in many resonators, a softening behaviour is observed at higher amplitudes due to geometric nonlinearity. Using the first approach stationary+frequency domain perturbation, presumably would be quicker than time domain solution. March 11, 2019

Hello Syed,