Numerical

main.c

#include <stdio.h>
#include <omp.h>

#define W 100
#define H 100
double diff(double backward, double material, double forward);
double one_side_diff(double wall, double material);

int main() {
    double Temperature = 10.0; // [℃]
    double delta = 3.0e-3; // [m]
    double delta_t = 1.0e-1; // [sec]
    int LoopNum = 6*60*60 / delta_t; // 6 [hour]
    double Q = 100 / 0.07; // [w/m3]
    static double T_map[H][W]; // [℃]
    static double L_map[H][W]; // [W/m K]
    static double C_map[H][W]; // [J/m3 K]
    static double T_diff[H][W];
    int i, j, iter;
    double av_temp = 0;

    FILE *fp;
    if ((fp = fopen("result.csv", "w")) == NULL) {
        printf("File open error!\n");
        return -1;
    }

    for (i = 0; i < H; ++i)
        for (j = 0; j < W; ++j) {
            T_map[i][j] = Temperature;
            if(i>=26 && j<W-26){
                if(i>=33 && j<W-33){ // 人
                    T_map[i][j] = 37.0;
                    L_map[i][j] = 0.582;
                    C_map[i][j] = 4184000;
                }else{ // 毛布
                    L_map[i][j] = 0.15;
                    C_map[i][j] = 80000;
                }
            }else{ // 掛け布団
                L_map[i][j] = 0.04;
                C_map[i][j] = 98000;
            }
        }

    printf("Start. (LoopNum: %d)\n", LoopNum);
    for (iter = 0; iter < LoopNum; ++iter) {
#pragma omp parallel
        {
#pragma omp for
            for (i = 0; i < H; ++i)
                for (j = 0; j < W; ++j) {
                    if(i!=H-1){
                        T_diff[i][j] = diff(
                            // 温度一定条件
                            i != 0 ? T_map[i - 1][j] : Temperature,
                            T_map[i][j],
                            T_map[i + 1][j]
                        );
                    }else{
                        T_diff[i][j] = one_side_diff(
                            // 断熱境界条件
                            T_map[i][j],
                            T_map[i-1][j]
                        );
                    }
                    if(j != 0){
                        T_diff[i][j] += diff(
                            T_map[i][j - 1],
                            T_map[i][j],
                            // 温度一定条件
                            j != W - 1 ? T_map[i][j + 1] : Temperature
                        );
                    }else{
                        T_diff[i][j] += -one_side_diff(
                            // 断熱境界条件
                            T_map[i][j+1],
                            T_map[i][j]
                        );
                    }
                     
                    T_diff[i][j] /= delta*delta;
                    T_diff[i][j] *= (L_map[i][j]*delta_t);

                    // 人体の発熱
                    if(i>=33 && j<W-33) {
                        T_diff[i][j] += Q*delta_t;
                    }

                    T_map[i][j] += T_diff[i][j] / C_map[i][j];
                }
        }
    }
    printf("Finished.\n");

    iter = 0;
    for (i = 0; i < H; ++i)
        for (j = 0; j < W; ++j) {
            fprintf(fp, "%f", T_map[i][j]);
            if(j!= W-1) fprintf(fp, ",");
            else fprintf(fp, "\n");

            // 対表面温度の平均を計測
            if(
                (i==33 && j<W-33) ||
                (j==W-33 && i>=33)
            ){
                av_temp += T_map[i][j];
                iter++;
            }
        }

    printf("surface tmp.: %f\n", av_temp/iter);

    return 0;
}

double diff(double backward, double material, double forward) {
    return (backward - 2 * material + forward);
}
double one_side_diff(double wall, double material) {
    return -2*(wall - material);
}

viewer.py

#! /usr/bin/env python
# coding: utf-8

import csv
import numpy as np
from matplotlib import pyplot as plt

def main():
    T_map = []
    with open("result.csv", "r") as f:
        reader = csv.reader(f)
        for row in reader:
            tmp = []
            for col in row:
                tmp.append(float(col))
            T_map.append(tmp)
    
    # print(T_map)
    plt.imshow(T_map)
    plt.colorbar()
    plt.show()

if __name__ == '__main__':
    main()