Hello! Yes, the Newmark-beta method can be used in conjunction with the finite difference method (specifically central difference) for the spatial nodes to accurately solve the Euler-Bernoulli beam equation with fixed boundary conditions. In fact, this is a common approach for solving dynamic structural problems, including those involving beams and other slender structures.

To provide a more detailed explanation, let's first review the Newmark-beta method and the Euler-Bernoulli beam equation.

The Newmark-beta method is a time integration scheme commonly used in structural dynamics to solve dynamic equations of motion. It is an implicit method, which means that the future state of the system is computed based on the current state and the previous state. The Newmark-beta method has two parameters: the beta parameter and the gamma parameter. These parameters control the accuracy and stability of the method. The algorithm for the Newmark-beta method is as follows:

The Euler-Bernoulli beam equation describes the bending behavior of a beam subjected to external forces. It is a fourth-order partial differential equation and is given by:

EI * w''''(x) + q(x) = 0

where EI is the flexural rigidity, w(x) is the transverse displacement of the beam at position x, and q(x) is the distributed load.

To solve the Euler-Bernoulli beam equation with the Newmark-beta method, we need to discretize the beam into a series of nodes and apply the finite difference method to approximate the second and fourth derivatives of the displacement with respect to the position. Specifically, we can use the central difference method to discretize the beam as follows:

w_{i-1} - 2*w_i + w_{i+1} = h^2 * (-q_i/EI)

where w_i is the transverse displacement at node i, q_i is the distributed load at node i, h is the distance between nodes, and EI is the flexural rigidity.

This equation can be rearranged to give an expression for the second derivative of the displacement at node i:

w''i = (w{i-1} - 2*w_i + w_{i+1}) / h^2 - q_i/EI

Similarly, we can use the central difference method to approximate the fourth derivative of the displacement:

w''''i = (w{i-2} - 4w_{i-1} + 6w_i - 4*w_{i+1} + w_{i+2}) / h^4

Substituting these expressions into the Euler-Bernoulli beam equation, we obtain a system of coupled second-order ordinary differential equations that can be solved using the Newmark-beta method.

As an example, consider a simply supported beam of length L subjected to a distributed load q(x) = 1000 N/m and with flexural rigidity EI = 10^7 Nm^2. The beam is discretized into 11 nodes, with a distance of h = L/10 between nodes. The boundary conditions are w(0) = 0 and w(L) = 0, representing fixed supports.

To solve for the transverse displacement of the beam at each time step, we need to first assemble the stiffness matrix K and the mass matrix M. For the Euler-Bernoulli beam equation, the stiffness matrix is given by:

K_ij = EI * (w''_i * w''_j) + (w''''_i * w''''_j)

And the mass matrix is given by:

M_ij = rho * A * h * (delta_ij + delta_{i,j+1} + delta_{i,j-1}) / 3

Where rho is the density of the beam material, A is the cross-sectional area, and delta_ij is the Kronecker delta.

Once we have assembled the matrices, we can use the Newmark-beta method to solve for the displacement and velocity at each time step. To do this, we need to compute the acceleration at the current time step, a_n, using the equation of motion:

a_n = M^(-1) * (F_n - K * u_n - C * v_n)

Where F_n is the external force vector at time step n, and C is the damping matrix (which we can assume is zero for this example).

For the distributed load q(x) = 1000 N/m, the external force vector at node i is given by:

F_i = q_i * h / 2

We can then use the Newmark-beta algorithm to compute the displacement and velocity at the next time step, u_{n+1} and v_{n+1}, based on the current displacement and velocity, u_n and v_n, and the acceleration at the current and next time steps, a_n and a_{n+1}.

Here is a simple Python code that implements the finite difference method and the Newmark-beta method to solve the Euler-Bernoulli beam equation with fixed boundary conditions:

========================

import numpy as np

import matplotlib.pyplot as plt

# Define parameters

L = 1.0 # Beam length

h = 0.1 # Node spacing

rho = 7800 # Density of steel

A = 0.01 # Cross-sectional area

EI = 1e7 # Flexural rigidity

q = 1000 # Distributed load

# Create node coordinates

x = np.linspace(0, L, int(L/h)+1)

# Assemble stiffness matrix

K = np.zeros((len(x), len(x)))

for i in range(1, len(x)-1):

K[i, i-1] = 1

K[i, i] = -2

K[i, i+1] = 1

K *= EI/h**3

for i in range(2, len(x)-2):

K[i, i-2] = 1

K[i, i+2] = 1

K *= EI/h

# Assemble mass matrix

M = np.zeros((len(x), len(x)))

for i in range(1, len(x)-1):

M[i, i-1] = 1

M[i, i] = 4

M[i, i+1] = 1

M *= rho*A*h/6

# Compute time step and total time

dt = 0.01

T = 1.0

# Set initial conditions

u = np.zeros(len(x))

v = np.zeros(len(x))

# Define Newmark-beta parameters

gamma = 0.5

beta = 0.25

# Compute constants

a0 = 1 / (beta*dt**2)

a1 = gamma / (beta*dt)

a2 = 1 / (beta*dt)

a3 = 1 / (2*beta) - 1

a4 = gamma / beta - 1

a5 = dt/2 * (gamma/beta - 2)

# Time integration loop

for t in np.arange(0, T, dt):

# Compute external force vector

F = np.zeros(len(x))

for i in range(1, len(x)-1):

F[i] = q*h/2

# Compute acceleration

a = np.linalg.solve(M, F - np.dot(K, u))

# Compute displacement and velocity at next time step

u_new = u + dt*v + dt**2*(a0*a - a2*v - a3*u)

v_new = v + dt*(a1*a - a4*v - a5*u)

# Apply boundary conditions

u_new[0] = 0

u_new[-1] = 0

v_new[0] = 0

v_new[-1] = 0

# Update state

u = u_new

v = v_new

# Plot displacement at current time step

plt.plot(x, u)

plt.xlabel('Beam length (m)')

plt.ylabel('Transverse displacement (m)')

plt.show()

===============================

This code first assembles the stiffness and mass matrices using the finite difference method for spatial discretization. It then uses the Newmark-beta method to compute the displacement and velocity at each time step based on the current state and the previous state, and applies fixed boundary conditions. Finally, the code plots the displacement of the beam at each time step over the full length of the beam.

Note that this code is a very simple example and does not include error checking, convergence criteria, or other advanced features. However, it should be enough to give you an idea of how to implement the Newmark-beta method and finite difference method for the Euler-Bernoulli beam equation in Python.

To run this code, simply copy and paste it into a Python environment such as Jupyter Notebook or Spyder, and run it. You should see a plot of the transverse displacement of the beam at each time step, which will show how the beam responds to the distributed load over time. You can adjust the parameters such as beam length, node spacing, and time step to see how they affect the behavior of the beam.

Alessandro Rizzo thank you for your reply. What if the beam equation is dynamic, i.e EI * w''''(x) + rho*A*w''(t) + q(x) = 0 where w''(t) is the second order time derivative of the deflection w. Would the same working principle for the Newmark-beta method in conjunction with the finite difference apply?

You are welcome Uchenna Ogunka,

Yes! the same working principle for the Newmark-beta method in conjunction with the finite difference method can be applied to the dynamic Euler-Bernoulli beam equation that you have provided:

EI * w''''(x) + rhoAw''(t) + q(x) = 0

This equation accounts for the inertial effects of the beam by introducing the term rhoAw''(t), where rho is the density of the material, A is the cross-sectional area, and w''(t) is the second-order time derivative of the deflection w.

To solve this dynamic beam equation using the Newmark-beta method and the finite difference method, you can follow a similar approach as in the previous explanation. Here's a brief summary of the steps:

M * w''(t) + K * w(x) = -q(x)

where w''(t) is the second-order time derivative of the deflection w, and w(x) is the displacement vector at each node.

The main difference between solving the static and dynamic beam equations is the inclusion of the inertial term and the use of the mass matrix M when computing the acceleration at each time step. The overall approach for using the Newmark-beta method in conjunction with the finite difference method remains the same.

Keep in mind that the provided Python code in the previous answer should be modified to account for the dynamic beam equation. Specifically, you will need to adjust the equation of motion to include the inertial term when computing the acceleration at each time step ;)

Hello, Alessandro Rizzo , and thank you for this detailed explanation!

Hope it is allowed to ask some followup questions here. I have tried to implement your code, but I get a singular matrix error. Does this code run for you? I notice that the first and last row of both K and M are zeros, and I assume that is causing the error. I'm having trouble understanding how the matrices are derived, and all I can find online about stiffness matrices is regarding Finite Elements. Could you explain a bit more how these are built, coming from your equations? Or is there some documentation you can point me the way to?

Once again, thanks for the explanation!

Hello Petter Nøst,

Yes, it is common to have a singular matrix when solving a boundary value problem with fixed boundary conditions, as the nodes corresponding to these boundary conditions have zero displacements. You can avoid the singular matrix error by adjusting the system of equations for the given boundary conditions.

Since the first and last rows of both K and M are zeros, you should apply the fixed boundary conditions by removing the rows and columns corresponding to the fixed nodes. In the case of the Euler-Bernoulli beam equation with fixed boundary conditions at both ends, you should remove the first and last rows and columns from the matrices K and M.

Here's a quick overview of how to modify the provided code to apply the fixed boundary conditions:

1. Remove the first and last rows and columns from the matrices K and M.

K = K[1:-1, 1:-1] M = M[1:-1, 1:-1]

2. Adjust the force vector q(x) accordingly to account for the reduced size of the matrices K and M.

q = q[1:-1]

3. Modify the boundary conditions in the solution array to include only the internal nodes, i.e., excluding the fixed nodes. w = w[1:-1]

After making these changes, you should be able to avoid the singular matrix error and solve the dynamic Euler-Bernoulli beam equation using the Newmark-beta method and the finite difference method.

As for the derivation of the stiffness and mass matrices, they are obtained from the discretization of the beam equation and the assembly of the local element equations. The approach to derive these matrices is different from the Finite Element Method, but the underlying concept of discretizing the governing equation and assembling the global system of equations remains the same.

In the finite difference method, the second and fourth spatial derivatives of the displacement are approximated using the central difference scheme. This discretization results in a set of local element equations that are assembled to form the global system of equations, which are then represented in the form of the stiffness matrix K and mass matrix M.

For more information on how to derive the stiffness and mass matrices for beam problems using the finite difference method, you can refer to the following resources:

I hope this could be helpful with your implementation, Good luck!

Hello Alessandro Rizzo

I understand that the concept of using a mass and stiffness matrix for finite element method. Can the mass matrix and stiffness matrix be used with finite difference method? Discretizing The non-dimensional Euler Bernoulli eqn: w''''(x) + w''(t) + F = 0 gives

(w[i-2] - 4w[i-1] + 6w[i] -4w[i+1] + w[i+2])/(h^4) + (w[i-1] - 2w[i] + w[i+1])/(k^2) + F[i] = 0. h is spatial step size and k is the time step size. How would the mass matrix and stiffness matrix be implemented in this case?

Hello Uchenna Ogunka,

Yes, the Newmark-beta method can be used in conjunction with the finite difference method to solve the 4th order Euler-Bernoulli beam equation. Here is a basic breakdown of how it can be done:

The Euler-Bernoulli beam equation describes the behavior of a slender beam under load and is given by:

EI * w''''(x) = q(x)

where E is the elastic modulus, I is the area moment of inertia, w is the deflection of the beam, and q is the distributed load along the beam.

First, you would discretize the beam into a finite number of nodes and use the central difference method to approximate the fourth derivative of the deflection, w''''(x), at each node. This gives you a system of second order ordinary differential equations.

The Newmark-beta method is used to solve these second order ordinary differential equations. It's a time integration scheme that allows you to solve for the deflection of the beam at each time step. It's an implicit method, which means the calculation for the new time step involves the solution at the new time step itself.

Here's a basic outline of the steps involved in implementing this method:

1. Begin by discretizing the beam into a series of nodes along its length.

2. Use the central difference method to create a system of second order ordinary differential equations that describe the behavior of the beam.

3. Choose initial conditions for the deflection and velocity of the beam at each node.

4. Use the Newmark-beta method to integrate these equations forward in time and solve for the deflection and velocity of the beam at each node for each subsequent time step.

Remember to ensure that the solution at each time step satisfies the boundary conditions of the beam. For a beam with fixed ends, the deflection and slope at the ends of the beam are zero, and these conditions should be enforced at each time step.

Keep in mind, the stability and accuracy of the solution can be influenced by the choice of time step and the parameters used in the Newmark-beta method, so you may need to experiment with these to get a solution that meets your needs.

Regards,

Thank you for your message Alessandro Rizzo

I have implemented the central difference scheme with the implicit Newmark beta method. However, I observed that it is conditionally stable after running tests with different time step sizes.

Hello Uchenna Ogunka,

It's great to hear that you've implemented the central difference scheme with the implicit Newmark-beta method.

The Newmark-beta method is a time integration method used primarily in structural dynamics and earthquake engineering for solving differential equations. It is a one-step method that uses information about the solution and its derivatives at the current and previous time steps.

In theory, when using the Newmark-beta method with certain parameters, it is unconditionally stable. This is one of its key advantages. Specifically, the method is unconditionally stable if the following conditions are met:

- β ≥ 0.25

- γ ≥ 0.5

where β and γ are the Newmark method parameters.

If you're finding that your implementation of the Newmark-beta method is conditionally stable, it's possible that the parameters you've chosen do not meet these conditions.

That being said, the central difference method is known to be conditionally stable. This means that the stability of the solution is dependent on the time step size, and this is probably the issue you're seeing.

In finite difference schemes, the stability condition (or Courant-Friedrichs-Lewy condition) often ties the time step size to the spatial step size. For your case, if you're seeing conditional stability, you might need to adjust your time step size based on your spatial discretization to ensure stability.

Furthermore, it's also possible that the source of your stability issues could be coming from your spatial discretization. If your mesh is not sufficiently fine, or if it's not correctly capturing the physics of the problem, this could also lead to instability.

Best of luck with your work!

Yours Alessandro