有限元求解器¶
单元类型¶
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}_{n+1}\right]\]
速度更新¶
\[\mathbf{v}_{n+1} = \mathbf{v}_n + \Delta t \left[(1-\gamma)\mathbf{a}_n + \gamma \mathbf{a}_{n+1}\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\) 为阻尼系数。