2014年西南交大数学建模新秀杯A题.docx

西南交通大学2014年新秀杯数学建模竞赛题目： A （填写A、B或C题） 组别： 大二组 （填写大一组或大二题） 参赛队员1参赛队员2参赛队员3姓名董吴卢学号学院信息学院信息学院机械学院专业通信工程通信工程机械设计制造及其自动化电话Email西南交通大学教务处西南交通大学实验室及设备管理处西南交通大学数学建模创新实践基地地热能源开发中热储层水温预测与钻井取水式系统评价模型摘要本文建立了地热能源开发中热储层水温预测与钻井取水式系统评价模型，从总体上就题给的三组钻孔深度与温度数据运用偏微分方程建模，再通过MATLAB软件处理，得到温度-深度变化关系，预测出约2000米深处的热储层温度并作了误差分析。查阅资料获得所需公式，依据公式分别估计出钻井取水至地表的水温度及流量，作了误差分析。最后评价并改进了模型，可以推广使用。对于问题一，首先，由附件1中所给出的数据利用matlab提供的数学工具，先拟合出在A、B、C三个钻口处40070米神的深度温度曲线，并求出“深度温度”方程。并作为特征方程用于以后的求解。为了便于对数据的分析，在认定岩石的热扩散率a1热导率与比热容c和密度之比：a=/（c其中：热导率（单位：W/(mK)）相同的情况下，我们利用热导率的偏微分方程及智能算法中的模拟退火算法求出在400-700各数据采样层面热泉通过位置的温度。最后估计出2000米处储热层温度大约为276.4000。对于问题二，首先，我们利用达西（Darcy）定律求出压强与水头差的关系，然后求出在地储层层面的水压。然后我们利用所查询的资料得知N80型号石油导管的特性，利用达西-魏斯巴赫公式求出水管内流动的沿程损失。然后再利用伯努利方程求得在地表泉口处的水流流速。最后，我们利用查询资料所找到的管道总传热计算公式求得水在管道中水温的变化规律，从而解得实际泉口处出水的温度。最后估计出地表泉口处的水流流量w=0.2959。而出水温度大约为=62。关键词：偏微分方程 最小二乘法 达西公式 MATLAB软件 LINGO软件 1 问题的重述1.1问题背景地热能源作为近几十年来的新兴能源越来越被世界所重视，目前，世界上已有80多国家直接利用地热能源。在我国，直接利用地热能源十分普遍，在开发地热能源前需要对热储层进行足够的分析，再结合实际情况决定是否值得开发以及怎样利用，因此准确预测地储层水温和评价获取热能的方式已经成为迫切解决的问题。1.2待解决问题1,2.1 问题1根据钻孔所得数据（见附件1）估计热储层水温，并分析物探测量埋深所带来的误差对水温估计值的影响；1.2.2 问题2通常可以通过打井的方式获得地热能。图1展示了一种在干热岩上打注水、出水两个井的方式。冷水从注水井进入，通过干热岩的加热，再从出水井流出，从而根据流出温度、流量决定如何利用。另一种更直接的办法是钻一个井直接打穿热储层，并在井中植入套管，让热水经由套管涌出地面。这种方式可以获得较大的热能，但有较大经济风险，因为实际中热储层位置并不容易把握，并且涌出热水的温度和流量是否能够达到预期要求也不好计算。现在想用后一种方式获得地热能，并采用型号为N80的石油套管（外径177.8毫米，壁厚8.05毫米），请估计出水温度和流量。2 问题的分析2.1 对两个问题的宏观分析 这两个问题大体上为数据的处理，拟合温度-深度间变化关系，最后得出数据的特征，然后由已有数据的特征对未来进行预测，并提前对相关部门提出建议。2.2 对两个问题的细化分析2.2.1 对问题一的分析 2.2.1.1 对深度-温度关系方程求解首先，由附件1中所给出的数据利用matlab提供的数学工具，先拟合出在A、B、C三个钻口处40070米神的深度温度曲线，并求出“深度温度”方程。并作为特征方程用于以后的求解。2.2.1.2 对热泉喷口处不同深度温度的求解 在这部分，为了便于对数据的分析，在认定岩石的热扩散率a1热导率与比热容c和密度之比：a=/（c其中：热导率（单位：W/(mK)）相同的情况下，我们利用热导率的偏微分方程及智能算法中的模拟退火算法求出在400-700各数据采样层面热泉通过位置的温度。2.2.1.3 对热储层温度的求解 在这里，我们利用在2.2.1.2中求出的热泉喷口处不同深度温度代入在2.2.1.1求出的“深度温度”特征方程中，求出在泉口喷口位置的深度温度方程。并代入泉口表面温度用来验证、修正所拟合曲线。最后代入埋深（物探等手段测得热储层深度）大致为2000米情况下的热储层温度。2.2.2 对问题二的分析2.2.2.1对地储层处水压的求解 首先，我们利用达西（Darcy）定律求出压强与水头差的关系，然后求出在地储层层面的水压。2.2.2.2 对泉口处水流速度的求解 然后我们利用所查询的资料得知N80型号石油导管的特性，利用达西-魏斯巴赫公式求出水管内流动的沿程损失。然后再利用伯努利方程求得在地表泉口处的水流流速。2.2.2.3 对泉口处水的温度的求解 最后，我们利用查询资料所找到的管道总传热计算公式求得水在管道中水温的变化规律，从而解得实际泉口处出水的温度。3 模型的假设3.1 在地储层上层的岩石的热物理特性相同。3.2 在地储层上层的岩石质密无空隙或孔隙大小对岩石的热物理性质无影响。3.3 岩石的热量仅来自热储层中水分的热量散失，不考虑水表面张力、黏度对能量传递造成的损失。23.4 储热层为无限大，其边界效应忽略不计。3.5 热储层的总能量保持恒定，且从热泉口流失的热能保持恒定。3.6 水的沿程水头损失可以近似视为04 符号说明_符号 代表含义x 相对深度t1 钻井A处各个深度的温度t2 钻井B处各个深度的温度t3 钻井C处各个深度的温度T 水的沸点温度h 所求位置的相对地表高度p 气体压强a 方程的特征系数t 在任意位置的岩层温度D（x） x-t方程的特征系数方差u(x,r) t关于x、r的关系方程 p 大气压强 g 重力加速度w 泉口水流流量v 泉口水流速度 地储层水温 用N80石油管道输送时地表水温d N80石油管道的内径_5 模型的建立与求解5.1 对问题一模型的建立与求解5.1.1 对问题一模型的建立对问题一我们主要建立用温度和深度两种指标的相关性来拟合此地区岩石的热物理特性。我们利用A、B、C三处的深度-温度数据，采取最小二乘法拟合出三处的深度-温度曲线方程，由拟合方程可以得出深度-温度的特征方程。以此特征方程可以代替热泉喷口处的深度-温度方程。然后我们代入在地表泉口处的修正温度（考虑到地表面大气压强，水与空气间的热量交换等因素后对温度进行逆求导）进行验证，修正方程。之后，我们采取先分段，再拟合的思路，利用二次中心差商公式以及热传导方程对泉口处400-700深度的各层面（每25米取一次数据点)的温度进行求解，然后利用先前求得的深度-温度特征方程得出泉口处温度-深度方程。在得出以上数据和方程后，利用岩石的热导率特性不变的条件，向泉口处温度-深度方程代入埋深（物探等手段测得热储层深度）大致为2000米情况下的热储层温度，并利用拟合求得的方程及其三维图像判断物探埋深对求解热储层温度的影响。5.1.2 数据预处理我们对A、B、C三点的深度-温度数据进行预处理，做成向量形式以便求取特征向量。处理前如下:1号钻孔2号钻孔3号钻孔钻孔深度（米）实测温度（摄氏度）钻孔深度实测温度钻孔深度实测温度40086400100400128.242591.6425102.442513245096.2450107450134.2475101.8475109475137500107500112500141.6525108525113525146550115.2550115.6550149575117.4575120.2575150.2600118600124600151.4625117625127625153650115.6650129650154.4675115.8675134675156700115.6700139.6700159处理后： 之后，我们对A、B、C三点的位置数据（A点的数据使用上文中筛选后的数据点）进行坐标化，建立了以泉口为原点的三维x-r圆柱坐标系。在此坐标系中得到了描述A、B、C三点深度、半径、温度的三个特征的向量，分别为：温度向量：tt=96.2 101.8 107 108 115.2 117.4 118 117 100 102.4 107 109 112 113 115.6 120.2 124 127 129 134 139.6 128.2 132 134.2 137 141.6 146 149 150.2 151.4 153 154.4 156 159深度向量：XX=450 475 500 525 550 575 600 625 400 425 450 475 500 525 550 575 600 625 650 675 700 400 425 450 475 500 525 550 575 600 625 650 675 700半径向量：RR=500 500 500 500 500 500 500 500 410 410 410 410 410 410 410 410 410 410 410 410 410 308 308 308 308 308 308 308 308 308 308 308 308 3085.1.3 对问题一的求解根据我们所查询到的资料和研究分析，我们得出：A、B、C三点的深度-温度拟合曲线：单位：depth(m);temperature()（程序详见附录）由matlab程序对三组曲线进行方程拟合，求得深度-温度曲线分别为：（详细过程见附录）拟合成一介线性方程的结果为： t1-x: t1=7.928x+-307 t2-x: t2=0.125x+ 49.15 t3-x: t3=0.1005x+ 90.25所得特征系数矩阵为： 由图像与方程可以分析出A的深度温度方程和曲线与B、C两点的不同，计算方差的公式为： 特征系数的方差D（x）=20.3595可见此时求得的方程的波动性较大，所以我们对A钻口的温度数据进行筛选，去掉了在400、650、675、700米深度的温度数据对A处的深度-温度数据进行二次拟合，求得方程如下：（详细过程见附录）t1-x: t1=0.1262x+42.25曲线如下：修正后所得特征系数矩阵为：此时的方差为D（x）=2.1036e-040.038528趋近于0；由均值公式为：求得= 0.1172由热传导的偏微分方程，我们可以解出在每一点温度关于半径、深度的限制条件：(此处为原始公式，x代表三维直角坐标系下x轴向)F（x,y,z）=c表示点(x,y,z)附近在单位体积与单位时间内产生的热量如果物体G自身任何位置都不产生热量，那么上式可简化为：由于我们建立圆柱坐标系，上式可以进一步化简为： （此处x为推导导出公式，x为符号说明中的相对深度）又因为除热源（热储层）外，其他物体都不产生热量，所以:F(r,y)=0 即根据泰勒级数展开，化简得中心差商公式：此处我们将以便于求导合并上面得到的两式，得到： 以这个式子作为约束条件并利用tt、XX、RR三组数据进行拟合曲线，用求得的（x的一介系数）和泉口地表温度=42.0摄氏度来修正曲线，得到方程如下（详细过程见附录）：所拟合图像如下：（XX：米；RR:米；tt：摄氏度）在所建立的圆柱坐标系z-r-中所拟合的图像如下： 向u(x,r)公式中代入数据为：x=2000.0 r=0.0解得：（详细过程见附录） = u(2000.0，0.0)= 276.4000 对于物探埋深的误差对结果估计的影响，由于在泉口位置的r=0，我们可以把u(x,r)函数简化为u(x,0)函数，即：拟合出u(x,0)的图像如下：由上述可知按照题目给出的物探埋深求出的地储层水温 =276.4000物探埋深的数值会直接影响对地储层温度的预测，且地储层预测温度与物探埋深深度（x）成正比，温度随深度的变化率为0.1172/m5.2 对问题二模型的建立与求解5.2.1对问题二模型的建立对于问题二，我们主要建立了在截面面积一定的情况下利用压强与流体速度的关系与管道内流体流动的沿程损失（达西-魏斯巴赫公式）运用matlab等数学工具计算出水流速度的变化关系，来求解地表表面的水流流速的方程，从而解得泉口流速。然后我们通过查阅资料与研究确定了管道总传热系数计算方法和在管道中流体流动的热损失方程。利用这两个公式，我们可以求得泉水温度变化与在中流动距离的关系方程。从而求得在泉口处泉水的水温。5.2.2 对问题二的求解根据我们所查询的资料可以得出水在沙土或砾石中的渗透规律：8达西（Darcy）定律：（其中为渗透率，为渗透路径长度，为水头差，为过水断面，为渗透系数）把公式进行进一步代换可得到： 以地处层面作为三维参考基准面。其中为水头差，为在地储层面产生的水柱高度，为在地表处产生的水柱高度在这里我们可以等效替代成在泉口处的大气压强所产生的高度。查询数据可知在一个标准大气压下产生水柱高度为10.34吗，即：=10.34m。由已知数据Q=50ml/s ; =2000m; k=0.5; A=0.78然后我们利用已经给出的数据Q(无管道时泉口处的水流速度)，（物探埋深）,k(岩层的渗透系数)代入上面化简得到的公式，求得：（详细过程见附录）=10.5964由达西-魏斯巴赫公式我们知道水在岩层中流动时有沿程水头损失，公式如下：其中：为流体的运动粘度，为流体的水头损失从而我们查询数据知道水的约为1.6*，所以水的水头损失很小，在此我们做理想化处理=0。由再假设能量守恒情况，水由管道从热储层流到地表表层后水的初速度与其最高能抵达的高度之间的关系满足关系如下：(重力势能)+(动能)=E E不变化，所以可以推导出如下关系： 代入上述数据可以解出： v=14.4114 m/s再代入水流速与面积和流量间公式： s= （其中w为流量，s为截面面积，v为截面流速）由查询到的N80管的数据可知其内径为：d=161.7 mm代入数据得到： w=0.2959然后根据查询到的管道总传热系数计算方法和热量平衡方程可得如下计算表达式：7 （1-1）【式中：总传热系数，W/(m2)；计算直径，m；（对于保温管路取保温层内外径的平均值，对于无保温埋地管路可取沥青层外径）；管道内直径，m；管道最外层直径，m； 油流与管内壁放热系数，W/(m2)； 管外壁与周围介质的放热系数，W/(m2)； 第层相应的导热系数，W/(m)； ，管道第层的内外直径，m，其中；结蜡后的管内径，m。】由上式，我们首先确定内部放热系数：在此模型中我们先假设在层流状态（Re t1=86 91.6 96.2 101.8 107 108 115.2 117.4 118 117 115.6 115.8 115.6t1 = Columns 1 through 9 86.0000 91.6000 96.2000 101.8000 107.0000 108.0000 115.2000 117.4000 118.0000 Columns 10 through 13 117.0000 115.6000 115.8000 115.6000 t2=100 102.4 107 109 112 113 115.6 120.2 124 127 129 134 139.6t2 = Columns 1 through 9 100.0000 102.4000 107.0000 109.0000 112.0000 113.0000 115.6000 120.2000 124.0000 Columns 10 through 13 127.0000 129.0000 134.0000 139.6000 t3=128.2 132 134.2 137 141.6 146 149 150.2 151.4 153 154.4 156 159t3 = Columns 1 through 9 128.2000 132.0000 134.2000 137.0000 141.6000 146.0000 149.0000 150.2000 151.4000 Columns 10 through 13 153.0000 154.4000 156.0000 159.0000 A=polyfit(x,t1,2);z=polyval(A,x);plot(x,y,r*,x,z,b)Undefined function or variable y. A=polyfit(x,t1,2);z=polyval(A,x);plot(x,t1,r*,x,z,b) hold on B=polyfit(x,t2,2);n=polyval(B,x);plot(x,t2,r*,x,n,b) hold on C=polyfit(x,t3,2);m=polyval(C,x);plot(x,t3,r*,x,m,b) p=polyfit(x,t1,1)p = 0.0996 53.3099 cftool-去掉部分数据后A的拟合 t11=91.6 96.2 101.8 107 108 115.2 117.4 118 117t11 = 91.6000 96.2000 101.8000 107.0000 108.0000 115.2000 117.4000 118.0000 117.0000 cftool xx=425 450 475 500 525 550 575 600 625xx = 425 450 475 500 525 550 575 600 625-求系数方差a=7.928 0.125 0.1005a = 7.9280 0.1250 0.1005 var(a) ans = 20.3595-求修正方差 b=0.1262 0.125 0.1005b = 0.1262 0.1250 0.1005 var(b)ans = 2.1036e-04-求系数a均值 average=mean(b)average = 0.1172-求解u(x,r)方程Linear model Poly22: f(x,y) = p00 + p10*x + p01*y + p20*x2 + p11*x*y + p02*y2Coefficients (with 95% confidence bounds): p00 = 42 (fixed at bound) p10 = 0.1172 (fixed at bound) p01 = 0.2939 (0.231, 0.3567) p20 = 0 (fixed at bound) p11 = -1.143e-05 (-9.517e-05, 7.23e-05) p02 = -0.0005804 (-0.0006754, -0.0004853)Goodness of fit: SSE: 1909 R-square: 0.8309 Adjusted R-square: 0.82 RMSE: 7.847-求解u（2000,0） r=0r = 0 x=2000x = 2000 u=42+0.1172*x+0.2939*r+(-1.143e-05)*x*r-0.0005804*r2u = 276.4000-求地壳不同深度的大气压强及水的沸点h=-2000p=exp(-h*1