🚙

空気抵抗と加速帯のシミュレーション

2023/10/23に公開


青い部分に入ると抗力を受けて減速し、赤い部分では加速させられる

公式の簡略化

実際の複雑な公式から「密度」や「面積」などを省くと単に「速度の二乗」×「抗力係数」になる。そこで時速 50 Km で抗力係数を -0.001 とした場合、

50 * 50 * -0.001  # => -2.5

となり時速 2.5 Km の空気抵抗を受けているとわかる。

遅いほど空気抵抗が減る

50 * 50 * -0.001  # => -2.5
40 * 40 * -0.001  # => -1.6
30 * 30 * -0.001  # => -0.9

減速するほど空気抵抗が減っていくのがわかる。

反対方向に移動しない

ここがおもしろいところで、普通は加速度に反対の力を加えていくとそのうち逆方向に進むが、空気抵抗の場合その力は自分自身の速度から生まれているので、逆方向に進むことはない。例えば速度 0 Km まで減速していくと加速度に加える力は 0 になる。

2 * 2 * -0.001  # => -0.004
1 * 1 * -0.001  # => -0.001
0 * 0 * -0.001  # => -0.0

わかりやすく言えば、自転車で停止していると風力を感じなくなるのと同じ。

流体抵抗は抗力が強くなっただけ

仮に空気抵抗の抗力係数を -0.001 から -0.01 に変更してみると、

50 * 50 * -0.01  # => -25.0

時速 50 Km で走っていた車が一瞬で 25 Km 減速した。これはもう海に落ちたとしか思えない。

抗力係数に符号を含める意味

Nature of Code では抗力係数に符号をつけていないが、抗力係数自体に符号を持たせた方が自由度が広がる。例えば、抗力係数の符号をプラスにしただけで、

50 * 50 * 0.001  # => 2.5

50 Km で走っている車はさらに加速していく。なので不思議な力でブーストさせたい場合などに特別なロジックはいらない。抗力係数の符号を反転させるだけでよい。

抗力係数を表す変数

Nature of Code では変数 c が使われているが、これじゃなんのことかわからない。c はおそらく抗力(drag)の係数(coefficient)の c だと思われるが単語にしたところで coefficient になるので何の係数かわからない。したがって drag_coefficient とそのまま書けばいいと思う。

計算式まとめ

物体の属性
location = V.zero
velocity = V[0.0, 50.0]
acceleration = V.zero
物体への空気抵抗
drag_coefficient = -0.001
speed = velocity.magnitude                           # => 50.0
drag_magnitude = speed * speed * drag_coefficient    # => -2.5
force = velocity.normalize_or_zero * drag_magnitude  # => (-0.0, -2.5)
物体の加速度に加算して動かす
acceleration += force     # => (0.0, -2.5)
velocity += acceleration  # => (0.0, 47.5)
location += velocity      # => (0.0, 47.5)
acceleration *= 0         # => (0.0, -0.0)
最小限のコード
class App < Base
  def initialize
    super

    a = window_wh * V.zero
    b = window_wh * V.one

    @location = a
    @velocity = (b - a).normalize * 40.0
  end

  def draw
    super

    drag_coefficient = -0.01
    drag_magnitude = @velocity.length**2 * drag_coefficient
    force = @velocity.normalize_or_zero * drag_magnitude

    @velocity += force
    @location.replace(@location + @velocity)

    vputs "force: #{force.round(8)}"
    vputs "velocity: #{@velocity.round(8)}"

    circle_draw(@location, 30)
  end

  show
end
実験用のコード
class LiquidArea
  attr_accessor :location

  def initialize
    @location = $app.window_wh * V[0.5, 0.3]
    @rect = $app.window_wh * V[1.0, 0.2]
    @drag_coefficient = -0.02
    @color = :blue_light
  end

  def update
    $app.balls.each do |ball|
      if contain?(ball)
        speed = ball.velocity.magnitude
        drag_magnitude = speed**2 * @drag_coefficient
        force = ball.velocity.normalize_or_zero * drag_magnitude
        ball.acceleration_add(force)
      end
    end
  end

  def contain?(target)
    s = @location
    t = target.location
    (s.x - @rect.x / 2) <= t.x && t.x <= (s.x + @rect.x / 2) && (s.y - @rect.y / 2) <= t.y && t.y <= (s.y + @rect.y / 2)
  end

  def draw
    $app.rect_draw(@location, @rect, color: @color)
  end
end

class AccelerationArea < LiquidArea
  def initialize
    super

    @location = $app.window_wh * V[0.5, 0.7]
    @drag_coefficient = 0.02
    @color = :red_light
  end
end

class Ball
  attr_accessor :location
  attr_accessor :velocity

  def initialize
    @location = V.rand * $app.window_wh
    @velocity = V[rand(-0.5..0.5), rand(0.0..1.0)]
    @acceleration = V.zero
    @mass = rand(1.0..2.0)
    @gravity = V[0.0, 0.005]
  end

  def update
    acceleration_add(@gravity * @mass)
    @velocity += @acceleration
    @velocity = @velocity.clamp_length(1.0, 20.0) if false
    @location.replace(@location + @velocity)
    @location.replace(@location.modulo($app.window_wh))
    @acceleration = V.zero
  end

  def draw
    $app.point_draw(@location, radius: @mass * 2, color: :grey)
  end

  def acceleration_add(force)
    @acceleration += force / @mass
  end
end

class App < Base
  attr_accessor :balls

  def initialize
    super

    reset
  end

  def reset
    @objects = []

    @balls = 100.times.collect { Ball.new }
    @objects.concat(@balls)

    @objects << LiquidArea.new
    @objects << AccelerationArea.new

    @points = @objects.collect(&:location)
  end

  def button_down(id)
    super

    case id
    when Gosu::KB_R
      reset
    end
  end

  def draw
    super

    @objects.each(&:draw)
    @objects.each(&:update)
  end

  show
end

参照

Discussion