真空一维热传导稳态问题数值求解-VS2019工程

/* * 一维热传导稳态问题。 * 一柱状物体，一侧接触恒温源，另一侧悬空，处于真空中。 * 根据能量守恒和傅里叶传热方程，可得到热传导方程 * du/dt = kappa/(rho c) d^2u/dx^2 - epsilon sigma/(rho c) Cr/S (u^4-u0^4) * 由于是稳态，左边为0，故方程可化简为 * d^2u/dx^2 - epsilon sigma/kappa Cr/S(u^4-u0^4) = 0 * 最终化简为一阶常微分方程组 * du/dx = v * dv/dx = ab(u^4-c) * u(0)=T0 * v(L)=-a[u^4(L)-c] * 其中 * a = epsilon sigma/kappa * b = Cr/S * c = u0^4 * L为柱长，T0为热源温度, u0为环境温度, Cr为横截面周长, S为横截面面积，epsilon为材料表面发射率，kappa为热导率， * sigma为玻尔兹曼辐射常数，sigma = 5.67e-8 W/m^2 K^4 （单位为瓦/平米 摄氏度^4） * 注意求解时对原方程做了变换，将发射面记为x=0，恒温面记为x=L */ #define WIN32 #define _CRT_SECURE_NO_WARNINGS #define GSL_DLL #include <stdio.h> #include <gsl/gsl_odeiv2.h> #include <gsl/gsl_roots.h> #include <gsl/gsl_errno.h> /**********************************************************/ /* ODE test system definitions */ /**********************************************************/ struct parameters { double a; double b; double c; double L; double T0; }; int func(double t, const double y[], double f[], void* params) { struct parameters* p = (struct parameters*)params; f[0] = y[1]; f[1] = p->a * p->b * y[0] * y[0] * y[0] * y[0] - p->a * p->b * p->c; return GSL_SUCCESS; } int jac(double t, const double y[], double* dfdy, double dfdt[], void* params) { struct parameters* p = (struct parameters*)params; dfdy[0] = 0.0; dfdy[1] = 1.0; dfdy[2] = 4 * p->a * p->b * y[0] * y[0] * y[0]; dfdy[3] = 0.0; dfdt[0] = 0.0; dfdt[1] = 0.0; return GSL_SUCCESS; } double diff_sys(double T, void* params) { struct parameters* p = (struct parameters*)params; gsl_odeiv2_system sys = { func, jac, 2, p }; gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0); double t0 = 0.0, tL = p->L; double y[2]; y[0] = T; y[1] = p->a * y[0] * y[0] * y[0] * y[0] - p->a * p->c; for (int i = 1; i <= 100; i++) { double ti = i * tL / 100.0; int status = gsl_odeiv2_driver_apply(d, &t0, ti, y); if (status != GSL_SUCCESS) { printf("error, return value=%d\n", status); break; } printf("%.5e %.5e %.5e\n", t0, y[0], y[1]); } gsl_odeiv2_driver_free(d); return y[0] - p->T0; } /**********************************************************/ /* Main function */ /**********************************************************/ int main(void) { struct parameters p; p.a = 0.85e-10; // epsilon*sigma/kappa 发射率epsilon, 辐射常数sigma, 热导kappa p.b = 4.0; // Cr/S 横截面周长除以横截面积 p.c = 16.0; // u0^4 p.T0 = 100; // 恒温面温度 p.L = 1.0; // 杆长 int status; int iter = 0, max_iter = 100; const gsl_root_fsolver_type* T; gsl_root_fsolver* s; double r = 0, r_expected = 100; double x_lo = p.T0 - 100.0, x_hi = p.T0 + 100.0; gsl_function F; F.function = &diff_sys; F.params = &p; T = gsl_root_fsolver_brent; s = gsl_root_fsolver_alloc(T); gsl_root_fsolver_set(s, &F, x_lo, x_hi); printf("using %s method\n", gsl_root_fsolver_name(s)); printf("%5s [%9s, %9s] %9s %10s %9s\n", "iter", "lower", "upper", "root", "err", "err(est)"); do { iter++; status = gsl_root_fsolver_iterate(s); r = gsl_root_fsolver_root(s); x_lo = gsl_root_fsolver_x_lower(s); x_hi = gsl_root_fsolver_x_upper(s); status = gsl_root_test_interval(x_lo, x_hi, 1e-3, 1e-3); if (status == GSL_SUCCESS) printf("Converged:\n"); printf("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n", iter, x_lo, x_hi, r, r - r_expected, x_hi - x_lo); } while (status == GSL_CONTINUE && iter < max_iter); gsl_root_fsolver_free(s); return status; }