DOE2

这是一段关于DOE-2建筑能耗模拟软件中太阳能反射计算的技术说明。它描述了三种处理外墙、门窗等表面太阳能增益的计算模式。

核心概念¶

反射因子 (Reflection Factors)：

• RF-AWSUN、RF-AWSKY等是预计算的数据
• 含义：表面获得的太阳辐射增益 ÷ 水平面的太阳辐照度
• 应用场景：外墙面、门、玻璃等围护结构
• 辐射来源：太阳直射和天空漫射

控制参数 ISREFL：

• 由输入文件中的 SOLAR-REFL-CALC 关键字设定
• 决定反射因子的使用方式

三种计算模式¶

1. SOLAR-REFL-CALC=NO-CALC (ISREFL=0)¶

• 行为：完全忽略反射因子数据
• 计算方式：采用DOE-2标准的太阳能增益计算流程
• 适用场景：不考虑周围表面反射贡献的常规分析

2. SOLAR-REFL-CALC=UNADJUSTED (ISREFL=1)¶

• 行为：直接用反射因子替代标准计算
• 计算方式：太阳能增益 = 反射因子 × 水平辐照度
• 特点：不做任何修正，完全依赖外部提供的反射数据
• 风险：可能重复计算某些辐射分量

3. SOLAR-REFL-CALC=ADJUSTED (ISREFL=2) - 最复杂的模式

• 核心思想：分离"直接辐射"和"反射辐射"的贡献
• 调整机制：
  • 在每月第一天，用反射因子减去DOE-2标准计算结果（已归一化到水平辐照度）
  • 得到纯反射分量的调整后反射因子
• 日常计算：
  • 调整后反射因子 × 水平辐照度 = 来自周围环境的反射辐射增益
  • 加上DOE-2标准计算的直接辐射增益
  • 总和 = 包含反射的净太阳能增益

通俗类比（C++背景）¶

可以把这理解为三种函数覆盖策略：

// ISREFL=0: 基类方法
class SolarCalc {
    double gain() { return standardDOE2Calc(); } // 只用基类实现
};

// ISREFL=1: 完全覆盖
class SolarCalc {
    double gain() { return reflectionFactor * irradiance; } // 彻底替换
};

// ISREFL=2: 装饰器模式（混合）
class SolarCalc {
    double gain() { 
        // adjustedRF = RF - standardDOE2Calc() (每月第一天计算)
        return standardDOE2Calc() + adjustedRF * irradiance; 
    }
};

ADJUSTED模式本质上是通过"做差法"分离出纯反射分量，避免重复计算，这是最物理真实的处理方式。

**核心子程序CALEXT（EXTERIOR CALCULATIONS），它负责每小时计算建筑外围护结构（墙体、门窗）的热工性能和太阳能得热。

整体架构（C++开发者视角）¶

// 伪代码：CALEXT的C++风格架构
class ExteriorSurfaceCalculator {
public:
    // 全局数据区（COMMON块的等价物）
    struct LocalData { double stalat, stalon, bazim, isrefl; /*...*/ };
    struct DaylightingData { /*...*/ };
    struct WallData { double soli, rgolge, filmu, /*...*/ };
    
    // 主入口
    void calculate(ZLSG_Array& zlsg, ZLSJM_Array& zlsjm, ZLSJY_Array& zlsjy) {
        initDaylightingVars();
        processWindowShading();      // 处理窗户遮阳
        processWalls();              // 处理墙体
        processDoors();              // 处理门
        // ...
    }
};

数据流特点：大量使用COMMON块（全局共享内存）和指针数组IA()/AA()（类似void*指针类型转换），这是FORTRAN时代的依赖注入。

核心计算模块¶

1. 窗户智能控制逻辑（第50-140行）¶

! 遮阳状态机：0=无遮阳, 1=打开, 2=关闭, 3=强制关闭
IA(MWI+96) = 1  ! 初始化

! 太阳能光学切换玻璃（电致变色）
IF(SVAR1 >= threshold) AA(MWI+129) = 0.0  ! 0=完全切换, 1=未切换

关键特征：

• 7种触发条件：直射辐射、总辐射、室外温度、负荷等
• 概率控制：引入随机数避免所有窗户同步动作
• 插值算法：在切换状态间线性插值可见光透过率

2. 太阳能反射计算（第170-240行）¶

实现了我们之前解释的三种模式：

IF(ISREFL.EQ.0) THEN
    ! NO-CALC: 只用标准DOE-2
    SOLABW = SOLI * AA(MP+5)  ! MP+5是表面吸收率
    
ELSE IF(ISREFL.EQ.1) THEN
    ! UNADJUSTED: 完全用预计算反射因子
    SOLABW = AA(MX+JJ+94)*RDNCC*RAYCOS(3) + AA(MX+119)*BSCC
    
ELSE IF(ISREFL.EQ.2) THEN
    ! ADJUSTED: 混合模式（物理最精确）
    SOLABW = SOLI * AA(MP+5) +  ! 标准DOE-2结果
             AA(MX+JJ+94)*RDNCC*RAYCOS(3) +  ! 反射贡献
             (AA(MX+119)-AWSKYD)*BSCC        ! 调整后的天空分量
ENDIF

关键公式：

• AWSKYD = SURABSO * (FFSW + FFG*GNDREF) - 基础反射分量
• RFSKYP = RF-AWSKY - AWSKYD - 纯反射增量

3. 响应因子法（RESPONSE FACTOR）¶

延迟墙体（热容材料）的核心算法（第2660-3400行）：

! 类似数字信号处理的IIR滤波器
XSXCMP = C_RATIO * (PREV_SUM_X) + SUMXDT + X(N)*DT
XSQCMP = C_RATIO * (PREV_SUM_Y) + SUMYDT + Y(N)*DT

! 热流计算
Q = (XSQCMP + C0*(T_surface - T_zone)) * AREA

物理意义：将墙体传热视为线性时不变系统，用预先计算的X/Y系数序列（响应因子）做卷积运算，避免有限差分法的时间步长限制。

4. Trombe墙专项处理（太阳能蓄热墙）¶

! 三层热阻网络
U1 = 1/(1/UCHAN + 1/UWW + 1/MP+14)  ! 玻璃到室内
U2 = 1/(1/UCHAN + 1/MP+14)          ! 玻璃到墙体
U3 = 1/(1/FILMW + 1/UWW)            ! 室外到墙体

关键数据表¶

现代C++重构建议¶

// 用现代C++概念封装
struct Surface {
    double azimuth, tilt;          // 方位角、倾角
    double area, absortivity;      // 面积、吸收率
    std::vector<double> rfFactors; // 反射因子
};

class SolarCalculator {
public:
    double computeAbsorbedSolar(const Surface& surf, 
                               const SolarData& solar,
                               int isreflMode) const {
        double baseGain = computeStandardDOE2(surf, solar);
        double reflected = computeReflectedContribution(surf, solar);
        
        switch(isreflMode) {
            case 0: return baseGain;
            case 1: return reflected;
            case 2: return baseGain + reflected; // ADJUSTED
        }
    }
};

调试支持¶

代码中遍布调试输出（如IREPRT(35)控制），可打印：

• 太阳几何参数（ETA, RAYCOS）
• 反射因子调整值
• Trombe墙内部温度分布
• 遮阳状态切换日志

这是80年代科学计算的典型风格：算法密集、全局状态、物理方程直接翻译。

AA(MX+19)在DOE-2代码中代表墙体的天空视角系数（Sky View Factor），是预计算的几何参数。出现NaN的情况通常源于输入数据异常或数值计算错误传播。

NaN产生的根本原因¶

1. 输入几何数据无效（最常见）¶

! AA(MX+19) 在输入阶段由建筑几何计算得出
! 如果输入文件中存在以下问题：
- 墙体面积 XSAREA = 0 或负数
- 墙体倾角 GAMMA = AA(MX+90) 超出 [0, π] 范围
- 方位角数据缺失或格式错误
- 墙体高度/宽度未正确定义

2. 计算错误传播链¶

! 典型传播路径示例：
GAMMA (NaN) → FFSIR (NaN) → FFGIR (NaN) → QIREW (NaN) 
           ↘ 存储到 AA(MX+19) 关联位置

关键代码位置：

第191行: IF (AA(MX+19).LT.0.) FFSIR = ABS(AA(MX+19))*0.5*(1.+GAMMA)
! 如果 GAMMA 为NaN，FFGIR = 1.-FFSIR 也会变成NaN

3. 指针/索引错误¶

MX = IA(MZ+5)  ! 墙体数据指针
! 如果 MZ 指针未正确初始化，MX可能指向无效内存区域
! 导致 AA(MX+19) 读取到未初始化的垃圾值

具体触发场景¶

诊断与调试方法¶

方法1：在CALEXT入口处添加NaN检查¶

! 在 SUBROUTINE CALEXT 开头添加
IF (isnan(AA(MX+19))) THEN
    WRITE(IOUTPT, *) 'ERROR: NaN detected in AA(MX+19)'
    WRITE(IOUTPT, *) 'Wall ID:', (IA(MX+II-1),II=1,4)
    WRITE(IOUTPT, *) 'XSAREA=', AA(MX+85), ' GAMMA=', AA(MX+90)
    CALL ABORT   ! 或赋默认值继续运行
ENDIF

方法2：启用调试输出¶

! 设置输入文件中的DUMP-OPTIONS
DUMP-OPTIONS = (DEBUG, LV-J)
! 这会激活第38-44行的调试WRITE语句
WRITE(IOUTPT,9300) ... AA(MX+19) ...

方法3：追踪数据来源¶

! AA(MX+19) 的来源是建筑几何定义阶段
! 检查输入文件中的EXTERIOR-WALL关键字定义：
EXTERIOR-WALL
    HEIGHT = 8
    WIDTH = 10
    TILT = 90  ! 倾角必须在0-180度之间
    AZIMUTH = 180
    ..

修复方案¶

临时修复（运行时防御）¶

! 在第210行附近添加保护代码
C------ CHECK FOR NaN IN SKY VIEW FACTOR
IF (isnan(AA(MX+19))) THEN
    AA(MX+19) = 0.5  ! 赋默认值（半球向天空）
    WRITE(IOUTPT, *) 'WARNING: MX+19 NaN replaced with 0.5'
ENDIF

根本修复（输入验证）¶

1. 检查所有墙体几何参数：

   • 确保XSAREA > 0
   • 确保0 ≤ GAMMA ≤ π (0°-180°)
   • 确保方位角在0-360度范围内

2. 验证输入文件完整性：

   • 使用DOE-2的BDL（Building Description Language）语法检查器
   • 确认所有EXTERIOR-WALL都有正确的几何定义

3. 更新到最新补丁：
   代码中已有FIXNAN补丁（如第24-27行），说明这是已知问题。检查是否应用了所有官方补丁。

C++开发者的经验类比¶

类似C++中的：

struct Wall {
    double skyViewFactor;  // AA(MX+19)
    double tilt;           // AA(MX+90)
};

// 问题等价于：
Wall* walls = (Wall*)malloc(...);  // 类似AA数组
// 如果walls[i].tilt = acos(invalid)，就会产生NaN
// 如果&walls[i]指向未初始化内存，就会读到垃圾值

核心教训：在读取外部输入的几何数据后，必须立即验证有效性，而不是等到计算时才暴露问题。