關于C語言的模擬物理模型

聲明：本文全部代碼效果基于C語言easyx圖形界面庫。

引言

關于很多游戲和模型的開發，都需要模擬真實的物理模型

比如：基本矢量運動模型（位移，速度，加速度），重力模型，碰撞模型，萬有引力模型。

本文我會一一介紹這幾種模型的開發。

既然是模擬顯示物理模型，那么引入的現實物理公式也是必要的。

本文所需要引用的公式有：其中TIKME_STEP是每幀變化的量，來彌補時間效果和逐幀效果的差

最后的總物理模型的話，可以選擇性給出每個球的狀態情況，比如每個方向的初速度，初加速度，以及是否受重力和球之間的萬有引力，對偏心碰撞，邊框碰撞。

代碼

位移公式： ▲x = v * t??

程序內使用：x += vx? y+= vy (逐幀使用）
void move()//移動
{for (int i = 0; i <= 20; i++) {if (a[i].live) {a[i].x += a[i].vx * TIME_STEP;a[i].y += a[i].vy * TIME_STEP;}}
}

速度公式：▲v = a * t?

程序內使用? vx += ax? vy += ay?(逐幀使用）
void speed()//速度變化
{for (int i = 0; i <= 20; i++) {if (a[i].live) {a[i].vx += a[i].ax * TIME_STEP;a[i].vy += a[i].ay * TIME_STEP;}}
}

加速度公式：a = F / m

程序內使用 (Fx+Fx1+Fx2...) / m = ax??(Fy+Fy1+Fy2...) / m = ay
void Fa()//受力
{for (int i = 0; i <= 20; i++) {if (a[i].live) {a[i].ax = (a[i].Fx + a[i].Fwx) / a[i].m;a[i].ay = (a[i].Fy + a[i].GN + a[i].Fwy) / a[i].m;}}
}

重力公式： G = mg

程序內使用 G = m*g? (g預定義數值） 然后G 是Fy的組成元素a[i].GN = a[i].m * g;

碰撞前后的動量守恒和能量守恒

m1*v1 + m2*v2 = m1*v3+m2*v4;

0.5*m1*v1*v1 + 0.5*m2*v2*v2 = 0.5*m1*v3*v3 + 0.5*m2*v3*v4;

然后再考慮是否對心碰撞分析左右偏移向量

// 檢測并處理小球碰撞
void collisionDetection()
{for (int i = 0; i < 20; i++) {if (a[i].live) {for (int j = i + 1; j <= 20; j++) {if (a[j].live) {double dx = a[i].x - a[j].x;double dy = a[i].y - a[j].y;double distance = sqrt(dx * dx + dy * dy);if (distance < 100) { // 兩個小球半徑之和為100double m1 = a[i].m;double m2 = a[j].m;double v1ix = a[i].vx;double v1iy = a[i].vy;double v2ix = a[j].vx;double v2iy = a[j].vy;// 計算碰撞法線方向的單位向量double n_x = dx / distance;double n_y = dy / distance;// 計算相對速度在法線方向和切線方向的分量double v1n = v1ix * n_x + v1iy * n_y;double v2n = v2ix * n_x + v2iy * n_y;double v1t = v1iy * n_x - v1ix * n_y;double v2t = v2iy * n_x - v2ix * n_y;// 應用動量守恒和恢復系數計算碰撞后的法向速度double v1nf = ((m1 - m2) * v1n + 2 * m2 * v2n + RESTITUTION_COEFFICIENT * m2 * (v2n - v1n)) / (m1 + m2);double v2nf = ((m2 - m1) * v2n + 2 * m1 * v1n + RESTITUTION_COEFFICIENT * m1 * (v1n - v2n)) / (m1 + m2);// 切線方向速度不變（假設無摩擦力）double v1tf = v1t;double v2tf = v2t;// 將法向和切向速度轉換回笛卡爾坐標系a[i].vx = v1nf * n_x - v1tf * n_y;a[i].vy = v1nf * n_y + v1tf * n_x;a[j].vx = v2nf * n_x - v2tf * n_y;a[j].vy = v2nf * n_y + v2tf * n_x;}}}}}
}

邊框碰撞模型

void peng()//碰撞影響一定放到位置更新之后
{for (int i = 0; i <= 20; i++) {if (a[i].live) {if (a[i].y - 50 < 0) {a[i].y = 50;  // 精確調整位置到邊界a[i].vy = -a[i].vy;  // 確保速度反向}if (a[i].x - 50 < 0) {a[i].x = 50;a[i].vx = -a[i].vx;}if (a[i].y + 50 > high) {a[i].y = high - 50;a[i].vy = -a[i].vy;}if (a[i].x + 50 > width) {a[i].x = width - 50;a[i].vx = -a[i].vx;}}}
}

檢測遇到邊界后，改變水平或者垂直方向的速度方向相反，保持速度大小不變，

實現碰撞墻壁效果。

萬有引力模型

F萬 =??G*m1*m2/(x*x)


void F_allocation()
{double sin_w = (a[0].y - a[1].y) / sqrt((a[0].y - a[1].y) * (a[0].y - a[1].y) + (a[0].x - a[1].x) * (a[0].x - a[1].x));double cos_w = (a[0].x - a[1].x) / sqrt((a[0].y - a[1].y) * (a[0].y - a[1].y) + (a[0].x - a[1].x) * (a[0].x - a[1].x));double Fw = Gk * a[0].m * a[1].m / ((a[0].y - a[1].y) * (a[0].y - a[1].y) + (a[0].x - a[1].x) * (a[0].x - a[1].x));a[0].Fwx = -Fw * cos_w;a[0].Fwy = -Fw * sin_w;a[1].Fwx = -a[0].Fwx;a[1].Fwy = -a[0].Fwy;
}

關于運行前的準備

預編譯和結構體

#include <graphics.h>
#include <conio.h>
#include <stdio.h>
#include <time.h>
#include <stdbool.h>
#include <math.h>#define g 98
#define width 1700
#define high 1000
#define TIME_STEP 0.05  // 固定時間步長
#define Gk 677770
#define RESTITUTION_COEFFICIENT 0.8 // 恢復系數，0 <= e <= 1，0為完全非彈性，1為完全彈性typedef struct stu {double x;double y;double ax;double ay;double vx;double vy;double Fx;double Fy;double Fwx;double Fwy;double m;double GN;bool live;
} G;

代碼是集成的，每個作用都在，如果要研究其中某種模型，只需要調用該模型函數即可。

比如研究萬有引力模型，需要將重力模型關掉，然后控制不要碰撞，這里數值我測出來了一個效果就是，用假設初速度完全滿足萬有引力提供向心力從而做圓周運動的模型，

我們需要將兩個球的質量全設置為10，然后鎖定右側球，兩球高度一致，賦予左球vy = 130便可以滿足萬有引力完全提供向心力的天體勻速圓周模型。

展示總函數

#include <graphics.h>
#include <conio.h>
#include <stdio.h>
#include <time.h>
#include <stdbool.h>
#include <math.h>#define g 98
#define width 1700
#define high 1000
#define TIME_STEP 0.05  // 固定時間步長
#define Gk 677770
#define RESTITUTION_COEFFICIENT 0.8 // 恢復系數，0 <= e <= 1，0為完全非彈性，1為完全彈性typedef struct stu {double x;double y;double ax;double ay;double vx;double vy;double Fx;double Fy;double Fwx;double Fwy;double m;double GN;bool live;
} G;G a[21];void Fa()//受力
{for (int i = 0; i <= 20; i++) {if (a[i].live) {a[i].ax = (a[i].Fx + a[i].Fwx) / a[i].m;a[i].ay = (a[i].Fy + a[i].GN + a[i].Fwy) / a[i].m;}}
}void speed()//速度變化
{for (int i = 0; i <= 20; i++) {if (a[i].live) {a[i].vx += a[i].ax * TIME_STEP;a[i].vy += a[i].ay * TIME_STEP;}}
}void move()//移動
{for (int i = 0; i <= 20; i++) {if (a[i].live) {a[i].x += a[i].vx * TIME_STEP;a[i].y += a[i].vy * TIME_STEP;}}
}void draw()
{for (int i = 0; i <= 20; i++) {if (a[i].live) {solidcircle((int)a[i].x, (int)a[i].y, 50);}}
}void peng()//碰撞影響一定放到位置更新之后
{for (int i = 0; i <= 20; i++) {if (a[i].live) {if (a[i].y - 50 < 0) {a[i].y = 50;  // 精確調整位置到邊界a[i].vy = -a[i].vy;  // 確保速度反向}if (a[i].x - 50 < 0) {a[i].x = 50;a[i].vx = -a[i].vx;}if (a[i].y + 50 > high) {a[i].y = high - 50;a[i].vy = -a[i].vy;}if (a[i].x + 50 > width) {a[i].x = width - 50;a[i].vx = -a[i].vx;}}}
}void init()
{for (int i = 0; i <= 20; i++) {a[i].x = 0;a[i].y = 0;a[i].ax = 0;a[i].ay = 0;a[i].vx = 0;a[i].vy = 0;a[i].Fx = 0;a[i].Fy = 0;a[i].Fwx = 0;a[i].Fwy = 0;a[i].m = 10;a[i].live = false;a[i].GN = a[i].m * g;}
}void F_allocation()
{double sin_w = (a[0].y - a[1].y) / sqrt((a[0].y - a[1].y) * (a[0].y - a[1].y) + (a[0].x - a[1].x) * (a[0].x - a[1].x));double cos_w = (a[0].x - a[1].x) / sqrt((a[0].y - a[1].y) * (a[0].y - a[1].y) + (a[0].x - a[1].x) * (a[0].x - a[1].x));double Fw = Gk * a[0].m * a[1].m / ((a[0].y - a[1].y) * (a[0].y - a[1].y) + (a[0].x - a[1].x) * (a[0].x - a[1].x));a[0].Fwx = -Fw * cos_w;a[0].Fwy = -Fw * sin_w;//a[1].Fwx = -a[0].Fwx;//a[1].Fwy = -a[0].Fwy;
}// 檢測并處理小球碰撞
void collisionDetection()
{for (int i = 0; i < 20; i++) {if (a[i].live) {for (int j = i + 1; j <= 20; j++) {if (a[j].live) {double dx = a[i].x - a[j].x;double dy = a[i].y - a[j].y;double distance = sqrt(dx * dx + dy * dy);if (distance < 100) { // 兩個小球半徑之和為100double m1 = a[i].m;double m2 = a[j].m;double v1ix = a[i].vx;double v1iy = a[i].vy;double v2ix = a[j].vx;double v2iy = a[j].vy;// 計算碰撞法線方向的單位向量double n_x = dx / distance;double n_y = dy / distance;// 計算相對速度在法線方向和切線方向的分量double v1n = v1ix * n_x + v1iy * n_y;double v2n = v2ix * n_x + v2iy * n_y;double v1t = v1iy * n_x - v1ix * n_y;double v2t = v2iy * n_x - v2ix * n_y;// 應用動量守恒和恢復系數計算碰撞后的法向速度double v1nf = ((m1 - m2) * v1n + 2 * m2 * v2n + RESTITUTION_COEFFICIENT * m2 * (v2n - v1n)) / (m1 + m2);double v2nf = ((m2 - m1) * v2n + 2 * m1 * v1n + RESTITUTION_COEFFICIENT * m1 * (v1n - v2n)) / (m1 + m2);// 切線方向速度不變（假設無摩擦力）double v1tf = v1t;double v2tf = v2t;// 將法向和切向速度轉換回笛卡爾坐標系a[i].vx = v1nf * n_x - v1tf * n_y;a[i].vy = v1nf * n_y + v1tf * n_x;a[j].vx = v2nf * n_x - v2tf * n_y;a[j].vy = v2nf * n_y + v2tf * n_x;}}}}}
}int main()
{init();a[0].live = true;a[0].x = 200;a[0].y = 500;/*a[0].vx = 50;a[0].vy = 20;*/a[1].m = 10;a[1].live = true;a[1].x = 600;a[1].y = 500;/*a[1].vx = 30;a[1].vy = 20;*/// 初始化圖形窗口，大小為 width x highinitgraph(width, high);// 設置填充顏色為綠色setfillcolor(GREEN);BeginBatchDraw();while (!_kbhit()) {cleardevice(); // 清除繪圖// F_allocation();Fa();speed();move();peng();collisionDetection(); // 檢測并處理碰撞draw();FlushBatchDraw();Sleep(8);}EndBatchDraw();// 關閉圖形窗口closegraph();return 0;
}

本文展現的是集成的物理模型框架，感興趣的可以嘗試調試實現各種效果。

我自己測試的是

1，萬有引力完全提供向心力的勻速圓周運動

2，對偏心碰撞的守恒

3，以及編寫初期的重力測試，和碰撞測試，加速度測試。

至于運行效果，這博客也不好發視頻，大家可以嘗試敲一下代碼運行試試。

當物理模型自己結合物理公式和編程寫出來展現在眼前的物理運動時，是真的相信課本里說的天體運動，也是切身感受到了，還有碰撞問題，感覺可以寫個桌球小游戲，不過這些交給大家發揮了，本文到此，感謝觀看。