有限元求解器¶
单元类型¶
PhysEngine支持以下单元类型：
| 单元类型 | 描述 | 自由度 |
|---|---|---|
| FEMS4R | 4节点壳单元 | 位移 + 旋转 |
| SPH2D4N | 4节点2D实体单元 | 双线性位移 |
| SPH3D8N | 8节点3D实体单元 | 三线性位移 |
形状函数¶
四边形单元（2D）¶
自然坐标系 $(\xi, \eta) \in [-1, 1]^2$
$$N_i(\xi, \eta) = \frac{1}{4}(1 + \xi_i \xi)(1 + \eta_i \eta)$$
节点编号与自然坐标：
| 节点 | $\xi_i$ | $\eta_i$ |
|---|---|---|
| 1 | -1 | -1 |
| 2 | +1 | -1 |
| 3 | +1 | +1 |
| 4 | -1 | +1 |
六面体单元（3D）¶
自然坐标系 $(\xi, \eta, \zeta) \in [-1, 1]^3$
$$N_i(\xi, \eta, \zeta) = \frac{1}{8}(1 + \xi_i \xi)(1 + \eta_i \eta)(1 + \zeta_i \zeta)$$
数值积分¶
高斯积分规则¶
2D四边形单元（2×2 Gauss点）¶
| Gauss点 | $\xi$ | $\eta$ | 权重 $w_i$ |
|---|---|---|---|
| 1 | -1/√3 | -1/√3 | 1 |
| 2 | +1/√3 | -1/√3 | 1 |
| 3 | -1/√3 | +1/√3 | 1 |
| 4 | +1/√3 | +1/√3 | 1 |
3D六面体单元（2×2×2 Gauss点）¶
8个Gauss点，权重均为1。
Newmark时间积分¶
结构动力学方程：
$$\mathbf{M}\ddot{\mathbf{u}} + \mathbf{C}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{F}$$
采用Newmark-$\beta$方法离散：
位移更新¶
$$\mathbf{d}{n+1} = \mathbf{d}_n + \Delta t \mathbf{v}_n + \frac{\Delta t^2}{2}\left[(1-2\beta)\mathbf{a}_n + 2\beta \mathbf{a}\right]$$
速度更新¶
$$\mathbf{v}{n+1} = \mathbf{v}_n + \Delta t \left[(1-\gamma)\mathbf{a}_n + \gamma \mathbf{a}\right]$$
参数选择¶
| 方法 | $\beta$ | $\gamma$ | 特性 |
|---|---|---|---|
| 平均加速度 | 1/4 | 1/2 | 无条件稳定 |
| 线性加速度 | 1/6 | 1/2 | 条件稳定 |
| Houbolt | 9/49 | 3/2 | 二阶精度 |
PhysEngine默认采用平均加速度法（$\beta = 0.25$，$\gamma = 0.5$），保证无条件稳定。
沙漏控制¶
在缩减积分单元中，为抑制零能模式（沙漏模式），引入沙漏粘性阻尼：
$$\mathbf{f}_{hourglass} = \alpha_h \mathbf{H} \mathbf{v}$$
其中 $\mathbf{H}$ 为沙漏模式矩阵，$\alpha_h$ 为阻尼系数。