Model and Simulations of Contact Between a Vibrating Beam and an Obstacle Using the Damped Normal Compliance Condition

1. Introduction

The interaction of vibrating structures via contact with other system parts and components is an important topic in mechanical and structural engineering, with applications ranging from robotics, transportation, automotive, to aerospace systems and MEMS. When a beam repeatedly strikes an obstacle, its dynamic response is governed in part by the condition used to describe the contact process. Classical contact or impact models, such as the Hertz condition, often assume rigid, instantaneous contact with a restitution coefficient, but such models fail to capture the finite stiffness and energy dissipation observed in real materials; see e.g., [1,2], for the mathematical aspects. For this reason, the normal compliance condition was introduced in [3,4,5]. The normal compliance condition introduces a finite stiffness that allows small penetrations; however, without damping, they still underestimate energy losses during collisions [6]. To address these shortcomings, the damped normal compliance (DNC) law has recently been proposed that augments the normal compliance contact law with a damping term [7]. This modification yields a more realistic contact boundary condition because it permits controlled penetration while explicitly accounting for energy loss or dissipation during the process and thus is likely to improve the agreement with experiments.

The present study uses the DNC condition to model the contact between a vibrating Bernoulli beam and a reactive obstacle located below its free end. Such settings can be found in many MEMS systems, in actuators and grippers. We construct a new partial differential equation model that incorporates beam bending, and the contact boundary conditions of elastic compliance of the obstacle and damping. The main interest here is to study the effects of the beam’s flexural stiffness, the obstacle’s stiffness and its damping on the vibrations amplitude and frequencies of the beam. This may shed more light into the various noise characteristics of such systems. We note that the obstacle problem for the Gao beam has been stydied in [8] where the normal compliance contact condition was used.

The rest of the article is structured as follows. Section 2 describes in detail the setting and constructs the classical mathematical model, with the usual 4th order beam equation. In particular, it concentrates on the new DNC contact condition at the right end of the beam where the obstacle is situated. Energy considerations can be found in Section 3, where the energy balance is derived. Section 4 studies the steady states of the system. The condition on the force f, the position of the obstacle $y_{-}$, and the beams stiffness ${\alpha }^2$, for existence of contact is obtained. Section 5 describes shortly, in Theorem 5.1 the existence of weak solutions to the model, presented in [9].

To explore the system dynamics numerically, we construct an implicit finite difference algorithm in Section 6. Numerical simulations reported in Section 7 illustrate how different combinations of obstacle stiffness and damping affect beam penetration, vibrations, and return to equilibrium. For very stiff obstacles with small damping, the behavior approaches rigid Signorini contact with noticeable rebounds, whereas moderate stiffness and larger damping lead to deeper penetration and faster decay of oscillations. These computational experiments underscore the wide range of behaviors captured by the DNC law and highlight its ability to bridge between purely elastic and highly dissipative contacts [4,7]. The paper ends with Section 8, which presents a short summary, some conclusions, and some unresolved problems, which we consider worth further study.

2. The Setting and the Model

This section presents the setting and constructs a classical model for the process of contact between a vibrating Bernoulli beam and a reactive object or obstacle, which is situated below the right end of the beam. The setting is illustrated in Figure 1.

The system consists of a Bernoulli beam that in its reference configuration occupies the (scaled with L – the original length of the beam) interval $[0,1]$. The beam’s displacement from this reference configuration is denoted by $u(x,t)$, scaled by L, chosen as positive when above the x-axis. The beam is acted upon by a distributed force $f=f(x,t)$. The right end of the beam, at $x=1$, can come in contact with a reactive obstacle, marked in Figure 1 as `obstacle.’ The position of the obstacle is at $y=y_{-}$, and it is where contact can take place. Our main interest lies in what happens during contact between the beam’s end and the obstacle. We denote partial derivatives with subscripts, so $u_x=\partial u/\partial x$, etc. The force applied to the beam is such that the beam can press against the obstacle, possibly leading to repeated contact.

We turn to the classical formulation of the model for the process.

Model 1.

Given the force $f=f(x,t)$, the excitation $\phi =\phi \left(t\right)$, and the coefficients $\alpha ,\kappa ,\beta$, all positive constants, find the displacement field $u:[0,1]\times [0,T]\to \mathbb{R}$, for $0<T$, such that:

$$\begin{array}{cc}\hfill & u_{tt}+{\alpha }^2{u}_{xxxx}=f\mathrm{in}{\Omega }_T=\{0<x<1,0<t<T\}\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill & u(0,t)=u_x(0,t)=0\mathrm{for}0\le t\le T,\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill & u_{xx}(1,t)=0,\mathrm{for}0\le t\le T,\hfill \end{array}$$

(1c)

$$\begin{array}{cc}\hfill & -{\alpha }^2{u}_{xxx}(1,t)=\left\{\begin{array}{cc}0,\hfill & u(1,t)>y_{-},\hfill \\ \kappa {(u(1,t)-y_{-})}_{-}-\beta u_t(1,t),\hfill & u(1,t)\le y_{-}.\hfill \end{array}\right.,\hfill \end{array}$$

(1d)

$$\begin{array}{cc}\hfill & u(x,0)=u_0\left(x\right),u_t(x,0)=v_0\left(x\right),0\le x\le 1.\hfill \end{array}$$

(1e)

The initial displacement $u(x,0)=u_0\left(x\right)$ and the initial velocity $v(x,0)=v_0\left(x\right)$ are given functions satisfying the compatibility conditions $u_0\left(0\right)=0$, $v_0\left(0\right)=0$, $u_{0xx}\left(1\right)=0$ and $u_0\left(1\right)>y_{-}$ and then $u_{0xxx}\left(1\right)=0$.   

Here, u and x are scaled with L, and so is $\phi$; ${\alpha }^2=E_{Y}I/\rho A{L}^4$, $\rho$ is the density of the beam material, A is the cross-sectional area, $E_Y$ is Young’s modulus (elasticity modulus), and I is the second moment of the area of the beam’s cross-section. The external force f is scaled with $LA$, so it is force per unit volume of the beam. In addition, ${(·)}_{-}$ denotes the negative part function, that is ${\left(r\right)}_{-}=\left|r\right|$ when $r<0$ and ${\left(r\right)}_{-}=0$ when $r>0$.

Remark 1.

A note on the DNC contact condition, (1d), [4]. When there is no contact, $u(1,t)>y_{-}$, the traction or stress at $x=1$ is zero since the end is free. When there is contact, $u(1,t)\le y_{-}$, the obstacle reaction is given by

$$R=\kappa {(u(1,t)-y_{-})}_{-}-\beta u_t(1,t),$$

where κ is the obstacle stiffness and β is its damping coefficient. Thus, the DNC condition explicitly takes into account the energy loss during contact. It is also noted that R may be positive or negative, although the stiffness term always positive, as it points up. However, since when the obstacle is soft (small κ), the damping is large (large β) and the velocity is large and in the upward direction, R may be negative. As we show below, the energy dissipation is always negative when the end is moving. Furthermore, $|u(1,t)-y_{-}|$ is the depth of the penetration of the beam into the obstacle, when in contact.

For the sake of simplicity, we introduce the cutoff function

$${\chi }_{*}\left(u\right)=\left\{\begin{array}{cc}0\hfill & u(1,t)>y_{-}\hfill \\ 1\hfill & u(1,t)\le y_{-}.\hfill \end{array}\right.$$

(2)

It vanishes when there is no contact and, therefore the obstacle resistance can be written concisely as

$$R=\left(\kappa (y_{-}-u(1,t))-\beta u_t(1,t)\right){\chi }_{*}\left(u\right).$$

(3)

3. Energy Balance and Dissipation

We derive the energy balance in the system, assuming that the solutions are sufficiently smooth for the manipulations to make sense. We multiply (1a) with $u_t$ and integrate over $[0,1]\times [0,t]$, for $0<t\le T$. We deal with the different terms separately. The first term, after integration by parts, yields

$${\displaystyle \frac{1}{2}}{\int }_0^1{u}_t^2{|}_0^{t}dx={\displaystyle \frac{1}{2}}{\int }_0^1{u}_t^2(x,t)dx-{\displaystyle \frac{1}{2}}{\int }_0^1{v}_0^2\left(x\right)dx=E_k\left(t\right)-E_k\left(0\right).$$

Here, the kinetic energy, at time t, is given by

$$E_k\left(t\right)={\displaystyle \frac{1}{2}}{\int }_0^1{u}_t^2(x,t)dx.$$

(4)

Next, we note that

$${\alpha }^2{\int }_0^t{\int }_0^1{u}_{xxxx}{u}_{t}dxd\tau ={\alpha }^2{\int }_0^t{u}_{xxx}{u}_t{|}_0^{1}d\tau -{\alpha }^2{\int }_0^t{\int }_0^1{u}_{xxx}{u}_{xt}dxd\tau .$$

We have, using (3),

$${\alpha }^2{\int }_0^t{u}_{xxx}{u}_t{|}_0^{1}d\tau =-{\int }_0^t(\kappa (y_{-}-u(1,\tau ))u_t(1,\tau )-\beta u_t^2(1,\tau )){\chi }_{*}\left(u\right)d\tau$$

$$-{\alpha }^2{\int }_0^t{u}_{xxx}(0,\tau ){\phi }^{\prime }\left(\tau \right)d\tau .$$

Moreover, using the initial condition $u_0\left(1\right)>y_{-}$,

$${\int }_0^t\left(\kappa (u(1,\tau )-y_{-})u_t(1,\tau )\right){\chi }_{*}\left(u\right)d\tau ={\displaystyle \frac{\kappa }{2}}{(u(1,t)-y_{-})}^2{\chi }_{*}\left(u\right).$$

Furthermore, using integration by parts and the boundary conditions and the facts that $u_x(0,t)=0$ and $u_{xx}(1,t)=0$, we obtain

$${\alpha }^2{\int }_0^t{\int }_0^1{u}_{xxx}{u}_{xt}dxd\tau ={\alpha }^2{\int }_0^t{u}_{xx}{u}_{xt}{{|}_0^{1}d\tau -{\displaystyle \frac{{\alpha }^2}{2}}{\int }_0^1{u}_{xx}^2|}_0^{t}dx$$

$$=-{\displaystyle \frac{{\alpha }^2}{2}}{\int }_0^1{u}_{xx}^2(x,t)dx+{\displaystyle \frac{{\alpha }^2}{2}}{\int }_0^1{u}_{0xx}^2\left(x\right)dx.$$

We now define the potential energy, at time t, of the system as

$$E_p\left(t\right)={\displaystyle \frac{{\alpha }^2}{2}}{\int }_0^1{u}_{xx}^2(x,t)dx+{\displaystyle \frac{\kappa }{2}}{(u(1,t)-y_{-})}^2{\chi }_{*}\left(u\right).$$

(5)

Then, using (4) and (5) and the above expressions, we obtain the following.

Result 1.

The system energy balance,

$$E\left(t\right)=E_k\left(t\right)+E_p\left(t\right)$$

$$={\displaystyle \frac{1}{2}}{\int }_0^1{u}_t^2(x,t)dx+{\displaystyle \frac{{\alpha }^2}{2}}{\int }_0^1{u}_{xx}^2(x,t)dx+{\displaystyle \frac{\kappa }{2}}{(u(1,t)-y_{-})}^2{\chi }_{*}\left(u\right)$$

$$=E_k\left(0\right)+E_p\left(0\right)+{\int }_0^t{\int }_0^{1}f(x,t)u_t(x,t)dxd\tau -\beta {\int }_0^t{u}_t^2(1,\tau ){\chi }_{*}\left(u\right)d\tau .$$

(6)

The first term on the right describes the initial kinetic energy and the second term the initial elastic potential energy. The third term describes the energy supplied or removed by the applied force, the last term is the energy dissipated by the obstacle. It follows, since the last term is negative when the end is in motion and in contact, that the obstacle’s resistance dissipates energy from the system.   

4. Steady States

This section studies the steady states of the system. We assume that the applied force f is a negative constant, since if $f>0$ then the end does not come into contact with the obstacle. The steady solutions, denoted by $\overline{u}\left(x\right)$, are found by setting the time derivatives in (1b)–(1d) to zero. Thus, we obtain the following ODEs model:

Model 2.

Given the force $f=const.$, $\alpha ,\kappa ,$ positive constants, find the displacement field $\overline{u}:[0,1]\to \mathbb{R}$, such that:

$$\begin{array}{cc}\hfill & {\alpha }^2{\overline{u}}_{xxxx}=f,0<x<1,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & \overline{u}\left(0\right)={\overline{u}}_x\left(0\right)=0,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & {\overline{u}}_{xx}\left(1\right)=0,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & -{\alpha }^2{\overline{u}}_{xxx}\left(1\right)=\kappa (\overline{u}\left(1\right)-y_{-}){\chi }_{*}\left(\overline{u}\left(1\right)\right).\hfill \end{array}$$

(10)

Integrating four times with respect to x gives the general polynomial solution for (7) and applying the boundary conditions in (8), we obtain

$$\overline{u}\left(x\right)={\displaystyle \frac{1}{2}}B{x}^2+{\displaystyle \frac{1}{6}}A{x}^3+{\displaystyle \frac{f}{24{\alpha }^2}}{x}^4.$$

Applying the boundary conditions at $x=1$, we have

$${\overline{u}}^{″}\left(1\right)=B+A+{\displaystyle \frac{f}{2{\alpha }^2}}=0,$$

hence

$$\overline{u}\left(x\right)=-\left({\displaystyle \frac{f}{4{\alpha }^2}}+{\displaystyle \frac{A}{2}}\right)x^2+{\displaystyle \frac{A}{6}}{x}^3+{\displaystyle \frac{f}{24{\alpha }^2}}{x}^4.$$

(11)

and

$$-{\alpha }^2{\overline{u}}^{‴}\left(1\right)=-{\alpha }^{2}A-f=\kappa \left(\overline{u}\left(1\right)-y_{-}\right){\chi }_{*}\left(\overline{u}\left(1\right)\right).$$

(12)

This implies that there are two cases to consider:

(i) No contact: When $\overline{u}\left(1\right)-y_{-}>0$ then ${\chi }_{*}\left(\overline{u}\left(1\right)\right)=0$, and we find $A=-{\displaystyle \frac{f}{{\alpha }^2}}$. Therefore, in this case, the solution is

$$\overline{u}\left(x\right)={\displaystyle \frac{f}{24{\alpha }^2}}\left(x^2-4x+6\right)x^2.$$

(13)

Evaluating this at $x=1$, we find $\overline{u}\left(1\right)={\displaystyle \frac{f}{8{\alpha }^2}}$, therefore:

Result 2.

[No contact] The condition for no contact is

$${\displaystyle \frac{f}{8{\alpha }^2}}>y_{-}.$$

(14)

Then ${\chi }_{*}\left(\overline{u}\left(1\right)\right)=0$ and the solution is just the usual one of a linear beam with a free right end, and is given in (13).   

Physically, this means that the (scaled) force is insufficient to bend the beam enough to contact the obstacle, since the beam stiffness opposes it.

(ii) Contact: When $\overline{u}\left(1\right)-y_{-}<0$, then $\left(\overline{u}\left(1\right)-y_{-}\right){\chi }_{*}\left(\overline{u}\left(1\right)\right)=y_{-}-\overline{u}\left(1\right)$, ${\chi }_{*}\left(\overline{u}\left(1\right)\right)=1$, and we find from (11) and (12) that

$$-{\displaystyle \frac{{\alpha }^2}{\kappa }}A-{\displaystyle \frac{f}{\kappa }}=y_{-}-\overline{u}\left(1\right)=y_{-}-\left(-{\displaystyle \frac{5f}{4{\alpha }^2}}-{\displaystyle \frac{A}{3}}\right),$$

$$=y_{-}+{\displaystyle \frac{5f}{24{\alpha }^2}}+{\displaystyle \frac{A}{3}}.$$

Then, after rearranging and some manipulations, we obtain

$$A={\displaystyle \frac{-3\kappa }{(\kappa +3{\alpha }^2)}}\left(y_{-}+f\left({\displaystyle \frac{5}{24{\alpha }^2}}+{\displaystyle \frac{1}{\kappa }}\right)\right),$$

(15)

which depends on ${\alpha }^2,f$ and $\kappa$. We note that $f<0,y_{-}<0$, hence $A>0$. Thus, with A given in (15), the solution of the steady problem when contact takes place is

$$\overline{u}\left(x\right)=-\left({\displaystyle \frac{f}{4{\alpha }^2}}+{\displaystyle \frac{A}{2}}\right)x^2+{\displaystyle \frac{A}{6}}{x}^3+{\displaystyle \frac{f}{24{\alpha }^2}}{x}^4.$$

(16)

In particular,

$$\overline{u}\left(1\right)={\displaystyle \frac{1}{(\kappa +3{\alpha }^2)}}\left({\displaystyle \frac{3}{8}}f+\kappa y_{-}\right).$$

(17)

We conclude that in the case with contact, we have:

Result 3.

[Contact] Assume that

$${\displaystyle \frac{f}{8{\alpha }^2}}\le y_{-}.$$

(18)

Then, there is contact between the beam and the obstacle and the steady solution $\overline{u}$ is given in (16).   

Physically, this means that the (scaled) force to the beam’s stiffness ratio is sufficient to bend the beam so it comes in contact with the obstacle.

Figure 2 shows the computed solutions of the steady states, using formulas (13) or (16), depending on the inequality (14). In the simulations ${\alpha }^2=1,\kappa =2$ and $y_{-}=-0.2$ and a constant force f are used. The figure shows the steady bending of the beam for various values of the force, and contact is established when the force is close to $f=-1.60$.

It is found that under a constant distributed load $f<0$ and the DNC contact condition with stiffness $\kappa >0$ and obstacle position $y_{-}$, the beam deflection naturally divides into two cases characterized by

$$\overline{r}={\displaystyle \frac{f}{8{\alpha }^2{y}_{-}}}.$$

When $\overline{r}<1$ there is no-contact, and the deflection of the free end is above the obstacle. When $\overline{r}\ge 1$ there is contact. Clearly, $\overline{r}=1$ is when the end of the beam just touches the obstacle, without penetrating it.

In the example, Figure 2, the transition case $\overline{r}=1$ is when $f=-1.60$, since ${\alpha }^2=1$ and $y_{-}=-0.2$. This can be clearly seen in the figure. The positions of the beam’s end is given by (17) for the different values of f.

5. Abstract Formulation and Existence

The detailed proof of the existence of a weak solution to the problem can be found in [9]. Here, we just provide a short description and state the main existence theorem. The proof is based on a sequence of regularized problems, and then passing to the limits. The main issue is that the boundary conditions for $v=u_t$ are not well defined when we only know that $u_t\left(t,·\right)$ has values in $L^2$. Also, the discontinuity of ${\chi }_{*}$ causes some difficulty. For mathematical reasons, we add viscosity to the beam equation in the form of $\epsilon u_{xxxxt}$, for small $\epsilon$, which allows us to obtain better estimates.

We first consider a problem in which ${\chi }_{\delta }$ is regularized by replacing it with a piecewise linear Lipschitz approximation given by

$${\chi }_{\eta }\left(r\right)=\left\{\begin{array}{cc}0\hfill & r>y_{-}\hfill \\ r\hfill & \left(y_{-}-\eta \right)\le r\le y_{-}\hfill \\ 1\hfill & r\le \left(y_{-}-\eta \right),\hfill \end{array}\right.$$

for a small $0<\eta$. Eventually, we pass to the limits $\epsilon ,\eta \to 0$.

We set the problem in terms of $v=u_t$, and then

$$u\left(x,t\right)=u_0\left(x\right)+{\int }_0^{t}v\left(x,s\right)ds.$$

We use the notation ${\Omega }_T=\{0<x<1,0<t<T\}$,  and let $H=L^2(0,1)$ and $V=\{w\in H^2(0,1);w\left(0\right)=w_x\left(0\right)=0\}$, with the usual inner products and norms.

The classical formulation of the regularized problem is:

Find $v:[0,T]\to V$ such that

$$v_t+\epsilon v_{xxxx}+{\alpha }^2{u}_{xxxx}=f\mathrm{in}{\Omega }_T,$$

(19)

$$u(0,t)=0,u_x(0,t)=0\mathrm{for}0\le t\le T,$$

(20)

$$u_{xx}(1,t)=0,\mathrm{for}0\le t\le T,$$

(21)

$$-{\alpha }^2{u}_{xxx}(1,t)={\chi }_{\eta }\left(u\left(1,t\right)\right)\left(\kappa (y_{-}-u(1,t))-\beta v(1,t)\right),$$

(22)

$$v(x,0)=v_0\left(x\right),v_0\in L^2\left(0,1\right),$$

(23)

$$u\left(x,t\right)=u_0\left(x\right)+{\int }_0^{t}v\left(x,s\right)ds.$$

(24)

Next, we set the problem in a variational form. The function space in which one expects to seek weak solutions is the Sobolev space.

$$V=\left\{u(·,t)\in H^2(0,1)|u(0,t)=0,u_x(0,t)=0,\right\}$$

and

$$\mathcal{V}=L^2(0,T,V).$$

For the time dependence on $t\in [0,T]$, we require: $u\in \mathcal{V}$ and $v=u_t\in L^2(0,T)\times H^2(0,1)$.

Denote such a space by $\mathcal{W}$,

$$\mathcal{W}=\left\{u:[0,T]\to \mathcal{V}|u_t\in L^2(0,T)\times \left(H^2(0,1)\right)\right\}.$$

We choose $\varphi (x,t)\in L^2(0,T;V)$ as a test function. Multiplying Equation (19) by $\varphi$ and integrating on $[0,1]$, using various manipulations, Green’s identities, and the boundary conditions, we obtain the following.

$$\begin{array}{cc}\hfill & {\int }_0^T{\int }_0^1{u}_{tt}\varphi dxdt+{\alpha }^2{\int }_0^T{\int }_0^1{u}_{xx}{\varphi }_{xx}dxdt+\epsilon {\int }_0^T{\int }_0^1{u}_{xxt}{\varphi }_{xx}dxdt\hfill \\ \hfill & -{\int }_0^T\left[\kappa (u(1,t)-y_{-})-\beta u_t(1,t)\right]{\chi }_{\eta }\left(u(1,t)\right)\varphi (1,t)dt={\int }_0^T{\int }_0^{1}f\varphi dxdt.\hfill \end{array}$$

(25)

The weak or variational formulation of the problem with viscosity, Model (19)–(24), is the following:

Model 3.

Given $f\in L^2\left(0,T)\times L^2(0,1\right)$, and positive constants $\alpha ,\kappa ,\beta ,\epsilon ,\eta$, find $u\in \mathcal{V}$, such that $u_t\in \mathcal{V}$, and for every test function $\varphi \in \mathcal{V}$ that satisfies $\varphi (x,T)=0$, (25) holds, together with

$$u(0,t)=0,u_x(0,t)=0,0\le t\le T,$$

and

$$u(x,0)=u_0\left(x\right),u_t(x,0)=v_0\left(x\right),0\le x\le 1.$$

As was noted in [9], we need to control the boundary term and hence, we introduce the truncation function.

$$m\left(r\right)=\left\{\begin{array}{c}0,0\le r,\hfill \\ r,-M\le r\le 0,\hfill \\ -M,r\le -M.\hfill \end{array}\right.$$

Then, we replace the term $\kappa (u(1,t)-y_{-})$ with

$$m\left(\kappa (u(1,t)-y_{-})\right),$$

which is bounded and Lipschitz continuous.

We have the following existence theorem, established in [9].

Theorem 1

([9]). There exists a variational solution u to Model 3 in which the term $\kappa (u(1,t)-y_{-})$ is truncated as $m\left(\kappa (u(1,t)-y_{-})\right)$, and ${\chi }_{\eta }$ is replaced with ${\chi }_{*}$. The solution with $\epsilon =0$ is weaker.

The proof uses a fixed-point argument and some regularization. Then, sufficient estimates are obtained that allow us to consider the case without viscosity, that is, the case $\epsilon \to 0$.

6. Finite Difference Algorithms

We use an implicit second-order finite-difference algorithm for the simulations of the problem. The beam is divided into M uniform length segments $\Delta x=L/M$. The simulations are performed over the time interval $[0,T]$, and the time step is given by $\Delta t=T/N$. Let $u_i^j$ be the approximation of the solution at the node $x_i=i\Delta x$ at time $t_j=j\Delta t$, that is, $u_i^j\approx u(x_i,t_j)$, for $i=0,1,\dots ,M$ and $j=0,1,\dots ,N$. First, we note that, thanks to the initial conditions in Equation (1e), we have

$$u_i^0=u_0\left(x_i\right),u_i^1=\Delta t{v}_0\left(x_i\right)+u_0\left(x_i\right),0\le i\le M.$$

(26)

Moreover, thanks to the Dirichlet boundary condition in Equation (1b), we have

$$u_0^j=0,\forall 1\le j\le N.$$

(27)

For $i=1,2,\dots ,M$ and $j=1,2,\dots ,N-1$, the implicit finite difference scheme is given by

$${\displaystyle \frac{u_i^{j+1}-2u_i^j+u_i^{j-1}}{\Delta t^2}}+{\alpha }^2{\displaystyle \frac{u_{i+2}^{j+1}-4u_{i+1}^{j+1}+6u_i^{j+1}-4u_{i-1}^{j+1}+u_{i-2}^{j+1}}{\Delta x^4}}=f_i^{j+1}.$$

(28)

Rearranging terms, and defining $\gamma ={\alpha }^2{\displaystyle \frac{\Delta t^2}{\Delta x^4}}$, we have

$$\gamma u_{i-2}^{j+1}-4\gamma u_{i-1}^{j+1}+(1+6\gamma )u_i^{j+1}-4\gamma u_{i+1}^{j+1}+\gamma u_{i+2}^{j+1}=\Delta t^2{f}_i^{j+1}+2u_i^j-u_i^{j-1}.$$

(29)

Equation (29) is used to construct a linear system at every time level which is solved to obtain the discrete solution at that level. We note that when $i=1$, we have the ghost node$u_{-1}^{j+1}$, for which we use the boundary condition in Equation (1b). In particular, using a second-order approximation of the derivative, this boundary condition gives us

$$u_x(0,t_{j+1})\approx {\displaystyle \frac{u_1^{j+1}-u_{-1}^{j+1}}{2\Delta x}}=0\Rightarrow u_{-1}^{j+1}=u_1^{j+1}.$$

Therefore, the finite-difference scheme Equation (29) for $i=1$ becomes

$$(1+7\gamma )u_1^{j+1}-4\gamma u_2^{j+1}+\gamma u_3^{j+1}=\Delta t^2{f}_1^{j+1}+2u_1^j-u_1^{j-1}+4\gamma \phi \left(t_{j+1}\right).$$

(30)

We also note that when $i=M-1$ and when $i=M$, there are two ghost nodes: $u_{M+1}^{j+1}$ and $u_{M+2}^{j+1}$. Using the boundary condition in Equation (1c), we have

$$u_{xx}(1,t_{j+1})\approx {\displaystyle \frac{u_{M+1}^{j+1}-2u_M^{j+1}+u_{M-1}^{j+1}}{\Delta x^2}}=0\Rightarrow u_{M+1}^{j+1}=2u_M^{j+1}-u_{M-1}^{j+1}.$$

Thus, the finite difference scheme for $i=M-1$ becomes

$$\begin{array}{cc}\hfill \gamma u_{M-3}^{j+1}-4\gamma u_{M-2}^{j+1}& +(1+5\gamma )u_{M-1}^{j+1}\hfill \\ \hfill & -2\gamma u_M^{j+1}=\Delta t^2{f}_{M-1}^{j+1}+2u_{M-1}^j-u_{M-1}^{j-1}.\hfill \end{array}$$

(31)

The contact boundary condition Equation (1d) is taken into account depending on whether the beam had contact with the obstacle at the previous time or not. We have the following two cases:

• No contact: $u_M^j>y_{-}$

  $$u_{xxx}(1,t_{j+1})\approx {\displaystyle \frac{-\frac{1}{2}{u}_{M-2}^{j+1}+u_{M-1}^{j+1}-u_{M+1}^{j+1}+\frac{1}{2}{u}_{M+2}^{j+1}}{\Delta x^3}}=0$$

  $$\Rightarrow u_{M+2}^{j+1}=u_{M-2}^{j+1}-2u_{M-1}^{j+1}+2u_{M+1}^{j+1}=u_{M-2}^{j+1}-4u_{M-1}^{j+1}+4u_M^{j+1}.$$
  Thus, in the case of no contact, the finite difference scheme is

  $$2\gamma u_{M-2}^{j+1}-4\gamma u_{M-1}^{j+1}+(1+2\gamma )u_M^{j+1}=\Delta t^2{f}_M^{j+1}+2u_M^j-u_M^{j-1}.$$

  (32)
• Contact: $u_M^j\le y_{-}$

  $$\begin{array}{cc}\hfill -{\alpha }^2{u}_{xxx}(1,t_{j+1})& \approx -{\alpha }^2{\displaystyle \frac{-\frac{1}{2}{u}_{M-2}^{j+1}+u_{M-1}^{j+1}-u_{M+1}^{j+1}+\frac{1}{2}{u}_{M+2}^{j+1}}{\Delta x^3}}\hfill \\ \hfill & =\kappa (y_{-}-u_M^{j+1})-\beta {\displaystyle \frac{u_M^{j+1}-u_M^j}{\Delta t}}\hfill \\ \hfill \Rightarrow u_{M+2}^{j+1}& =-{\displaystyle \frac{2\kappa }{{\alpha }^2}}\Delta x^3(y_{-}-u_M^{j+1})+{\displaystyle \frac{2\beta }{{\alpha }^2}}\Delta x^3{\displaystyle \frac{u_M^{j+1}-u_M^j}{\Delta t}}\hfill \\ \hfill & +u_{M-2}^{j+1}-4u_{M-1}^{j+1}+4u_M^{j+1}\hfill \end{array}$$
  Thus, in the case of no contact, the finite difference scheme is

  $$\begin{array}{c}\hfill 2\gamma u_{M-2}^{j+1}-4\gamma u_{M-1}^{j+1}+\left(1+2\gamma \left({\displaystyle \frac{\Delta x^3}{{\alpha }^2}}\left(\kappa +{\displaystyle \frac{\beta }{\Delta t}}\right)+1\right)\right)u_M^{j+1}\\ \hfill =\Delta t^2{f}_M^{j+1}+2\left(1+\gamma {\displaystyle \frac{\beta \Delta x^3}{{\alpha }^2\Delta t}}\right)u_M^j-u_M^{j-1}+2\gamma {\displaystyle \frac{\Delta x^3}{{\alpha }^2}}\kappa y_{-}.\end{array}$$

  (33)

The algorithm is as follows.

Algorithm 1 Implicit finite difference scheme for the vibrating beam

Initialization: 1: for $i=0$ to M do 2:     Set $u_i^0$ and $u_i^1$ using Equation (26). 3: end for 3: Time Marching: 4: for $j=1$ to $N-1$ do    (1) Set Dirichlet boundary condition on the left using Equation (27).    (2) Construct the $M\times M$ matrix A and $M\times 1$ vector b. 5:     Assign the first row of A and first entry of b according to Equation (30). 6:     for $i=2$ to $M-2$ do 7:         Assign i-th row of A and i-th entry of b according to Equation (29) 8:     end for 9:     Assign $M-1$-th row of A and $M-1$-th entry of b according to Equation (31) 10:     Assign M-th row of A and M-th entry of b. 11:     if $u_M^j>y_{-}$ then                   ▹ No-contact scenario 12:         Use Equation (32) 13:     else                        ▹ Contact scenario 14:         Use Equation (33) 15:     end if    (3) Solve the linear system with matrix A and right hand side b 16: end for

7. Simulations

The implicit finite difference Algorithm 1 was implemented in MATLAB, and the backslash function was used to solve the linear system. In all simulations $v_0\left(x\right)\equiv 0$, the beam length was $L=1$, the final time was $T=5$, and $\Delta x=\Delta t=0.01$. Moreover, the force was taken as $f=-0.25$ and ${\alpha }^2=1.296$. The obstacle was located at $y_{-}=-0.02$. In this section, we report some of the more interesting computer experiments. Typically, each simulation took about 10-15 seconds, running in MATLAB.

The following are example simulations for various obstacle stiffnesses, $\kappa$, and damping coefficients $\beta$. Typically, in Figure 3 we show the position of the right end of the beam over time for different values of $\kappa$ and $\beta$. It is seen that, progressively, the penetration decreases, and substantially changes between the values $\kappa =10$ and $\kappa =1000$, in the latter case there is virtually no penetration.

Result 4.

We observe the following:

• For large $\kappa$ ($>10$) and small $\beta$ ($\approx 0.1)$, the contact process behaves similarly to that with a Signorini (rigid) contact condition: very little penetration, high contact force and repeated bounces when the initial energy is sufficient.
• For moderate stiffness $\kappa$ ($\approx 1$), one obtains a normal compliance-like penetration, with partial rebounds damped out, depending on the values of $\beta$.
• For small $\kappa$ ($\approx 0.1$) and large $\beta$ ($\approx 10)$, the obstacle effectively behaves like a “soft, dissipative" material, such as tissue, allowing for deeper penetration and higher energy dissipation.

These results confirm that the DNC law can model a whole spectrum of contact behaviors, bridging between purely elastic and purely dissipative extremes.

8. Conclusions and Unresolved Questions

This work studies the vibrations of a Bernoulli beam that can come into contact with an obstacle, which is situated under its right end. The contact process is described by the Damped Normal Compliance boundary condition. The model consists of the usual 4th-order dynamic beam equation, together with the initial and boundary conditions. The interest lies in the obstacle’s resistance R, (3), which includes stiffness and damping terms, and is active once the right end is in contact.

The system energy balance is derived in Section 3, showing that the system loses energy because of damping, once the beam is in contact. The steady states of the system are found in Section 4. It is found that the ratio $\overline{r}$, where

$$\overline{r}={\displaystyle \frac{f}{8{\alpha }^2{y}_{-}}}$$

controls when there is contact. Indeed, when $\overline{r}<1$ there is no contact in the steady state and when $\overline{r}\ge 1$ the beam is in contact. When $\overline{r}=1$ the beam’s end just touches the obstacle. The variational formulation of the problem is provided in Section 5, that includes the functional setting. It is noted that the term $v=u_t(1,t)$ needs a special treatment, since for $v=u_t\in V$ the trace $v(1,t)$ is not well defined. We circumvent this problem by considering regularized problems with viscosity, the solutions of which are sufficiently regular to allow for the trace to exist. The existence of the variational solutions to the regularized problems are stated in Theorem 5.1. Moreover, the solutions of the problem without viscosity are obtained as limits and can be found in the weaker space $\mathcal{W}$. The proofs can be found in [9].

A finite difference scheme and algorithm for the simulations of the problem can be found in Section 6. The algorithm is implicit, and special attention is paid to the discretized contact conditions. Section 7 provides the simulation results. It is found that for a sufficiently large force, the beam penetrates the obstacle. Then, depending on the stiffness and damping, the simulations indicate that when the stiffness decreases and the damping increases, the penetration increases and the dying out of the oscillations is faster. These are summarized in the Results 4.

These results confirm that the DNC law can model a whole spectrum of contact behaviors, bridging between purely elastic and purely dissipative extremes.

We turn to some unresolved problems. More simulations and a numerical study of the dependence of the frequencies of vibration on the beam rigidity and the obstacle stiffness may be of practical interest. The case when there are two obstacles, one above and the other below, is of interest, especially numerically, since the vibration frequencies of the beam depend on the beam’s rigidity and the stiffness of the obstacles.

Further analysis of the problem is of interest.

Another aspect that may be of interest is the possible addition of friction to the contact, so that when the beam’s end, when in contact, undergoes friction and frictional energy dissipation, [6,10,11,12].

8.1. Steady States – Numerical Solution

In addition to the explicit formulas for the steady states of the beam, we implemented a numerical method, based on the dynamic case above, and compared the numerical and the closed-form solutions, to estimate the numerical errors.

Steady states do not involve time, so we discretize the beam into M equal segments of length $\Delta x={\displaystyle \frac{1}{M}}$ such that $x_i=i\Delta x,i=0,1,...,M$, and let ${\overline{u}}_i\approx \overline{u}\left(x_i\right)$ denote the displacements at $x_i$.

For interior points $i=2,3,...,M-2$, we can approximate the fourth derivative by the standard central difference scheme. Thus the discretized equation for ${\overline{u}}_{xxxx}\left(x_i\right)=f$ becomes

$${\alpha }^2{\displaystyle \frac{u_{i+2}-4u_{i+1}+6u_i-4u_{i-1}+u_{i-2}}{{(\Delta x)}^4}}=f$$

(34)

At the left end $(x=0),i=0$, ${\overline{u}}_0=0$ and ${\overline{u}}_x=0$, so we use the second-order approximation of ${\overline{u}}_x\left(0\right)\approx 0$ to be

$${\overline{u}}_x\left(0\right)\approx {\displaystyle \frac{-3{\overline{u}}_0+4{\overline{u}}_1-{\overline{u}}_2}{2\Delta x}}=0⟹-3{\overline{u}}_0+4{\overline{u}}_1-{\overline{u}}_2=0.$$

But ${\overline{u}}_0=0$, which implies that ${\overline{u}}_2=4{\overline{u}}_1$.

Now, at the right end $(x=1),i=M$, we have

$$-{\alpha }^2{\overline{u}}_{xxx}\left(1\right)=\kappa {\left(\overline{u}\left(1\right)-y_{-}\right)}_{-}.$$

We approximate the third derivative with the second-order backwards difference as

$$G\equiv -{\alpha }^2{\displaystyle \frac{{\overline{u}}_M-3{\overline{u}}_{M-1}+3{\overline{u}}_{M-2}-{\overline{u}}_{M-3}}{{(\Delta x)}^3}}=\kappa {\left({\overline{u}}_M-y_{-}\right)}_{-}.$$

(35)

Then, we apply the DNC condition:

$${\overline{u}}_M>y_{-}\Rightarrow G=\kappa {\left({\overline{u}}_M-y_{-}\right)}_{-}=0;$$

(36)

$${\overline{u}}_M<y_{-}\Rightarrow G=\kappa \left(y_{-}-{\overline{u}}_M\right).$$

(37)

Figure 4. Comparing the steady state numerical solution for the beam without or with contact, for $f=-2.0$ with the exact solution. The obstacle is at $y_{-}=-0.2$.

Figure 4. Comparing the steady state numerical solution for the beam without or with contact, for $f=-2.0$ with the exact solution. The obstacle is at $y_{-}=-0.2$.

The following is the algorithm we used for the numerical solutions of the steady states. Comparison with the closed-form solution (13) or (16) provides confidence in our numerical simulations. Indeed, the comparison in Figure 3 shows numerical results that are very close to the closed-form results. Furthermore, Table 1 provides a more precise comparison, which is depicted in Figure 5.

Algorithm 2 Case-Based Solver for the Steady-State Beam with DNC Condition

Input: number of intervals M, beam parameters $(\alpha ,\kappa )$, external force f, obstacle position $y_{-}$. Output: steady-state displacement array $\overline{u}$ at grid points $x_i,i=0,\dots ,M$. 1: Discretization: Compute $\Delta x\leftarrow 1/M$ and set $x_i\leftarrow i\Delta x$ for $i=0,\dots ,M$. Initialize an array $\overline{u}$ (all zeros). 2: Case 1: No-Contact Solve Construct the system matrix $\mathbf{A}$ and right-hand side $\mathbf{rhs}$ by imposing: (a) Beam equation ${\alpha }^2{\overline{u}}_{xxxx}=f$ for $i=2,\dots ,M-2$, using Equation (34) (b) Left boundary conditions at $x=0$, (c) No-contact condition ${\overline{u}}_{xxx}\left(1\right)=0$ using the formula in Equation (35) and the condition Equation (36). Solve the linear system $\mathbf{A}\overline{u}=\mathbf{rhs}$ to obtain ${\overline{u}}_A$. Check consistency: if ${\overline{u}}_A\left(M\right)<y_{-}$, then set successA= false, else true. 3: Case 2: Contact Solve Rebuild the system matrix $\mathbf{A}$, adjusting only the contact boundary condition using the formula in Equation (35) and the condition Equation (37). Solve $\mathbf{A}\overline{u}=\mathbf{rhs}$ to obtain ${\overline{u}}_B$. Check consistency: if ${\overline{u}}_B\left(M\right)\ge y_{-}$, then set successB= false, else true. 4: Compare and Select 5: if successA is true then 6:     ${\overline{u}}_{\mathrm{sol}}\leftarrow {\overline{u}}_A$                                    ▹ No-contact scenario is consistent 7: else 8:     ${\overline{u}}_{\mathrm{sol}}\leftarrow {\overline{u}}_B$                                    ▹ Contact scenario is consistent 9: end if 9: return ${\overline{u}}_{\mathrm{sol}}$ as the final solution.

The results of the numerical computations were found to be very close to those of the closed-form solutions, and are shown in Figure 3. Furthermore, the following table provided numerical comparison.