Vehicle-Bridge Interaction System
Introduction to IDE (1)Yamamoto Kyosuke, Asst. Prof.
ver. 2.2022.1011

Road unevenness induces
vibrations on a travelling vehicle

The simplest model of the travelling vehicle isa mass-spring modelSuspensions are usually modeled byspring and damping

Equation of motion of the vehicle: (𝑚𝑎 = 𝐹)
𝑚: Mass𝑚𝑧 = −𝑐 𝑧 − 𝑟 − 𝑘 𝑧 − 𝑟
Vehicle Vibration

damping force

restoring force

𝑘: stiffness
𝑐: damping
Linear Differential Equation
Road Profile

System Parameters: 𝑚, 𝑐, 𝑘

Output: 𝑧 𝑡

Input: 𝑟 𝑡𝑑𝑑𝑡𝑑=𝑑𝑡=Road Profile: 𝑟(𝑡)Displacement (m)Time (s)System Parameters:

Displacement (m)Mass: 𝑚Damping: 𝑐Stiffness: 𝑘

10000 [kg]2000 [kg/s]

10000 [N/m]

Vehicle Vibration: 𝑧 𝑡Time (s)%-- set vehicle parametersfprintf('setting vehicle parameters...¥n')m = 10000; %-- mass (kg)c = 2000; %-- damping (kg/s)
k = 10000; %-- stifness (N/m)=(kg/s/s)
v = 10.0; %-- speed (m/s)
dt= 0.001; %-- time increment;
T = 60;%-- end time%-- calculate input
fprintf('calculating vehicle input...¥n')
R = load('profile-1.txt'); %-- load: R(x)x = 0:0.20:(length(R)-1)*0.20; %-- dx=0.20(m) for R(x)t = 0:dt:T; %-- time, dt=0.001(s)xx= v*t; %-- vehicle position: x(t)r = spline(x,R,xx); %-- spline interpolation: R(x) to r(t)=R(x(t))
dr= 1/dt*gradient(r); %-- numerical differential
r=r-r(1); %-- reference value setting: r(0)=0
dr=dr-dr(1); %-- reference value setting: dr(0)=0
f=c*dr+k*r; %-- right-hand-side of equation of motion of the vehicle%-- simulate output
fprintf('vehicle vibration simulation: newmark-beta method...¥n')G=1/2; %-- gammma coefficient for newmark-beta method
B=1/4; %-- beta coefficient for newmark-beta method[ t, z, dz, ddz ] = newmark_beta( m,c,k,f,1/1000,zeros(1,3), G, B );function [ t, x, dx, ddx ] = newmark_beta( m,c,k,f,dt,X, G, B )%% initialization%-- memory space
x=zeros(size(f));
dx=zeros(size(f));
ddx=zeros(size(f));%-- input initial valuesx(:,1)=X(:, 1);dx(:,1)=X(:, 2);ddx(:,1)=X(:, 3);%-- time vector
T=length(f);t=(0:T-1)*dt;%% solving
A = m + dt*G*c + dt^2*B*k; %-- global matrix
M = inv(A); %-- inverse matrix of A
for tt=2:T%-- solving preparation
b1 = -c * ( dx(:,tt-1) + (1-G)*dt*ddx(:, tt-1) );b2 = -k * ( x(:,tt-1) + dt*dx(:,tt-1) + (1/2-B)*dt^2*ddx(:,tt-1) );b = f(:,tt) + b1 + b2; %-- right hand side%-- visualization
fprintf('outputting the simulation results...¥n')figure(1)subplot(2,1,1); plot(t,r,'-b')subplot(2,1,2); plot(t,z,'-r')%-- animation
figure(2);set(gcf, 'Unit', 'pixel', 'Position',[120, 280, 800, 300]);
writerObj = VideoWriter('ex_01.avi','Motion JPEG AVI'); %-- setting movie
classwriterObj.FrameRate = 20; %-- flame rate (slide/sec)open(writerObj); %-- initialization of movie objectfor tt=1:100:length(t)plot(xx(1:tt)',[r(1:tt);z(1:tt)+0.01]');hold onplot(xx(tt),r(1,tt),'bo');plot(xx(tt),z(1,tt)+0.01','ro');hold off
xlim([0 max(xx)]);ylim([-0.004 0.020]);
frame = getframe(figure(2)); %-- get the figure as a frame of theanimation
writeVideo(writerObj,frame); %-- add the gotten frame to the movie object
end
close(writerObj); %-- end of animation%-- calculation of acceleration at t=t+1ddx(:,tt) = M*b;%-- calculation of velocity and displacement responsesdx(:,tt) = dx(:,tt-1) + dt*(1-G)*ddx(:,tt-1) + dt*G*ddx(:,tt);x(:,tt) = x(:,tt-1) + dt*dx(:,tt-1) + dt^2*(1/2-B)*ddx(:,tt-1) + dt^2*B*ddx(:,tt);
end

end

Equation of motion of the vehicle: (𝑚𝑎 = 𝐹)𝑚 𝑠 = −𝑐 𝑠 − 𝑧 − 𝑘 𝑠 − 𝑧

𝑚 : sprung-mass
Sprungmass
Vibration𝑘 : sprungstiffness𝑚 : unsprungmass𝑘 : unsprungstiffness𝑐 : sprungdamping

𝑚 𝑧 =𝑐 𝑠−𝑧 +𝑘 𝑠−𝑧 −𝑘 𝑧−𝑟

𝑚 𝑠+𝑐 𝑠+𝑘 𝑠 =𝑐 𝑧+𝑘 𝑧
𝑚 𝑧+𝑐 𝑧+ 𝑘 +𝑘 𝑧 =𝑘 𝑟+𝑐 𝑠+𝑘 𝑠

UnsprungmassVibration𝑚𝑚𝑐𝑠+ −𝑐𝑧−𝑐𝑐

𝑘𝑠+𝑧−𝑘−𝑘𝑘 +𝑘Road Profile

𝐌𝒛 + 𝐂𝒛 + 𝐊𝒛 = 𝒇0𝑠=𝑘 𝑟𝑧

Matrix and Vector
System of Equations:

Using Matrix:

𝑎 𝑥 +𝑎 𝑥 =𝑏𝑎𝑎𝑎𝑎𝑥𝑥

𝑏=𝑏𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑥𝑥𝑥𝑎 𝑥 +𝑎 𝑥 =𝑏𝑎 𝑥 +𝑎 𝑥 +𝑎 𝑥 =𝑏𝑎 𝑥 +𝑎 𝑥 +𝑎 𝑥 =𝑏𝑎 𝑥 +𝑎 𝑥 +𝑎 𝑥 =𝑏𝑏= 𝑏𝑏Solution:
Inverse matrix

1) Download: http://www.kz.tsukuba.ac.jp/~yamamoto_k/material/matlab_for_ide.zip
2) Decompress the downloaded zip file3) Launch Matlab and set the decompressed directory as current directory4) Run the following script:>> sample_02

Two quarter-cars can be easily
extended to a half-car model𝑑𝑚𝑠𝑠

𝑑𝑠 = 𝑟𝜃+𝑠𝑟
𝑃Definition:𝑠 = 𝑟𝜃+𝑠𝑚𝜃

𝑠 = 𝑟𝜃+𝑠𝑚 =𝑚

The Total mass 𝑚 is the summationof mass 𝑚 (small particles)

𝑚 𝑟𝑔 = 0

The arm length of G is zero𝑃Equation of motion of mass point 𝑚𝑚𝑠 =𝑓internal force
𝑓 =𝑃 +𝑃𝑚𝑚𝑚

action and reaction will disappear if you calculate
the summation of internal forces𝑑𝑚𝑠𝑠𝑑𝑠 = 𝑟𝜃+𝑠𝑟
𝑃Definition:𝑠 = 𝑟𝜃+𝑠𝑚

𝜃𝑠 = 𝑟𝜃+𝑠𝑚 =𝑚

The Total mass 𝑚 is the summationof mass 𝑚 (small particles)

𝑚 𝑟𝑔 = 0

The arm length of G is zero

𝑃Equation of motion of mass point 𝑚

𝑚𝑠 =𝑓
Law ofAction andReaction

internal force
𝑓 =𝑃 +𝑃𝑟 ×𝑓 =𝑑 ×𝑃 +𝑑 ×𝑃

The summation of internal forces becomes
only external forces𝑑𝑚𝑠𝑠𝑑𝑟𝑃𝑚𝜃𝑚𝑠 =𝑚 𝑟𝜃+𝑠=

𝑚𝑟 𝜃+𝑃=𝑚 𝑠𝑓 =𝑃 +𝑃𝑚𝑠𝑑𝑚𝑠𝑠𝑑𝑟𝑃𝑚𝜃

Summation of Equation of motion

𝑚𝑠 =𝑃

Equation of motionof translation: 𝑚𝑎 = 𝐹𝑓𝑚 𝑠 =𝑃 +𝑃𝑑𝑚𝑠𝑠

𝑑𝑟𝑃𝑚𝜃𝑟 ×𝑚 𝑠 =𝑚 𝑟 𝑟𝜃+𝑠𝑃=

𝑚𝑟=𝐼𝜃𝑟 ×𝑓 =𝑑 𝑃 +𝑑 𝑃

𝜃+𝑚𝑟 𝑠𝑑𝑚𝑠𝑠𝑑𝑟

𝑃𝑚𝑃

Summation of Equation of motion × arm length𝑚 𝑟𝑠 =𝜃𝑟𝑓Equation of motionof rotation: 𝐼𝜃 = 𝑇𝑚𝑟𝜃 =𝑑 𝑃 +𝑑 𝑃

TorqueInertia
MomentAngularAcceleration

Equation of Motion of Translation𝑚𝑑𝑑 𝑚𝑑 𝑚𝑠 +𝑠 = −𝑐 (𝑠 − 𝑧 )𝑑 +𝑑𝑑 +𝑑−𝑐 (𝑠 − 𝑧 )𝑑−𝑘 (𝑠 − 𝑧 )𝜃−𝑘 (𝑠 − 𝑧 )Equation of Motion of Rotation

𝑘𝑐𝑚𝑘𝑐𝐼𝐼𝑠 −𝑠 = −𝑑 𝑐 (𝑠 − 𝑧 )𝑑 +𝑑𝑑 +𝑑+𝑑 𝑐 (𝑠 − 𝑧 )−𝑑 𝑘 (𝑠 − 𝑧 )𝑚+𝑑 𝑘 (𝑠 − 𝑧 )Equation of Motion of Unsprung-mass

𝑚 𝑧 =𝑐𝑚 𝑧 =𝑐𝑠 −𝑧+𝑘𝑠 −𝑧

𝑠 −𝑧−𝑘+𝑘𝑧 −𝑟𝑠 −𝑧−𝑘𝑧 −𝑟

𝑑 𝑚
𝑑 +𝑑𝐼
𝑑 +𝑑𝑑 𝑚𝑑 +𝑑−𝐼𝑑 +𝑑𝑠𝑠𝑧𝑧

𝑚𝑐𝑑 𝑐+ −𝑐𝑐−𝑑 𝑐−𝑐−𝑐−𝑑 𝑐𝑐−𝑐𝑑 𝑐𝑐𝑠𝑠𝑧𝑧𝑘𝑑 𝑘+−𝑘𝑘−𝑑 𝑘−𝑘−𝑑 𝑘𝑘 +𝑘−𝑘𝑚

Newmark-β method

Solving the system of equations−𝑘𝑑 𝑘𝑘+𝑘𝑠𝑠𝑧𝑧

=𝑘 𝑟
𝑘 𝑟

Consider the next time step

Newmark-β method𝒛 𝑡 + Δ𝑡 = 𝒛 𝑡 + Δ𝑡 1 − 𝛾 𝒛 𝑡 + Δ𝑡𝛾𝒛 𝑡 + Δ𝑡
1− 𝛽 𝒛 𝑡 + Δ𝑡 𝛽𝒛 𝑡 + Δ𝑡𝒛 𝑡 + Δ𝑡 = 𝒛 𝑡 + Δ𝑡𝒛 𝑡 + Δ𝑡
2𝒇 𝑡 + Δ𝑡 − 𝐂 𝒛 𝑡 + Δ𝑡 1 − 𝛾 𝒛 𝑡
1− 𝐊 𝒛 𝑡 + Δ𝑡𝒛 𝑡 + Δ𝑡−𝛽 𝒛 𝑡2

𝐀𝒛 𝑡 + Δ𝑡 = 𝒃

Once 𝒛 𝑡 , 𝒛 𝑡 , 𝒛 𝑡 and 𝒇 𝑡 + Δ𝑡 are obtained,
we can solve 𝒛 𝑡 + Δ𝑡 , 𝒛 𝑡 + Δ𝑡 and 𝒛 𝑡 + Δ𝑡 .

1) Download: http://www.kz.tsukuba.ac.jp/~yamamoto_k/material/matlab_for_ide.zip
2) Decompress the downloaded zip file3) Launch Matlab and set the decompressed directory as current directory4) Run the following script:>> sample_03

Consider the effect ofBridge Vibrations

Vehicle and Bridge can be modeled by
Half-car model(Rigid Body Spring Model)

andBeamFinite Element Methodusing 1D Finite Beam Element ModelVehicle Inputs include Bridge vibrations
as well as Road unevenness

road unevenness

bridge vibration

The Vehicle Vibrations affect on
the Contact Forces acting on the Bridge𝑑 𝑚𝑃 𝑡 =𝑔−𝑠 𝑡 +𝑚 𝑔−𝑧 𝑡
𝑑 +𝑑𝑑 𝑚=𝑔+𝑚 𝑔+𝑘 𝑧 −𝑢𝑑 +𝑑𝑑 𝑚
𝑔−𝑠 𝑡 +𝑚𝑑 +𝑑𝑑 𝑚=
𝑔 + 𝑚 𝑔 + +𝑘𝑑 +𝑑𝑃 𝑡 =weight

𝑔−𝑧 𝑡
𝑧 −𝑢

restoring forceof unsprung-stiffness

Bridge Deflection Vibrations can be calculatedby Structure Mechanics
𝑥 𝑃 𝑡𝑥

𝑃 𝑡𝑦 𝑥, 𝑡𝜌𝐴𝑦 𝑥, 𝑡 +𝜕𝜕𝐸𝐼𝑦 𝑥, 𝑡𝜕𝑥𝜕𝑥𝑥=𝛿 𝑥−𝑥 𝑡 𝑃 𝑡 +𝛿 𝑥−𝑥 𝑡 𝑃 𝑡
𝛿 𝑥 = 0 , when 𝑥 ≠ 0

Mass perunit length

FlexuralRigidity

DeltaFunction

𝛿 𝑥 = ∞ , when 𝑥 = 0
𝛿 𝑥 d𝑥 = 1

The equation of motion of Vehicle System:𝑑 𝑚𝑧
𝑑 +𝑑𝐼𝑧𝑑 +𝑑𝑑 𝑚𝑧
𝑑 +𝑑𝐼𝑡 −𝑧
𝑑 +𝑑𝑡 +𝑡 = −𝑘𝑧𝑡 −𝑧𝑡 = −𝑑 × 𝑘𝑧

𝑡−𝑘𝑧

𝑡 −𝑧

𝑡 −𝑧

𝑡+𝑑 ×𝑘𝑡𝑧−𝑐𝑧

𝑡 −𝑧

𝑡 −𝑧

𝑡−𝑑 ×𝑐𝑚𝑧

𝑡 =𝑘𝑧

𝑡 −𝑧

𝑡+𝑐𝑧

𝑡 −𝑧

𝑡−𝑘𝑧

𝑡 −𝑢 𝑡𝑚𝑧𝑡 =𝑘𝑧𝑡 −𝑧𝑡+𝑐𝑧𝑡 −𝑧𝑡−𝑘𝑧𝑡 −𝑢 𝑡𝑑𝑧𝑑𝑚𝑧𝑡−𝑐𝑧𝑧𝑡 −𝑧𝑡 −𝑧𝑡Input Profile:𝑢 𝑡 =𝑅 𝑥 𝑡𝑢 𝑡 =𝑅 𝑥 𝑡+ 𝑦 𝑥 𝑡 ,𝑡+ 𝑦 𝑥 𝑡 ,𝑡𝑧𝑢

𝑘𝑐 𝑘𝑚𝑚𝑐𝑑 𝑚𝑑 +𝑑𝑑 𝑚
𝑃 𝑡 =𝑑 +𝑑𝑧
𝑘𝑦 𝑥, 𝑡𝑢
𝐸𝐼 𝜌𝐴

The equation of motion of Bridge System:𝜌𝐴𝑦 𝑥, 𝑡 +𝜕𝜕
𝐸𝐼𝑦 𝑥, 𝑡𝜕𝑥𝜕𝑥=𝛿 𝑥−𝑥 𝑡 𝑃 𝑡 +𝛿 𝑥−𝑥 𝑡 𝑃 𝑡−𝑑 ×𝑐𝑧𝑡 −𝑧𝑡Contact Forces:
𝑃 𝑡 =𝑘𝑡𝑔−𝑧𝑡+𝑚𝑔−𝑧𝑡𝑔−𝑧𝑡+𝑚𝑔−𝑧𝑡FEM is a method for discretizing
a Continuum Body𝑥𝑦 𝑡⋮𝒚 𝑡 = 𝑦 𝑡⋮The equation of motion of Bridge can be also
expressed by matrices and vectors𝑥𝑦 𝑡
𝐌 𝒚 𝑡 +𝐊 𝒚 𝑡 =𝑳 𝑡 𝑷 𝑡
MassMatrix

Stiffness
MatrixEquivalent Nodal ForceDistribution Matrix

Contact Force

𝐊

126𝑙
𝐸𝐼 6𝑙4𝑙=𝑙 −12 −6𝑙6𝑙2𝑙−12 6𝑙−6𝑙 2𝑙
12 −6𝑙−6𝑙 4𝑙

𝐌 𝒛 𝑡 +𝐂 𝒛 𝑡 +𝐊 𝒛 𝑡 =𝐊𝑧𝑘𝑧𝑢𝑑𝑚𝑐 𝑘𝑚𝑚𝑐𝑑𝑘𝑟 𝑡𝑦 𝑡𝒓 𝑡 +𝑳 𝑡 𝒚 𝑡𝑧𝑧𝑘𝑢𝑟 𝑡𝑦 𝑡

𝐌 𝒚 𝑡 +𝐊 𝒚 𝑡 =𝑳 𝑡 𝐌 𝒈−𝒛 𝑡

𝐌
𝑳 𝑡 𝐌

𝐌𝒛
𝐂𝒚 +𝐂𝒛𝐊𝒚 +𝐊 𝑳 𝑡𝐊𝒛
𝐊 𝒓 𝑡
𝒚 = 𝑳 𝑡 𝐌 𝒈

Equation of Motion of VBI system
is NOT a linear differential equation
𝐌𝑳 𝑡 𝐌𝐌𝒛
𝐂𝒚 +𝐂𝒛𝐊𝒚 +Coefficients temporary changes

𝐊 𝑳 𝑡
𝐊𝒛
𝐊 𝒓 𝑡
𝒚 = 𝑳 𝑡 𝐌 𝒈

Newton-Raphson Methodfor solving a non-linear problem

Convergence calculation

Road ProfileVehicle
Sprung-mass:𝑚 (kg)
Front: 𝑖=1
Rear: 𝑖=2

9000

Unsprung-mass:𝑚 (kg)

Sprung-damping:𝑐 (kg/s)Sprung-stiffness:𝑘 (N/m)Unsprungstiffness:𝑘 (N/m)Distance from G:
𝑑 (m)

500

2000

4500

60000

1.75

500

2000

4500

60000

1.75

Vehicle Vibration

Input Profile

Contact Force

Bridge Vibration
BridgeFlexural Rigidity:
𝐸𝐼 (GNm2)

Mass per unit Length:𝜌𝐴 (kg/m)

Span Length:
𝐿 (m)Element Number

156

3000

100

20

Travel speed:
𝑣 (m)
10.0

Bridge

Input Profile

Road Profile

Sprung-mass vibration

Unsprung-mass vibration

1) Download: http://www.kz.tsukuba.ac.jp/~yamamoto_k/material/matlab_for_ide.zip
2) Decompress the downloaded zip file3) Launch Matlab and set the decompressed directory as current directory4) Run the following script:>> sample_04>> sample_04_application

VBI simulation is applied toBridge Design and MaintenanceTraffic Load for Bridge Design
Vibration-based Structure Health Monitoring
Future Technology (Watching Logistics)

Traffic Loadfor Bridge Design

Maximum Value of 𝑃 of 10t vehicle is about 50 (kN)
(100kN)(NOTICE: Only considering very smooth road profile)Vibration-based
Structure Health Monitoring
Bridge Acceleration Vibrations
𝐿/4

𝐿/2

3𝐿/4

Drive-by Monitoring(Simultaneous Inspection of Vehicle-Bridge-Road)

Only MeasuringVehicle Acceleration Vibrations and Position

Identifying
Vehicle Parameters, Bridge Parameters
and Road Profile

Assignment
1) Download the Matlab with Campus-wide License
2) Execute the sample programs of numerical simulations
You can change the parameters: heavy or light vehicles, rough pavements, strong or weak bridges, etc.
Discuss about the obtained results.