Numerical simulations of the Dynamic Beam Equation in discontinuous media

(1)

MAT-VET-F 20006

Examensarbete 15 hp Juni 2020

Numerical simulations of the Dynamic Beam Equation in discontinuous media

David Niemelä

Valter Wagner Zethrin

Niklas Wik

(2)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Numerical simulations of the Dynamic Beam Equation in discontinuous media

David Niemelä, Valter Wagner Zethrin, Niklas Wik

The study examines the Projection method and the simultaneous-

approximation-term (SAT) method as boundary treatment for the dynamic beam equation using summation-by-parts (SBP) operators for handling the inner domain. The methods are examined for both the homogeneous constant coefficient case, and the inhomogeneous piecewise constant coefficient case with a coupled interface. The outer boundaries are handled by SAT or Projection, the coupled interfaced is handled by Projection or a mix between Projection and SAT. Solutions are integrated in time with finite central difference schemes and compared to analytical solutions. A convergence study is conducted with respect to the spatial discretization to measure the accuracy, and the

stability is examined by numerical simulations of the CFL-condition.

The study shows that Projection has the same accuracy as SAT for most boundary conditions while allowing for a larger timestep. A

discontinuity in the medium is found to be handled equally accurate by Projection and the Projection and SAT mixture for all but one case studied, where the mixture was slightly more accurate.

(3)

Populärvetenskaplig sammanfattning

Numerisk simulering är ett viktigt verktyg för att studera fysikaliska händelseförlopp som är omöjliga eller ekonomiskt ohållbara att studera i experiment. I denna studie undersöks två olika numeriska metoders noggrannhet i att simulera den dynamiska balkekvationen som beskriver hur vibrationer propagerar i bland annat balkar. Den dynamiska balkekvationen är viktig vid design av t.ex. broar och liknande strukturer som är sårbara mot naturfenomen som starka vindar eller jordbävningar. Ekvationen kan även används för att beskriva hur havsvågor kan skapa vibrationer som propagerar i havsisar.

Det här projektets syfte är att undersöka två olika simuleringsmetoder för att se ifall någon av dem är bättre att använda än den andra. Målet med en simulation är att så nära som möjligt återskapa vad som sker i naturen. För detta krävs att simuleringsmetoden har hög noggrannhet. Andra krav på metoder kan vara hur många beräkningar som krävs för att simulera till en viss tidpunkt, det vill säga hur stora tidssteg man kan ta medan lösningen fortfarande överensstämmer med verkligheten. Varför detta är intressant är för att det kan ta flera timmar eller till och med dagar för en dator att simulera ett visst scenario. Om man med en metod kan ta dubbelt så stora tiddssteg som en annan så kan man anta att den totala beräkningstiden halveras.

Då man ibland kräver en viss noggrannhet, till exempel att man är inom några centimeter när man bygger en bro, finns anledning till att undersöka dessa metoders noggranhet och komplexitet. Det kan visa sig att metoderna är olika noggranna, med den ena mycket noggrannare än den andra men att den är mer kom- plex och därför kräver fler beräkningar. Det kan alltså skilja sig vilken metod som kommer att vilja användas.

Studien visar att den största skillnaden mellan metoderna är hur stora tidssteg som tillåts, där Projektions- metoden tillåter störst steg. Noggrannheten för metoderna skiljer sig även något beroende på vilka fall som simuleras, där det varierade vilken av metoderna som var noggrannast.

(4)

1 Introduction

The dynamic beam equation (DBE), also known as the Euler-Bernoulli beam equation, is a standard model in engineering applications. The equation describes the deflection of rigid beams and can be used when mod- elling everything from ice-shelves [4] to the Eiffel tower. The DBE is a linear partial differential equation (PDE) that is second order in time and fourth order in space. It shares similarities with the wave equation but with dispersive characteristics, meaning the propagation speed of waves depend on their wavelengths.

While usually used to describe the deflection of homogeneous 1D beams, the DBE can be applied to piecewise combinations of different materials such as welded metal beams.

See Figure 1 for four different applications involving (or modeled by) DBE:

Figure 1: Examples where DBE is employed. Top left: Eiffel tower, Top middle: Ferris Wheel (London Eye), Top right: Truss bridge (Denmark), Bottom: Ocean wave interactions with floating ice shelves (Flexural gravity waves)

1.1 The dynamic beam equation

The governing equation of the 1-D DBE is given for a beam of length L with its axis along the x-direction, denoting the deflection of the beam from its axis as u(x, t), as

µ(x)∂²u(x, t)

∂t² = − ∂²

∂x²

E(x)I(x)∂²u(x, t)

∂x²

+ F (x, t), 0 ≤ x ≤ L, t ≥ 0 u(x, 0) = f1(x), d

dtu(x, 0) = f2(x), 0 ≤ x ≤ L

(1)

where µ(x) is the mass per unit length, E(x) is the elastic modulus of the beam, I(x) is the second moment of area of the cross section of the beam, F (x, t) is a forcing term and fl,r are initial data.

(6)

By introducing a(x) = E(x)I(x) and a different notation for derivatives in x and t to simplify further notation, (1) becomes,

µ(x)u_tt= − (a(x)u_xx)_xx+ F (x, t), 0 ≤ x ≤ L, t ≥ 0

u(x, 0) = f1(x), ut(x, 0) = f2(x), 0 ≤ x ≤ L (2)

For a homogeneous beam, µ(x) = µ and a(x) = a are independent of x. This important special case was analyzed numerically in [6] using SBP-SAT approximations in space and a finite central difference scheme in time.

1.2 Purpose

The present study intends to extend the numerical analysis to a non-homogeneous beam, in particular when µ(x), E(x) and I(x) have a jump at x = l (for example when combining beams of different materials, i.e, where a(x) and µ(x) are piece-wise constant with a jump at x = l.

(µ(x) = µ1, a(x) = a1, 0 ≤ x ≤ l

µ(x) = µ2, a(x) = a2, l ≤ x ≤ L (3)

Combining summation-by-part (SBP) operators with either Simulatneous Approximation Term (SAT) method or the Projection method is a well-proven, stable methodology for well-posed initial boundary value problems. The SAT and the Projection methods are implemented to impose the boundary conditions and the SBP-operator approximating the inner domain.

The main focus in the present study is to analyse well-posedness for (2) for the case with a jump in the coefficients (3), and implement stable and high-order accurate SBP-Projection and SBP-Projection-SAT, finite difference approximations, and to determine if one of the methods is preferable over the other for handling the inner and outer boundary.

2 Theory

2.1 Well-posedness

A well posed problem is a term coined by the french mathematician Jacques Hadamard, which implied that a problem is well-posed if there exists an unique solution, who’s behaviour continuously changes with the initial conditions. Problems that are not well posed, i.e. ill-posed, are hard to treat numerically, since the solution has to be stable under small changes of the initial data. One example of an ill-posed problem is the inverse heat equation. In fact, a broad variety of inverse problems that arise in physics are ill-posed [8].

α = k

h² (4)

2.2 The energy method

The energy method is used to analyze well-posedness for initial and boundary value problems. By finding an expression for the time dependent energy it is possible to determine the necessary boundary conditions (BC) for a stable solution to exist. In order for a PDE to be well-posed, the energy estimate must be greater or equal to zero. Furthermore the PDE must be energy-stable, i.e. the energy of the system must diminish or remain constant over time. The energy method is a strong tool to show well-posedness for PDE that also has a physical intuition and interpretation.

(7)

2.3 SBP method

The Summation By Parts method is a high order finite difference method, which uses a central finite difference stencil and cautiously chosen one-sided difference stencils to close of the boundaries. The operators are designed to mimic the integration by parts formula in a discrete diagonal norm H, and can be designed to approximate different derivatives. Stable and high order accurate SBP operators are already known and well-proven for well-posed initial boundary value problems, where boundary conditions are imposed by combining the SBP operator with either the simultaneous approximation term (SAT) or the Projection method.

The SBP operators used in this project are referred to by the accuracy of the central scheme. The SBP operators used in the present study were provided by Prof. Ken Mattsson at the Department of Information Technology, Division of Scientific Computing at Uppsala University, and are listed in Appendix. Read [1], [2], [3] for more about Summation By Parts operators and how to construct them for an approximation of desired derivative.

Note that the SBP operators derived in [3], approximating 4^th derivatives, are not the same as the once used in the present study, as these have been updated but not yet published.

2.4 SAT method

Simultaneous approximation term (SAT) is a method to weakly impose boundary conditions to a spatially discretized PDE whilst still preserving its SBP property and energy estimate. The method was first introduced in [10], and shown to be time-stable for any energy-bounded hyperbolic system. The boundary condition is achieved by adding a penalty term to the SBP operator, which is proportional to the difference between the discrete boundary value and the boundary data g(t). In contrast to directly imposing e.g.

a Dirichlet BC u_N = g(t), the SAT method solves a derivative equation on the whole domain, where the extra term accounts for the information at the boundary without altering the accuracy of the internal scheme.

The SAT term is chosen appropriately such that energy stability is retained. Depending on the BC and the PDE itself, the term added varies in complexity. Note that in the present study, SAT is not applied to every BC at the inner boundary of the piece-wise constant coefficient DBE, for the sole reason of its complexity. This was derived at the Department of Information Technology at Uppsala University, however never published. The complexity of the SAT derived brought the idea of the present study, as the Projection method is known to be less complex to derive and to analyse for certain BC and PDE. Recent examples of the SBP-SAT approach can be found in [7] and specifically for the DBE in [4], [6].

2.5 Projection method

The Projection method is a method of imposing boundary conditions on spatially discretizied PDE. The operator projects the numerical solution to a vector space V = {v|Lv = g(t)} where the BC are satisfied, and as such, the Projection operator imposes boundary values completely analogous to the BC in the continuous PDE, allowing it to preserve the SBP property and the energy estimate.

The method was first introduced in [11], [12], with an improved Projection method proposed in [5]. The Projection operator P is defined by

P = I − ¯H⁻¹L^T(L ¯H⁻¹L^T)⁻¹L (5)

where H is the same norm as the SBP operators are based on, L is the boundary operator which represents a discretization of the analytic BC, and where I is the identity matrix.

The Projection method is used in [9] to impose BC for the piece-wise constant coefficient wave equation.

(8)

2.6 Definitions and properties

To simplify the notations later in the study, the inner boundary x = l is placed such that is is equidistant from each outer boundary, x = 0, and x = L. Each domain, left and right of the inner boundary, are discretized using the same number of equidistant grid points,

x^l_i= (i − 1)h, i = 1, 2, ..., N

x^r_j = (i − 1)h, j = N, N + 1, ..., 2N − 1, h = L − l N − 1 = l

N − 1

(6)

note that x^l_N = x^r_N denotes the same grid point at the interface x = l. The approximate solution left of the inner boundary at grid point x^l_i, is denoted ui, and the discrete solution vector left of the inner boundary is denoted u^T = [u₁, u₂, ..., u_N]. Similarly for the grid points right of the inner boundary, the approximate solution for grid point x^r_j is denoted v_j, and the discrete solution vector v = [v_N, v_n+1, ..., v_{2N −1}]^T. The solution for the whole interval 0 ≤ x ≤ L can therefore be written as w^T = [u^T, v^T].

Let e₁ = [1, 0, ..., 0]^T and e_N = [0, ..., 0, 1]^T be (N × 1) vectors. It follows that e^T₁u = u₁, e^T₁v = v_N, e^T_Nu = uN and e^T_Nv = v2N −1.

The difference operator,

D4= H⁻¹(−M + BS) = H⁻¹(−M − e^T₁d3:1+ e^T_Nd3:N+ d^T_1:1d2:1− d^T_1:Nd2:N) (7) approximating ∂⁴/∂x⁴, using a pth-order accurate narrow-stencil in the interior, is said to be a pth order diagonal-norm fourth-derivative SBP operator if H is a diagonal matrix that defines a discrete norm, M is positive semi-definite, e1v ' u|^l, eNv ' u|^r, d^T_1:1v ' ux|^l, d^T_1:Nv ' ux|^r, d^T_2:1v ' uxx|^l, d^T_2:Nv ' uxx|^r, d^T_3:1v ' uxxx|^l, d^T_3:Nv ' uxxx|^r, are finite difference approximations of the first, second and third derivatives at the left and right boundaries on the domains 0 ≤ x ≤ l and l ≤ x ≤ L of the exact solution.

The dissipative part M of a diagonal-norm fourth-derivative SBP operator has the following properties:

v^TM V = hα2((d2:1v)²+ (d2:Nv)²) + v^TM˜2V

v^TM V = h³α3((d3:1v)²+ (d3:Nv)²) + v^TM˜3V (8) where ˜M2,3 are symmetric and positive semi-definite, and α2,3 are positive constants, independent of the discretization step size h. The values of α2,3 were numerically derived in [6] for the second-, fourth- and sixth-order accurate finite difference SBP operators

α^(2nd)₂ α^(2nd)₃ α^(4nd)₂ α^(4nd)₃ α^(6nd)₂ α^(6nd)₃ 1.250 0.4 0.548 1.088 0.322 0.156

3 Analysis

The focus of the analysis is towards a case of the DBE in which there is an inner boundary, i.e where there is a jump in the coefficients. The section contains analysis of well-posedness of the DBE, where inner and outer boundary conditions are found through the use of the energy method. The outer BC used are listed in (12), and the general inner BC in (13).

Moreover, analytical solutions are derived for the DBE as formulated in (2) where the forcing term is set to zero and the coefficients are piece-wise constant. The solutions are found using the separation of variables

(9)

Finally, stability of the semi-discrete DBE with piece-wise constant coefficients for the different finite difference methods used in the paper is shown by the use of the energy method. The combinations analysed include the semi-discrete SBP-Projection approximation, where all BC are treated with Projection, the semi- discrete SBP-Projection-SAT approximation, where the higher order BC terms at the inner boundary are treated with SAT, and the two lower order terms are treated with Projection. Here, similar analysis can be made for the different mixtures of Projection and SAT at the inner boundary. The mixtures left out from the analysis are where only the highest or the second highest order term is treated with SAT, and the three remaining terms are treated with Projection. The outer BC are analysed for both SAT and Projection.

Note that similar analysis for stability of the semi-discrete DBE and derivation of analytical solutions also can be made for the case where the coefficients are constant over the whole spacial domain. These are left out of the analysis since they easily follow from the case with piece-wise constant coefficients.

3.1 Well-posedness of the DBE

Let the inner product for two real valued function u, v ∈ L²[0, L] be defined by (u, v) = RL

0 dx, and let (u, u) be the corresponding norm ||u||². Well-posedness of (2), with piece-wise continuous and differentiable coefficients al(x) ∈ C¹[0, l], ar(x) ∈ C¹[l, L], with an internal boundary at x = l, is analysed using the energy method for the homogeneous case, where the forcing term F (x, t) and the outer boundary data are set to zero.

Multiplying by u_t and piece-wise integrating by parts twice on the right hand side, known as the energy method, gives

(u_t, µ(x)u_tt) = − u_t(a_l(x)u_xx)_x

l

0− ut(a_r(x)u_xx)_x

L l

+ u_txa_l(x)u_xx

l 0

+ u_txa_r(x)u_xx

L

l

− (utxx, a(x)u_xx)

(9)

Further, adding the transpose leads to d

dtE =2

−ut(a(x)uxx)x

l

0− ut(a(x)uxx)x

L

l + utxa(x)uxx

l

0+ utxa(x)uxx

L l

(10) where the continuous energy E is defined by

E = ||µ(x)u_t||²+ ||a(x)u_xx||² (11)

For E to be an energy, E ≥ 0 is required, implying that µ(x) and a(x) have to be greater or equal to zero.

This corresponds to the physical interpretation that the mass per unit length µ, as well as the second moment of area of the cross section, and the elastic modulus of the beam are positive quantities. For the problem to be well-posed the energy is required to be conserved or decreasing over time, i.e _dt^dE ≤ 0. Therefore (10) results in restrictions to all boundaries, where 2 boundary conditions are needed at x = 0 and x = L respectively. The most commonly used BC, and which can be physically interpreted, are listed below,

u = 0, ux= 0, ”Clamped”

u = 0, u_xx= 0, ”Hinged”

ux= 0, uxxx= 0, ”Sliding”

uxx= 0, uxxx= 0, ”F ree”

(12)

(10)

Furthermore, well-posedness for (2) is satisfied with continuity at the inner boundary x = l. The boundary conditions for continuous media therefore reads,

lim

x−→lu = lim

x₊→lu lim

x−→lux= lim

x+→lux

lim

x−→la(x)u_xx= lim

x₊→la(x)u_xx lim

x−→l(a(x)uxx)x= lim

x+→l(a(x)uxx)x

(13)

Introducing an internal boundary at x = l is only relevant for a case where there is a jump in the coefficient a(x). Such a case is relevant to many different applications and has been done before for the acoustic wave equation, see [9].

Choosing a piece-wise constant coefficient a(x) with a discontinuity at x = l, i.e, a(x) = a₁, − 1 ≤ x ≤ 0

a(x) = a2, 0 ≤ x ≤ 1 (14)

The boundary conditions at the outer boundary are unchanged, and the inner boundary conditions become, lim

x−→lu = lim

x+→lu lim

x−→lu_x= lim

x₊→lu_x a1 lim

x−→luxx= a2 lim

x+→luxx

a₁ lim

x−→lu_xxx= a₂ lim

x₊→lu_xxx

(15)

3.2 Analytical solution to the DBE

Analytical solutions to (2) are found with the separation of variables method. However, only analytical solutions for special cases of boundary data and initial conditions are derived. For every case studied, the forcing term F (x) and the outer boundary data are set to zero. Furthermore, the material dependent constant a(x) is chosen to be piece-wise constant, with a discontinuity at x = l as in (14), and the mass per unit length µ(x) to be constant over the whole domain. The BC considered for the analytical solutions are the commonly used BC mentioned in (12).

To simplify the calculations, the domain of (2) is changed and the equation is rewritten as

utt= −a1uxxxx, −1 ≤ x ≤ 0, t ≥ 0

u(x, 0) = f₁(x), u_t(x, 0) = f₂(x), −1 ≤ x ≤ 0

vtt= −a2vxxxx, 0 ≤ x ≤ 1, t ≥ 0

v(x, 0) = g1(x), vt(x, 0) = g2(x), 0 ≤ x ≤ 1

(16)

where u and v are the solutions to the DBE on −1 ≤ x ≤ 0 and 0 ≤ x ≤ 1 respectively, a1 and a2 are constants, and f1(x), f2(x), g1(x) and g2(x) are initial data.

Separation of variables gives; u(x, t) = Tu(t)X_u(x) and v(x, t) = T_v(t)X_v(x), where T_u(t), T_v(t) and X_u(x),

(11)

Xv(x) are functions of only t and x respectively. Inserting in (16), gives the following equations:

(T_u)_tt Tu

= −a₁(X_u)_xxxx Xu

= λ₁, −1 ≤ x ≤ 0, t ≥ 0 (Tv)tt

T_v = −a₂(Xv)xxxx

X_v = λ₂, 0 ≤ x ≤ 1, t ≥ 0 (17)

where λ1and λ2are constants, since the right and left hand sides are functions of independent variables and therefore has to be constant. Rewriting (17) results in the following system of ODE:s

(Tu)tt− Tuλ1= 0, t ≥ 0

− a₁(X_u)_xxxx− X_uλ₁= 0 −1 ≤ x ≤ 0 (Tv)tt− Tvλ2= 0, t ≥ 0

− a2(Xv)xxxx− Xvλ2= 0 0 ≤ x ≤ 1

(18)

The general solutions to (18) are:

T_u= A₁cosp λ₁t

+ A₂sinp λ₁t

, t ≥ 0

Xu= B1sinh(b1x) + B2sin(b1x) + B3cosh(b1x) + B4cos(b1x) −1 ≤ x ≤ 0 Tv = A3cosp

λ2t

+ A4sinp λ2t

, t ≥ 0

X_v= B₅sinh(b₂x) + B₆sin(b₂x) + B₇cosh(b₂x) + B₈cos(b₂x) 0 ≤ x ≤ 1

(19)

where A1, A2, A3, A4, B1, B2, B3, B4, B5, B6, B7, B8 are constants, b1 =q

λ1

√a₁ and b2 = q

λ2

√a₂. As- suming f2(x) = g2(x) = 0, it follows that A2 = A4 = 0. Furthermore, for simplicity, f1(x) = Xu(x) and g1(x) = Xv(x) is assumed, which implies that A1= A3= 1.

Applying the continuity conditions (15), at the inner boundary for (19) gives the following correlations;

√a1b²₁−√

a2b²₂= 0 B₃+ B₄− A7− A8= 0 b1(B1+ B2) − b2(B5+ B6) = 0 a1b²₁(B3− B4) − a2b²₂(B7− B8) = 0 a₁b³₁(B₁− B2) − a₂b³₂(B₅− B6) = 0

(20)

Depending on the outer BC chosen, further correlations between the unknown constants vary. For BC

"clamped" at x = −1 and x = 1 the following correlations are found;

− B1sinh(b1) − B2sin(b1) + B3cosh(b1) + B4cos(b1) = 0 b1(B1cosh(b1) + B2cos(b1) − B3sinh(b1) + B4sin(b1)) = 0 B5sinh(b2) + B6sin(b2) + B7cosh(b2) + B8cos(b2) = 0 b2(B5cosh(b2) + B6cos(b2) + B7sinh(b2) − B8sin(b2)) = 0

(21)

For BC "free" at x = −1 and x = 1 the following correlations are found;

b²₁(−B1sinh(b1) + B2sin(b1) + B3cosh(b1) − B4cos(b1)) = 0 b³₁(B₁cosh(b₁) − B₂cos(b₁) − B₃sinh(b₁) − B₄sin(b₁)) = 0 b²₂(B5sinh(b2) − B6sin(b2) + B7cosh(b2) − B8cos(b2)) = 0 b³₂(B5cosh(b2) − B6cos(b2) + B7sinh(b2) + B8sin(b2)) = 0

(22)

(12)

For BC "sliding" at x = −1 and x = 1 the following correlations are found;

b₁(B₁cosh(b₁) + B₂cos(b₁) − B₃sinh(b₁) + B₄sin(b₁)) = 0 b³₁(B1cosh(b1) − B2cos(b1) − B3sinh(b1) − B4sin(b1)) = 0 b₂(B₅cosh(b₂) + B₆cos(b₂) + B₇sinh(b₂) + B₈sin(b₂)) = 0 b³₂(B₅cosh(b₂) − B₆cos(b₂) + B₇sinh(b₂) + B₈sin(b₂)) = 0

(23)

For BC "hinged" at x = −1 and x = 1 the following correlations are found;

− B1sinh(b₁) − B₂sin(b₁) + B₃cosh(b₁) + B₄cos(b₁) = 0 b²₁(−B1sinh(b1) + B2sin(b1) + B3cosh(b1) − B4cos(b1)) = 0 B5sinh(b2) + B6sin(b2) + B7cosh(b2) + B8cos(b2) = 0 b²₂(B₅sinh(b₂) − B₆sin(b₂) + B₇cosh(b₂) − B₈cos(b₂)) = 0

(24)

Combining (20) with either one of (21)-(24) multiple non-trivial solutions can be found. Listed in Tables3 and4in the Appendix are a few different solutions found by solving f (x) = 0 in MATLAB with the function fsolve, where x = [b1, b2, B1, B2, B3, B4, B5, B6, B7, B8], and f (x) = 0 corresponds to all the correlations found between the unknown constants.

(13)

(a) a1 = 1 and a2 = 100, clamped, T = 0. (b) a1 = 1 and a2 = 100, clamped, T = 0.09.

(c) a1 = 1, a2 = 4, free, T = 0. (d) a1 = 1, a2 = 4, free, T = 0.14.

Figure 2: Analytic solution for BC "clamped" and "free" for different time T and constants a1 and a2

3.3 Stability analysis using pure Projection

The semi-discrete SBP-Projection approximation of the DBE as formulated in (2), with forcing term F (x) and outer boundary data fl(x) and fr(x) set to zero, and with constant coefficient µ, can be written as,

µw_tt= − ¯D₄^(a)w, t ≥ 0 Llw = gl(t), t ≥ 0 L_rw = g_r(t), t ≥ 0 Lcw = gc(t), t ≥ 0 w = f1, wt= f2, t = 0

(25)

where L_l, L_rand L_c are the semi-discrete boundary operators, g_l, g_r and g_c are the boundary data for the left, right and centre boundaries respectively, f is the initial data, and ¯D4

(a)is,

D¯^(a)₄ = ¯H⁻¹(− ¯M^(a)+ ¯BS^(a)) =

H⁻¹ 0 0 H⁻¹

−a1M 0 0 −a2M

+

a1BS 0 0 a₂BS

(26) Note that a non-constant coefficient µ that depends on x is allowed. However, for simplicity, µ was chosen to be constant.

(14)

Introducing the full boundary operator from the corresponding boundary data,

g(t) =



 g_l(t) g_r(t) g_c(t)



 (27)

let’s the semi-discrete boundary operator to be written as,

Lw = g(t) (28)

Applying the semi-discrete approximations of the continuity conditions (15) found for the continuous case gives the central boundary operator,

Lc =







e^T_N −e^T₁ d^T_1,N −d^T_1,1 a1d^T_2,N −a2d^T_2,1 a₁d^T_3,N −a₂d^T_3,1







(29)

The left and right semi-discrete boundary operators with the semi-discrete approximation of the other boundary conditions (12) are given by,

Ll=

e^T₁ 0 d_1:1 0

, Lr=

0 e^T_N 0 d^T_1:N

, ”Clamped”

Ll=

e^T₁ 0 d^T_2:1 0

, Lr=

0 e^T_N 0 d^T_2:N

, ”Hinged”

Ll=

d^T_1:1 0 d_3:1 0

, Lr=

0 d^T_1:N 0 d^T_3:N

, ”Sliding”

Ll=

d^T_2:1 0 d^T_3:1 0

, Lr=

0 d^T_2:N 0 d^T_3:N

, ”F ree”

(30)

The Projection matrix is defined by (5), where the boundary operator L is on the form L = [L^T_l , L^T_c, L^T_r]^T. The improved Projection method, see [5], with homogeneous boundary data gives,

µwtt= −P ¯D^(a)₄ P w, t ≥ 0

Lw = 0 t ≥ 0

w = f1, wt= f2, t = 0

(31)

Multiplying by w^T_tH from the left and adding the transpose gives,¯

w_t^THµw¯ tt+ wtHµw¯ ^T_tt= w_t^THP ¯¯ H⁻¹(− ¯M^(a)+ ¯BS^(a))P w + bwtH(P ¯¯ H⁻¹(− ¯M^(a)+ ¯BS^(a))P w)^T (32) simplifying this further using the properties of the matrices ¯H and P , ( ¯HP = P^TH), gives¯

d

dt||µwt||H¯ = (wtP )^T(− ¯M^(a)+ ¯BS^(a))P w + wtP (− ¯M^(a)+ ¯BS^(a))(P w)^T (33) d

dt

||µw_t||H¯ + (P w)^TM¯^(a)P w

= (P w_t)^TBS¯ ^(a)P w + ( ¯BS^(a)P w)^TP w_t (34)

Numerical simulations of the Dynamic Beam Equation in discontinuous media

Examensarbete 15 hp Juni 2020

Numerical simulations of the Dynamic Beam Equation in discontinuous media

David Niemelä

Valter Wagner Zethrin

Niklas Wik

Abstract

Numerical simulations of the Dynamic Beam Equation in discontinuous media

Populärvetenskaplig sammanfattning

Contents

1 Introduction

2 Theory

3 Analysis