力学的システムの動的シミュレーション
Dynamic Simulation for Mechanical System

車両-橋梁相互作用 01
Vehicle-Bridge Interaction

システム情報系 助教

山本亨輔 Kyosuke
YAMAMOTO

yamamoto_k@kz.tsukuba.ac.jp

作成：2018年6月5日

Today’s Goal

Numerical Solutions of Equation of Motion（運動方程式の数値解）

 Equation of Motion:




External Force Vector ( ∈ ): input, usually given
Mass Matrix ( ∈): system parameter
Displacement Vector ( ∈ ): output, or solution
Considering a mechanical “system”, force and displacement fields are input and output of the system, respectively
A matrix acting on the output is system parameters: e.g.) Mass Matrix and Stiffness Matrix ()is DOF (degree of freedom, 自由度): “point” “direction”
is the “solution“ of the equation of motion as “differential equation”
The displacement field
So-called “Simulation” is a method to calculate solutions numerically (numerical solution)

2

Skill you can learn

Model, Solve and Visualize It!
（モデル化し、解いて、視覚化する）

 To simulate the motions ofthe mass-points(Basic Elements of System Parameters): mass: damping: spring stiffness3

SDOF: Single Degree of Freedom

 Equation of Motion:sum ofexternal force

inertiaterm

 Force is an agent that changes theshape and motion of an objectmotion is velocity, and the change of motion is acceleration

 Load: becomesacceleration = the change of motion:
displacement = the change of shape:

4Special Conditions of Equation of Motion:

 Equilibrium:0 The sum of external force equals to zero:0 The acceleration (the change of motion) becomes zero from Inertia: “object keeps its state of motion” = “no force, same velocity”

 Free Vibration:0 If there is no Load: , the system is in the state of free vibration
homogeneous The equation of motion becomes LHDE: Linear
Differential Equation

05

Analysis Analysis is “Mathematics to calculate exact solution”
 Free Vibration:
0Assuming that

sin
cos
sin
and

Substitution of

into the given equation:

0000

AngularFrequency

0 Not Requested

sin
Solution
6

Time Differential for Numerical Simulation

Taylor Expansion for approximate formula（テイラー展開による近似式の導入）

 Discretization of time Differential∆∆

Central Differential2⋅∆
 Taylor Expansion
∆∆

∆

∆

∆

∆

21212∆∆∆

161626∆⋅⋅⋅∆⋅⋅⋅∆⋅⋅⋅∆∆

Forward Scheme∆∆

Backstep Scheme

7

Example 01:

Numerical Simulation of Free Vibration for SDOF by Forward Differential
 Governing Equation:
 Discretization:

 Substitution:

1
∆
1
∆
1
∆

2∆

∆
2∆

22

∆

∆

 Discretized Equation:
2∆

2

∆
known

∆

1known Matlab Codingx(tt,:) = 2*x(tt-1,:)-(k/m*dt^2+1)*x(tt-2,:);8

Example 01:

Numerical Simulation of Free Vibration for SDOF by Forward Differential

 Function File:function [ t, x ] = ex_01_sdof( m, k, X, dt ,T )%-- preparationt = (0:dt:T)’; %-- time
TT = length(t); %-- length of solution

time:1∆2%-- initializationx = zeros(TT,size(X,2));x([1 2],:) = X([1 2],:);

2∆

3∆

InitialCondition

10⋮⋮TT%-- upwind scheme
for tt = 3:TT
x(tt,:) = 2*x(tt-1,:)-(k/m*dt^2+1)*x(tt-2,:);end

solution:∆∆

The length oftime-series data

end

9

Example 01:

Numerical Simulation of Free Vibration for SDOF by Forward Differential Launcher Script:
1.01.0∆00.0001

11.01.0[ t, x ] = ex_01_sdof( m, k,X,dt , T )[ t, x ] = ex_01_sdof( 1, 1, [0;0.0001], 0.0001 ,100 );figure(1); plot(t,x)
set(gca,'XLim', [0 100], 'XGrid', 'on')set(gca,'YLim', [-1.5 1.5], 'YGrid', 'on')set(gca,'FontSize',18)
print('ex_01_sdof_01','-dmeta','-r1200') print out

∆

00.01

1010∆visualization[ t, x ] = ex_01_sdof( 1, 1, [0;0.01], 0.01 ,100 );figure(2); plot(t,x)
set(gca,'XLim', [0 100], 'XGrid', 'on')set(gca,'YLim', [-1.5 1.5], 'YGrid', 'on')set(gca,'FontSize',18)
print('ex_01_sdof_02','-dmeta','-r1200')
10

Example 01:

Numerical Simulation of Free Vibration for SDOF by Forward Differential

 Results and Discussion
1.01.0∆00.0001

11.01.0

∆

00.01

1

 Larger ∆ is Low Computation Cost and Low Accuracy
 For not enough small ∆ , the explosion of solution occurs

11Time Differential for Numerical Simulation

Newmark- scheme（ニューマークβ法による高精度な計算）

 Discretization
1) The unknowns∆∆∆

are replaced within the governing equation

2) All factors consisting the obtained equation are known, except∆∆∆∆ ⋅1

∆
∆∆

12

∆

12

SDOF: Single Degree of Freedom

 Equation of Motion:
sum ofexternal force

inertiaterm

 The Governing Equation:

 MCK-system

13

Example 02:

Numerical Simulation of SDOF by Newmark-β method Governing Equation:
 Discretization:

∆∆∆

∆ ⋅1

∆∆

12

∆

 Substitution:
∆∆

∆ 1

∆

∆12

∆∆∆

∆

 Discretized Equation:∆ 1∆12∆∆∆∆∆∆known14

MDOF: Multi-Degree of Freedom

 Equation of Motion:

15

MDOF: Multi-Degree of Freedom
 Governing Equation:
 Discretization:

∆∆

∆

∆ ⋅1

∆
∆

12

∆

 Substitution:
∆∆

∆ 1∆

∆

12

∆∆

∆

∆

∆

12

∆

∆

∆

∆

 Discretized Equation:∆GlobalMatrix∆∆ 1Solution RHS

Right-Hand Side

The problem ofInverse Matrix
16Newmark-β method0 Function File∆00
the inputs

function [ t, x, dx, ddx ] = newmark_beta( m,c,k,f,dt,X, G, B )%% initialization
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);∆∆ 1∆
t=(0:T-1)*dt;%% solving
1
A = m + dt*G*c + dt^2*B*k; %-- global matrix
∆
∆
∆
∆
∆
2
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%-- 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

17

Example 03:

Numerical Simulation of MDOF by Newmark-β methodGoverning Equation:

010.05

0.01

0.01

-0.5

0.5

System Parameters:

1.0


-1

2.0

1.0

4.0

8.0

Input and Initial Condition:

0018

Example 03:

Numerical Simulation of MDOF by Newmark-β method Launcher Script:%-- System Parameters
m1 = 1.00; m2 = 2.00;
c1 = 0.05; c2 = 0.01; c12= 0.01;
1.0
2.0
0.05 0.01 0.01
k1 = 1.00; k2 = 4.00; k12= 8.00;%-- Simulation Parameters∆G = 1/2; B = 1/6;
1/2
1/60.01100dt = 0.01; T = 100;t = 0:dt:T; TT = length(t);%-- Input
f = zeros(2,TT);%-- Initial Condition
000-0.50X = [-0.5 0 0;0000.500.5 0 0];%-- Making Matrices
M = [ m10;0m2];
C = [ c1+c12 -c12;-c12c1+c12];
K = [ k1+k12 -k12;-k12k2+k12];%-- Simulation[ t, x, dx, ddx ] = newmark_beta( M,C,K,f,dt,X, G, B );1.0

4.08.00019

Example 03:

Numerical Simulation of MDOF by Newmark-β method

 Visualization%-- Visualization
figure(1);plot(t, [1;-1]*ones(1,TT),'k:')hold on;plot(t,[-1+x(1,:); 1+x(2,:)]);
hold offset(gca,'XLim', [0 100], 'XGrid', 'on')set(gca,'YLim', [-2 2], 'YGrid', 'on')set(gca,'FontSize',18)
print('ex_03_nmb_01','-dmeta','-r1200')%-- Animation
figure(2)writerObj = VideoWriter('ex_03_nmb_sample.avi','Motion JPEG AVI');
writerObj.FrameRate = 20;
open(writerObj);set(gcf, 'Unit', 'pixel', 'Position',[240, 120, 600, 600]);
fs=100; %-- flame sizefor tt=1:TT
figure(2);subplot('Position',[0.1 0.1 0.2 0.8])
BALL_PLOT=plot(0,[-1+x(1,tt); 1+x(2,tt)],'o');set(BALL_PLOT(1), 'FaceColor', [0,0,1])
set(BALL_PLOT(2), 'FaceColor', [1,0,0])ylim([-2 2])
subplot('Position',[0.4 0.1 0.5 0.8])
GRID_PLOT=plot(t, [-1;1]*ones(1,TT),'k:');hold on
CURVE_PLOT=plot(t(1:tt),[-1+x(1,tt); 1+x(2,tt)],'-');
hold offset(CURVE_PLOT(1), 'Color', [0,0,1])
set(CURVE_PLOT(2), 'Color', [1,0,0])xlim([0 T])
frame = getframe(figure(2));
writeVideo(writerObj,frame);
endclose(writerObj); %-- end of animation

20

Example 03:

Numerical Simulation of MDOF by Newmark-β method21

22