美文网首页
利用蒙特卡罗法求解热传导方程

利用蒙特卡罗法求解热传导方程

作者: 深雨露 | 来源:发表于2020-01-17 21:09 被阅读0次

    蒙特卡罗方法简介

    蒙特卡罗（Monte Carlo, abbr. MC）方法是利用独立重复的统计实验来对物理及数学问题求解的方法。一个简单的关于MC方法的应用即求解图形的面积。

    例：求解(x-\frac{1}{2})^2+(y-\frac{1}{2})^2=\frac{1}{2}的面积？

    %投点次数
M = 1000;
A = rand(M,2);
r = sqrt((A(:,1) - 0.5).^2 + (A(:,2) - 0.5).^2);
total = sum(r < 0.5);
pr = total / M;
disp(pr)

%绘图
xt = @(t) 0.5*cos(t) +0.5;
yt = @(t) 0.5*sin(t) + 0.5;
fplot(xt,yt)
hold on
plot(A(:,1),A(:,2),'.');
grid on
axis equal 
axis([0,1,0,1])

    计算结果：0.7630


    计算结果

    MC求解面积的数学原理便是弱大数定理的方法

    设X_1,X_2,\cdots相互独立，服从同一分布的随机变量序列，且具有数学期望E(X_k) = \mu \quad (k=1,2,\cdots)，作前n个变量的算术平均\frac{1}{n}\sum_{k=1}^n X_k，则对任意的\varepsilon> 0，有
    \lim_{n \rightarrow \infty} \Pr \left\lbrace \left| \frac{1}{n} \sum_{k=1}^n X_k -\mu \right|<\varepsilon \right\rbrace =1.

    可以简单的理解为若统计的数据足够大，则事物的频率\frac{1}{n} \sum_{k=1}^n X_k可无限接近于其期望\mu。

    利用MC方法求解一维热传导方程的边值问题

    下面以2018年数学建模国赛A题为例，利用MC方法进行求解其第一问。

    题目：热防护服的设计

    人们在温度较高的环境下工作时，通常都要穿上专门的服装来避免灼伤。一般的避免灼伤专用服装是由三层的织物材料构成，分布标记为 I、II、III 层，而其中的 I 层与外界环境接触，III 层与皮肤之间还存在着空隙，所以我们将此空隙标记为 IV 层。为了降低研发专用服装成本、缩短研发设计专用服装的周期，首先把体内温度控制在 37ºC 的假人放置在实验室的高温环境中，然后建立数学模型来测量假人皮肤外侧的温度，并且解决下列问题：

    （1）附件1给出了专用服装材料的一些参数值，附件 2 给出了在环境温度为 75 ℃、II 层的厚度为 6 毫米、IV 层的厚度为 5 毫米、工作时间为 90 分钟的条件下假人皮肤外侧的温度的变化情况。请建立适当的数学模型来计算温度的分布，最后把计算好的温度生成为 Excel 文件。

    附件1. 专用服装材料的参数值

    分层 密度(kg/m^3) 比热 (J/(kg\cdotºC)) 热传导率 (W/(m·℃)) 厚度 (mm)
    I层 300 1377 0.082 0.6
    II层 862 2100 0.37 0.6-25
    III层 74.2 1726 0.045 3.6
    IV层 1.18 1005 0.028 0.6-6.4

    附件2. 假人皮肤外侧的测量温度

    时间 (s) 温度 (℃) 时间 (s) 温度 (℃)
    0 37.00 5396 48.08
    1 37.00 5397 48.08
    2 37.00 5398 48.08
    3 37.00 5399 48.08
    4 37.00 5400 48.08

    解析题目

    对于这种方法能量传递类问题，一般应考虑利用热传导方程来进行求解
    \frac{\partial u}{ \partial t} = k \cdot \nabla^2u,
    其中k是热传导系数。对于如题的一维问题，可简化上述方程为
    \frac{\partial u}{\partial t} = k_{\iota} \cdot \frac{\partial^2 u}{\partial x^2}, \quad x \in [0,L],\iota \in \{1,2,3,4\}
    其中k_i代表不同层次的热传导系数（如第 I 层的热传导系数k_1 = 0.082）。

    完善初始条件和边界条件

    初始条件
    u(x,0) = \varphi(x),
    边界条件
    u(0,t)=\mu_1(x), \quad u(L,t)=\mu_2(x).
    至于在两层的界面之间的联系则依据热流量密度k_\iota \frac{\partial u}{\partial x}连续的原理进行定义
    (k_\iota \frac{\partial u}{\partial x})^- = (k_{\iota+1}\frac{\partial u}{\partial x})^+
    下面使用MC方法求解热传导方程^{[1]}。取差分距离点：对距离x进行差分
    x_j = j \Delta x, \quad j = 1,2,\cdots;
    对时间t进行差分
    t_n = n \Delta t, \quad n=1,2,\cdots
    利用六点半隐式差分方法对原方程差分化
    \frac{{u_j^{n + 1} - u_j^n}}{{\Delta t}} = \frac{k}{{2(\Delta {x)^2}}}[(u_{j + 1}^{n + 1} - 2u_j^{n + 1} + u_{j - 1}^{n + 1}) + (u_{j + 1}^n - 2u_j^n + u_{j - 1}^n)]

    其中拉布拉斯算子采用(n+1)层和n层的中央差分和的平均

    (n+1)层
    \frac{\partial u}{\partial x} \approx \frac{(u_{j + 1}^{n + 1} - 2u_j^{n + 1} + u_{j - 1}^{n + 1})}{(\Delta x)^2}
    n层
    \frac{\partial u}{\partial x} \approx \frac{(u_{j + 1}^{n } - 2u_j^{n } + u_{j - 1}^{n })}{(\Delta x)^2}
    合并两层

    \frac{\partial u}{\partial x} \approx \frac{1}{2} \cdot \left[ \frac{(u_{j + 1}^{n + 1} - 2u_j^{n + 1} + u_{j - 1}^{n + 1})}{(\Delta x)^2} + \frac{(u_{j + 1}^{n } - 2u_j^{n } + u_{j - 1}^{n })}{(\Delta x)^2} \right]

    令 \lambda = k \frac{\Delta t}{(\Delta x)^2} 得

    - \frac{\lambda }{2}u_{j + 1}^{n + 1} + (1 + \lambda )u_j^{n + 1} - \frac{\lambda }{2}u_{j - 1}^{n + 1} = \frac{\lambda }{2}u_{j + 1}^n + (1 - \lambda )u_j^n + \frac{\lambda }{2}u_{j - 1}^n

    考虑 \lambda \leqslant 1
    \begin{align} u_j^{n + 1} &= \frac{1}{{1 + \lambda }}\left[ {\frac{\lambda }{2}u_{j + 1}^{n + 1} + u_{j - 1}^{n + 1} + \frac{\lambda }{2}u_{j + 1}^n + u_{j - 1}^n + 1{\rm{ - }}\lambda u_j^n} \right] \\ &={\alpha _1}u_{j + 1}^{n + 1} + {\alpha _2}u_{j - 1}^{n + 1} + {\alpha _3}u_{j + 1}^n + {\alpha _4}u_{j - 1}^n + {\alpha _5}u_j^n \end{align}
    其中\sum_{i=5}^5 \alpha_i =1 .

    依据上式可以得知，u_j^{n+1}的求得与u_{j + 1}^{n + 1} ,u_{j - 1}^{n + 1} ,u_{j + 1}^n ,u_{j - 1}^n ,u_j^n这五项及其相应对u_j^{n+1}的贡献权重/比例（即\alpha_i）密切相关。

    那么我们可设想人随机走动的模式：位于(j,n+1)处的人分别以概率\alpha_1,\alpha_2,\alpha_3,\alpha_4,\alpha_5向临近点移动。确定待移位置后，人到达待移位置，继续以类似的形式移动，直到人到达边界处（n=1或者边界点Q_{k1}）.将其赋值为f(Q_{k1}).

    原理图

    对一个u(x_j,t_{n+1})的子样f(Q_{k1})重复M次掷点，从而得到f(Q_{k1}),f(Q_{k2}),\cdots,f(Q_{kM})。那么可认为u(x_j,t_{n+1})的无偏估计为
    \widehat{u_M}=\frac{1}{M} \sum_{i=1}^{M} f(Q_{ki})
    依据上面的思路，即可求解出整个的热量传输过程的温度变化

    计算机实现步骤

    为加快求解速度，采用MATLAB中MATLAB Coder应用生成mex文件后（使用JIT编译器）进行生成。上述应用只能对函数进行编译，因而下文均在函数——MC_PDE中进行书写。下文以第IV层计算为例予以说明。

    运行实现

    定义函数的参数和返回值

    function u = MC_PDE(temp_90,M)
% temp_90 ：90分钟外界边界条件
% M = 10：投点次数

    初始条件设置

    %热传导方程系数
k = 0.028;
%内层边界条件（皮肤）
edge_first = repmat(37,1,5401);
%长度维起始行
edge_location = 1;
edge_location = int32(edge_location);
% 长度维终止点
l_end = 50;
l_end = int32(l_end);

%% x方向
% x方向划分间隔
h = 0.1;
% x方向：值域，x ∈ ( 0, x_max)
x_max = 15.2;
%x_initial = 0 : h :1;
% x方向迭代次数
x_times = x_max / h - 1;
x_times = int32(x_times);

%% t方向
% t方向划分间隔
tau =1;
% 时间值域：t ∈ ( 1, t_max)
t_max = 5401;
% 时间迭代次数
t_tims = t_max / tau;
t_tims = int32(t_tims);

temp_total = 0;

%系数求解
r = k * ( tau / (h^2) );
coef = 1 / (1 + r);

%%u初始化
u = zeros(x_times+1,t_tims) * nan;
% 初值条件
u(:,1) = 37;

% 边值条件
u(edge_location,:) = edge_first;
u(x_times+1, :) = temp_90(:,2)';

    随机走动部分

    for j_value = (edge_location + 1) : l_end
    disp(j_value)
    for n = 2:5401
        % 投掷点
        temp = zeros(1,M);
        for point_num = 1 : M
            t = int32(n);
            l = int32(j_value);
            while(1)
                xi = rand();
                if xi > 0 && xi < coef * (1/2 * r)
                    l = l + 1;
                    
                    %检验函数
                    if l == 152
                        basicvalue = u(152,t);
                        temp(point_num) = basicvalue;
                        signal = 1;
                    elseif   l == 1
                        basicvalue = u(1,t);
                        temp(point_num) = basicvalue;
                        signal = 1; 
                    elseif t == 1
                        basicvalue = u(l,1);
                        temp(point_num) = basicvalue;
                        signal = 1;
                    else
                        temp(point_num) = nan;
                        signal = 0;
                    end

                    if signal == 1
                        break
                    end

                elseif xi >= coef * (1/2 * r) && xi < coef * ( r)
                    l = l - 1;

                    %检验函数
                    if l == 152 || l == 1 || t == 1 
                        basicvalue = u(l,t);

                        temp(point_num) = basicvalue;
                        signal = 1;
                    else
                        temp(point_num) = nan;
                        signal = 0;
                    end

                    if signal == 1
                        break
                    end
                elseif xi >= coef * ( r) && xi < coef
                    t = t - 1;

                    %检验函数
                    if l == 152 || l == 1 || t == 1 
                        basicvalue = u(l,t);

                        temp(point_num) = basicvalue;
                        signal = 1;
                    else
                        temp(point_num) = nan;
                        signal = 0;
                    end

                    if signal == 1
                        break
                    end
                elseif xi >= coef && xi < coef * (1 + (1/2)*r)
                    l = l + 1;
                    t = t - 1;

                    %检验函数
                    if l == 152 || l == 1 || t == 1 
                        basicvalue = u(l,t);

                        temp(point_num) = basicvalue;
                        signal = 1;
                    else
                        temp(point_num) = nan;
                        signal = 0;
                    end

                    if signal == 1
                        break
                    end
                elseif xi >= coef * (1 + (1/2)*r) && xi < 1
                    t = t - 1;

                    %检验函数
                    if l == 152 || l == 1 || t == 1 
                        basicvalue = u(l,t);

                        temp(point_num) = basicvalue;
                        signal = 1;
                    else
                        temp(point_num) = nan;
                        signal = 0;
                    end

                    if signal == 1
                        break
                    end
                end
            end    
        end
        
        %期望计算部分
        temp_total = sum(temp);
        u(j_value, n) =  temp_total / M;
    end
end

    运行结果分析

    下面结论依据掷点数为1000的条件得到。运行时间为10803.82s（约3h）。

    三维图

    不可否认，由于MC方法求解的特殊性，温度随着时间和厚度的变化并不均匀。这一点可以从上图不断凸起的低矮“小山”这一现象中可知。

    二维热图

    二维热图可清楚可见三道明显的分界线，进热源点随时间温度变高，符合构建模型预期。

    总结

    1. 从计算复杂度看，复杂度略显偏大。求解速度呈现出求解区域靠近边界附近运算速度快，中间较为慢的特点，这个特点与这个算法的机制密切相关（因为中间的点难以足够快的速度移动到边界）。
    2. 从求解方法看，MC法相比隐式求解求解而言，对某个多维点而言，不必求解出各个层次的值，只需对这个多维点取足够多的子样即可求解。

    [1]刘春光. 偏微分方程边值问题的蒙特卡罗解法[D]. 吉林大学, 2004.

    相关文章

      网友评论

          本文标题：利用蒙特卡罗法求解热传导方程

          本文链接：https://www.haomeiwen.com/subject/kfshzctx.html

          延伸阅读

          深度阅读

          • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

          栏目导航

          热点阅读

          关于我们|服务条款|联系我们|利用蒙特卡罗法求解热传导方程|投稿指南|网站地图|RSS订阅|排版工具|手机版

          提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

          Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

          备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

          本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

          所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！