fix: trapped particle flux

This commit is contained in:
liuyihui 2022-05-14 21:51:24 +08:00
parent a0fea724d1
commit 16e6a77aea
4 changed files with 97 additions and 7 deletions

View File

@ -137,7 +137,7 @@ $$
<div align=center><img src="docs/ExFit-trapped-electron.png" style="max-width: 80%;"><br /></div>
2. 质子
与电子保持一致,额外考虑幂律平移模型$j=j_0(E-\beta)^{-\alpha}$同电子进行`BIC test`,最终选用幂律平移模型$j_0=559.76377,\ \alpha=2.40873,\ \beta=-1.53424$。
与电子保持一致,幂律模型修改为$j=j_0(E-\beta)^{-\alpha}$,进行`BIC test`,最终选用幂律模型$j_0=559.76377,\ \alpha=2.40873,\ \beta=-1.53424$。
<div align=center>
@ -151,19 +151,19 @@ $$
<div align=center><img src="docs/ExFit-trapped-proton.png" style="max-width: 80%;"><br /></div>
## 辐射源设置
使用`Geant4`的`General Particle Source``GPS`)生成源,考虑在`World`球面上均匀随机发射粒子(各向同性),保证入射方向为球内,则球面上的总通量为:
使用`Geant4`的`General Particle Source``GPS`)生成源,考虑在`World`球面上均匀随机发射粒子(各向同性),保证入射方向为球内,则球面上的计数率为:
$$
F = A\int_{E_a}^{E_b} \phi(E)\mathrm{d}E,\ A=4\pi R^2
$$
`GCR`的单位为$\mathrm{m^{-2}sr^{-1}s^{-1}}$,乘以$4\pi A$即可得到通量;俘获辐射的单位为$\mathrm{cm^{-2}s^{-1}}$,乘以$A$即可得到通量,最终结果如下表所述。
`GCR`的单位为$\mathrm{m^{-2}sr^{-1}s^{-1}}$,乘以$4\pi A$即可得到计数率;俘获辐射的单位为$\mathrm{cm^{-2}s^{-1}}$,乘以$A$即可得到计数率。经过检验模型中使用的外壳3层2mm 铝 + 10mm 芳纶 + 5mm 铝)能够抵挡能量在$50\ \mathrm{MeV}$以下的质子和能量在$100\ \mathrm{keV}$以下的电子,最终结果如下表所述。
<div align=center>
| Partical | Model | $E_{min}(\mathrm{MeV})$ | $E_{max}(\mathrm{MeV})$ | Flux($\mathrm{s^{-1}}$) |
|:--------:|:--------------------------------------------------------------------:|:-----------------------:|:-----------------------:|:-----------------------:|
| $e^{-}$ | $j=j_0e^{-\frac{E}{E_0}}$ | 0.04 | 10 | 1.91227E+12 |
| `p` | $j=j_0(E-\beta)^{-\alpha}$ | 0.1 | 1000 | 5623595363 |
| $e^{-}$ | $j=j_0e^{-\frac{E}{E_0}}$ | 0.1 | 10 | 9.23242E+11 |
| `p` | $j=j_0(E-\beta)^{-\alpha}$ | 50 | 1000 | 42853994.15 |
| `H` | $\frac{C_1\beta^{\alpha-1}}{pc}\left(\frac{pc}{pc+C_2}\right)^{C_3}$ | 200 | 100000 | 12854310.92 |
| `He` | $C_1e^{-C_2E}(1-e^{-C_3E+C_4})$ | 200 | 100000 | 1718828.12 |
| `Li` | $C_1e^{-C_2E}(1-e^{-C_3E+C_5})$ | 200 | 100000 | 12982.36284 |
@ -195,6 +195,41 @@ $$
</div>
考虑到`GCR`和俘获辐射计数率相差太大,分为两部分分别进行模拟,同时按比例进行放缩。
### GCR
电子计数率与质子计数率之比约为`21544`,设置质子的源强为`1`,模拟`923285`个入射事件(数值上为$F\times10^{-6}\ \mathrm{s}$),因此最终的剂量需要乘以$10^6\ \mathrm{s^{-1}}\times30\ \mathrm{d}\times86400\ \mathrm{s/d}$倍。
```macro
/gps/source/add 21544
/gps/particle e-
/gps/pos/type Surface
/gps/pos/shape Sphere
/gps/pos/radius 15 m
/gps/ene/type Exp
/gps/ene/min 100 keV
/gps/ene/max 10 MeV
/gps/ene/ezero 0.0824
/gps/ang/type iso
/gps/ang/maxtheta 90 deg
/gps/source/add 1
/gps/particle proton
/gps/pos/type Surface
/gps/pos/shape Sphere
/gps/pos/radius 15 m
/gps/ene/type Pow
/gps/ene/min 200 MeV
/gps/ene/max 100 GeV
/gps/ene/alpha 2.40873
/gps/ang/type iso
/gps/ang/maxtheta 90 deg
/run/beamOn 923285
```
### 俘获辐射
假定`Co`的计数率为1设定其他粒子的计数率模拟`707013`个入射事件,因此最终的剂量需要乘以$20.9\ \mathrm{s^{-1}}\times30\ \mathrm{d}\times86400\ \mathrm{s/d}$倍。
## 空间站建模
1. 尺寸与分区[^6]
<div align=center><img src="docs/size.webp" style="max-width: 50%;"></div>

28
debug.mac Normal file
View File

@ -0,0 +1,28 @@
# 多线程设置
/run/numberOfThreads 1
# verbose
/control/verbose 1
/run/verbose 0
/event/verbose 0
/tracking/verbose 0
/gps/verbose 0
# 初始化
/run/initialize
# 电子射程
/gps/particle e-
/gps/energy 100 keV
/gps/position 0 0 5 m
/gps/ang/type iso
/gps/ang/maxphi 60 deg
/gps/ang/maxtheta 20 deg
# 质子射程
/gps/particle proton
/gps/energy 50 MeV
/gps/position 0 0 5 m
/gps/ang/type iso
/gps/ang/maxphi 60 deg
/gps/ang/maxtheta 20 deg

View File

@ -13,7 +13,7 @@
/gps/source/clear
# Trapped electron
/gps/source/add 1
/gps/source/add 21544
/gps/particle e-
/gps/pos/type Surface
/gps/pos/shape Sphere
@ -25,5 +25,19 @@
/gps/ang/type iso
/gps/ang/maxtheta 90 deg
/gps/source/add 1
/gps/particle proton
/gps/pos/type Surface
/gps/pos/shape Sphere
/gps/pos/radius 15 m
/gps/ene/type Pow
/gps/ene/min 200 MeV
/gps/ene/max 100 GeV
/gps/ene/alpha 2.40873
/gps/ang/type iso
/gps/ang/maxtheta 90 deg
/run/beamOn 923285
# list
/gps/source/list
/gps/source/list

View File

@ -229,6 +229,19 @@ G4VPhysicalVolume* DetectorConstruction::Construct() {
G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld, G4Material::GetMaterial("Vacuum"), "World");
G4VPhysicalVolume* physicsWorld = new G4PVPlacement(0, G4ThreeVector(), logicWorld, "World", 0, false, 0, true);
// Debug for particle range
// G4Box* solidPlane1 = new G4Box("Plane1", 10 * m, 10 * m, 2 * mm);
// G4LogicalVolume* logicPlane1 = new G4LogicalVolume(solidPlane1, G4Material::GetMaterial("AluminumAlloySeries5"), "Plane1");
// new G4PVPlacement(0, G4ThreeVector(), logicPlane1, "Plane1", logicWorld, false, 0, true);
// G4Box* solidPlane2 = new G4Box("Plane2", 10 * m, 10 * m, 10 * mm);
// G4LogicalVolume* logicPlane2 = new G4LogicalVolume(solidPlane2, G4Material::GetMaterial("Taparan"), "Plane2");
// new G4PVPlacement(0, G4ThreeVector(0, 0, -6), logicPlane2, "Plane2", logicWorld, false, 0, true);
// G4Box* solidPlane3 = new G4Box("Plane3", 10 * m, 10 * m, 5 * mm);
// G4LogicalVolume* logicPlane3 = new G4LogicalVolume(solidPlane3, G4Material::GetMaterial("AluminumAlloySeries5"), "Plane3");
// new G4PVPlacement(0, G4ThreeVector(0, 0, -14.5), logicPlane3, "Plane3", logicWorld, false, 0, true);
// 节点舱
ConstructSectionSphere(logicWorld, (5.18 + 0.33) * m);
// 过渡段 1