理工学基礎実験 熱の移動 学籍番号 61414050

sim.c

#define _USE_MATH_DEFINES

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define LENGTH_COPPER 0.15
#define RADIUS_COPPER 0.015

// The number of divisions in the vertical direction
#define VERTICAL_RESOLUTION 100
// The number of steps per second
#define TEMPORAL_RESOLUTION 1000
// Seconds of simulation
#if defined(EXP_A) || defined(EXP_B)
#define SIMULATION_LENGTH 450
#elif EXP_C
#define SIMULATION_LENGTH 960
#endif
// The initial temperature in degrees Celsius
#define INITIAL_TEMPERATURE 18.0

// Specific heat
#define SPECIFIC_HEAT 384.0
// Thermal conductivity
#define THERMAL_CONDUCTIVITY 401.0
#define DENSITY 8940

typedef struct copper {
	// index 0 is top.
	double temperature[VERTICAL_RESOLUTION];
} copper;

copper *c1;
copper *c2;
double mesh_length;
double mesh_volume;
double cross_section;
double mesh_weight;
double heat_capacity;

void prepare()
{
	int i;
	if ((c1 = (copper*)malloc(sizeof(copper))) == NULL ||
		(c2 = (copper*)malloc(sizeof(copper))) == NULL) {
		printf("malloc error\n");
		exit(1);
	}
	for (i = 0; i < VERTICAL_RESOLUTION; i++) {
		c1->temperature[i] = INITIAL_TEMPERATURE;
	}
	mesh_length = LENGTH_COPPER / VERTICAL_RESOLUTION;
	cross_section = M_PI * RADIUS_COPPER * RADIUS_COPPER;
	mesh_volume = mesh_length * cross_section;
	mesh_weight = mesh_volume * DENSITY;
	heat_capacity = mesh_weight * SPECIFIC_HEAT;
}

double heater(double time)
{
#ifdef EXP_A
    if (time < 180.0) {
        return 17.0;
    } else {
        return 0.0;
    }
#elif EXP_B
    return 17.0;
#elif EXP_C
    if (time < 240.0 || (480 < time && time < 720)) {
        return 17.0;
    } else {
        return 0.0;
    }
#endif
}

double cooler(double time)
{
#ifdef EXP_A
    return 0;
#elif EXP_B
    return -17.0;
#elif EXP_C
    return -8.5;
#endif
}

void simulate(double time, double time_length,
	copper *copin, copper *copout)
{
	int i;
	copout->temperature[0] = copin->temperature[0] + // original temperature
		(heater(time) + (copin->temperature[1] - copin->temperature[0])
			* cross_section * THERMAL_CONDUCTIVITY / mesh_length) // energy balance in Watt
		* time_length / heat_capacity;
	for (i = 1; i < VERTICAL_RESOLUTION - 1; i++) {
		copout->temperature[i] = copin->temperature[i] +
			((copin->temperature[i - 1] + copin->temperature[i + 1]
				- 2 * copin->temperature[i])
				* cross_section * THERMAL_CONDUCTIVITY / mesh_length)
			* time_length / heat_capacity;
	}
	copout->temperature[VERTICAL_RESOLUTION - 1] =
		copin->temperature[VERTICAL_RESOLUTION - 1] +
		(cooler(time) +
			(copin->temperature[VERTICAL_RESOLUTION - 2]
				- copin->temperature[VERTICAL_RESOLUTION - 1])
			* cross_section * THERMAL_CONDUCTIVITY / mesh_length)
		* time_length / heat_capacity;
}

int main()
{
	int i;
	copper *tmp;
	double time = 0.0;
	double time_length = 1.0 / TEMPORAL_RESOLUTION;
	prepare();
	while (time < SIMULATION_LENGTH) {
		simulate(time, time_length, c1, c2);
		tmp = c1;
		c1 = c2;
		c2 = tmp;
#ifdef EXP_C
        if (fmod(time, 1.0) < time_length) {
            printf("%d,%f,%f,%f\n",
                    (int)time,
                    c1->temperature[0],
                    c1->temperature[VERTICAL_RESOLUTION / 2 - 1],
                    c1->temperature[VERTICAL_RESOLUTION - 1]);
        }
#endif
		time += time_length;
	}
#if defined(EXP_A) || defined(EXP_B)
	printf("top: %f, middle: %f, bottom: %f\n",
		c1->temperature[0],
        c1->temperature[VERTICAL_RESOLUTION / 2 - 1],
        c1->temperature[VERTICAL_RESOLUTION - 1]);
#endif
	return 0;
}