Matlab数值分析

Matlab数值分析（精选7篇）

Matlab数值分析 篇1

0 引言

天线设计过程中,通常采用单个天线就能满足需求,但在一些特殊的情况下会对天线结构及参数提出更高的要求,例如高增益、低旁瓣,波束可控性等,此时就需要采用阵列天线的形式[1]。目前,阵列天线在军事通信、飞行遥测等方面得到了广泛的应用,随着通信的不断发展,如何完善阵列天线的结构以及提升天线的性能成为了研究的热点[2]。

Matlab是一款功能强大的数学软件,高效的数值计算功能使人们从繁杂的数学运算分析中解脱出来,同时可实现计算结果和编程的可视化[3]。Matlab以其卓越的数值计算能力和强大的绘图功能,近年来被广泛应用在天线的分析和设计中[4]。

本文将Matlab引入阵列天线设计中,利用其高效数值计算功能为天线设计提供依据。借助Matlab可以绘制出阵列天线的二维和三维方向图,直观地从方向图中看出天线性能随各参数的变化情况,从而为阵列天线设计提供依据。

1 阵列天线理论基础

阵列天线是由许多形式相同并按一定规律排列的辐射单元所组成的辐射系统,辐射单元称为阵元或单元天线。单元天线可以是贴片天线、印刷振子天线、缝隙天线等不同天线形式。阵列天线特性主要由阵元性质、阵元个数、阵元位置、阵元排列方式以及电流幅度和相位分布等因素决定[5]。阵列天线的形式多样,针对不同的工程需要。可以设计不同的阵列天线结构,阵列天线结构参数多,天线设计比较复杂。本文设计的阵列天线中天线单元如图1所示,图1是最常用的矩形微带天线及其辐射场,其所示的天线结构是阵列天线中最常用的单元天线形式[6]。

阵列天线的设计比较灵活,例如为了提高天线的增益,可以不断增加阵元的个数。阵列天线的设计过程中必须注意一些规范,例如阵列天线中阵元间的距离一般要满足在0.5个波长到1个波长之间,如果间距过小,阵元间的互耦会增大,如果间距大于1个波长,阵列天线又会出现栅瓣[7]。因此,如何设计一个理想的阵列天线结构,使其可以满足特定的设计需求,需要一个仔细调试的过程。天线的经典理论是天线设计过程中必须注意的部分[8],天线谐振时,谐振长度为:

式中:εe为等效的相对介电常数;h为基片厚度;L为贴片长度;W为贴片宽度。

以往人们在阵列天线设计过程中,阵列天线结构复杂、参数众多,天线模型的建立以及天线的调试都是基于设计者的经验,因此存在一定的弊端,很难综合考虑天线整体的参数性能。本文在经典天线理论的基础上,利用Matlab数值计算进行辅助设计,归纳天线的性能与各项参数之间的关系,从而高效地进行天线的设计。

2 阵列天线数值分析

2.1 天线单元分析

实际的天线设计过程中,天线性能与天线结构各项参数密切相关,随着天线领域的不断发展,对天线的尺寸、性能等都提出了更高的要求,但天线的各项性能与结构参数之间是相互制约的,因此结构参数与天线性能之间往往需要折中考虑[9]。阵列天线由于是许多天线单元组合的一个辐射系统,因此阵列天线的设计更加复杂。

图1中,W的尺寸影响天线整体尺寸的同时影响天线的辐射电阻、输入阻抗以及方向性函数,因此影响天线的频带宽度和辐射效率。适当增大W尺寸有利于天线的频带、效率及阻抗匹配,但当W尺寸大于半波长时将产生高次模,从而引起场的畸变。利用Matlab对天线理论公式进行数值计算,可以将复杂的函数关系直观的显示出来。

902~928 MHz频段是目前射频领域常用的频段,根据要求设计一个四单元阵列天线,要求天线尺寸不超过360 mm×360 mm的范围。根据上述设计需求,天线单元的贴片宽度必须满足w

mm

通过计算结果可以看出,在εr=1.0的情况下,W对L的影响不大,h值较小时,L值偏大。在考虑天线性能的同时还要考虑天线的尺寸限制,选取h为18 mm或19 mm时较为妥当。综合比较,h取18 mm,W为130 mm,L为140 mm,此时的天线单元结构比较理想。

2.2 阵列参数分析

阵列天线是多个天线单元共同作用的结果,考虑单个天线单元性能的同时必须考虑整体的性能。阵列天线结构复杂、参数繁多,因此阵列天线的设计过程中,如何在各个参数与性能之间折中考虑,往往花费大量的时间[10]。通过Matlab辅助分析阵列天线辐射方向图,可以将复杂的函数以二维、三维图形直观的展示出来,极大地方便了对阵列天线辐射场空间分布的理解,从而为阵列天线的设计提供参考思路。

在满足可生产性要求的前提下,为了使天线调试便捷、馈线网络简单,阵列天线采用等幅同相馈电网络。

通过运用Matlab进行阵列天线远场辐射方向图分析时,把单元天线当做点源来考虑。图2为四单元点源分布图,相位中心为天线阵列的几何中心,横向和纵向间距分别为2dx和2dy。

四个单元等幅并且同相馈电的情况下,四单元阵列天线的辐射远区总场可以表示为:

其中,E0为点源在无限远区形成的辐射场。

把单元天线当作点源来考虑,主要考虑天线单元的排列方式对天线性能的影响,排列方式是指天线单元之间的横向间距和纵向间距,通过Matlab对式(4)进行数值计算,运算的结果如图3所示。

数值计算结果直观的显示出了阵列天线辐射场的空间分布,由数值计算结果可以看出,间距dx,dy由80 mm增大到110 mm过程中,波瓣变窄,旁瓣电平随之抬升。在考虑天线性能以及天线尺寸的情况下,间距dx,dy取90~100 mm较为妥当。

3 天线仿真结果

根据Matlab对阵列天线结构的运算结果,在HFSS中按照运算结果进行阵列天线的模型设计,HFSS中仿真的四单元阵列天线图如图4所示。

HFSS中所建立的阵列天线仿真模型,单元天线为矩形微带天线,贴片与反射面之间为空气介质。通过对天线结构参数的适当调整,天线性能达到理想状态,此时天线结构参数h为18 mm,W为128 mm,L为140 mm,单元天线的间距dx=dy=96mm,此时天线增益达到12.94 dBi。天线的三维方向图以及回波损耗如图5所示。

由于仿真环境存在差别,对于Matlab数值计算的结果需要适当调整,由天线仿真结果来看,介质厚度h对天线的带宽影响比较大,带宽与h值成正比例关系。W的尺寸影响着天线的方向性函数,随着W尺寸的变大,天线的频带宽度和辐射效率都有所改善,相比之下,L对天线的影响较小。L,W共同决定天线的尺寸;在保持单元天线尺寸不变的情况下,随着单元天线间距dx,dy的变大,天线增益变大,带宽减小。

由仿真结果来看,HFSS仿真过程中天线的结构参数、天线性能的变化趋势与Matlab理论计算一致,Matlab的数值计算给天线设计提供了快速有效的指导意义。

4 结语

本文在阵列天线设计中引入Matlab进行辅助设计,利用其高效的数值计算能力和强大的绘图功能,从而把复杂的理论函数关系以数值曲线、方向图的形式直观的显示出来,因此可直观地了解阵列天线性能随着天线结构参数的变化情况,最终为阵列天线设计提供快速有效的指导意义。实验结果表明,利用Matlab进行辅助设计可以满足特定条件下阵列天线高效、快速的设计要求。

参考文献

[1]POZAR D M,SCHAUBERT D H.Microstrip antennas:Theanalysis and design of Microstrip Antennas and arrays[M].[S.l.]:Wiley IEEE Press,1995.

[2]邵一鹏,白明,马慧瑾,等.ETC用5.835GHz微带阵列天线的设计[J].现代电子技术,2012,35(4):148-151.

[3]陈天禄,郭燕红.阵列天线方向图的Matlab实现[J].西藏大学学报:自然科学版,2010,25(1):103-107.

[4]肖磊,蔡君.Matlab在相控阵天线中的应用[J].青海大学学报:自然科学版,2004,27(4):42-48.

[5]汪茂光,李周淼,吴黎明.阵列天线分析与综合[M].成都:电子科技大学出版社,1999.

[6]阮成礼.超宽带天线理论与技术[M].哈尔滨:哈尔滨工业大学出版社,2006.

[7]方大纲.天线理论与微带天线[M].北京:科学出版社,2006.

[8]林昌禄,聂在平.天线工程手册[M].北京:电子工业出版社,2002.

[9]刘燕.阵列天线方向图综合算法研究[D].西安:西北工业大学,2007.

[10]魏宏亮,段文涛,李思敏.运用Ansoft HFSS设计圆极化微带天线阵[J].现代电子技术,2008,31(1):71-73.

Matlab数值分析 篇2

通常求解结构动力响应的方法主要有:反应谱分析方法和时程分析方法。用反应谱理论来计算结构地震反应的分析方法, 其优点是计算简单且可以同时考虑多个地震地面运动激励的影响, 但其缺点是所求得的内力与位移的最大值并不是真实的, 无法求得内力与位移随时间发展的全过程, 同时很难应用于非线性分析。

时程分析法可以进行线性或非线性分析, 并且可以对结构的内力与位移进行全过程分析, 可以观察到结构内力与位移达到峰值点的时刻。

基于振型分析的时程分析方法, 求解多自由度结构体系的动力反应时, 在求出结构的自振频率和振型的基础上, 需要用到振型叠加原理。通常, 时程分析法在求解结构弹性地震反应时, 第一种方法是采用振型分析法求出结构的自振频率与振型, 再通过叠加原理求出结构的地震反应。第二种方法是直接从运动微分方程出发, 用数值方法直接进行积分, 得出结构的地震反应[1]。

本文通过运用几种不同的数值方法, 对结构的地震反应进行分析, 并通过Matlab语言编写相应的程序将两种做法的计算结果进行对比。最后, 通过一算例说明运用Matlab语言程序的方便性和有效性。

2 振型分解法基本原理[2]

由于多自由度结构体系在动力作用下是一组耦联的微分方程组, 所以振型分解法就是利用各振型相互正交的特性, 将耦联的微分方程组解耦为相互独立的微分方程组, 从而使求解原来多自由度的结构体系变为求解若干个广义单自由度体系的问题, 最后对解进行组合可得到多自由度结构体系的反应。

层间剪切型多自由度结构体系的运动方程为:

式中, [M], [C]和[K]分别为体系的质量、阻尼和刚度矩阵。

由于结构的阻尼是非常复杂的, 包括材料的阻尼、连接的阻尼等, 有的能量耗散还不一定是黏滞阻尼, 因而常常给不了阻尼矩阵。因此, 在实际工作中也常常是通过阻尼比, 求出结构的阻尼矩阵。本文采用瑞雷阻尼矩阵, 即:

式中, α, β分别为未知常数。

将u (t) 写为:

式中, U, q (t) 分别为振型矩阵和广义坐标。

将式 (2) 和式 (3) 代入式 (1) , 且在两边都左乘UjT, 则可得:

式中, γj, ωj为第j个振型的振型参与系数和频率。

取前两阶频率由式 (5) 可以确定两未知常数α和β:

结构的反应就可以通过式 (6) 、式 (7) 求得:

3 单自由度结构体系弹性反应分析

根据达朗贝尔原理, 物体运动的任一瞬时, 作用在物体上外力与惯性力相互平衡, 可得:

式中, 为地面运动加速度。

3.1 Duhamel积分计算方法

在整个荷载作用下, 任一时刻的位移反应, 可以看成是各微小位移反应之和, 可得:

对式 (9) 微分后, 可得速度与加速度:

分别对式 (9) 、式 (10) 及式 (11) 进行积分, 可得结构的位移、速度及加速度反应。

3.2 Newmar k-bet a法计算原理

Newmark-beta法是线性加速度法的一种推广, 可按式 (12) 、式 (13) 进行计算:

式中, α与δ为参数, 可调节。当α=1/6, δ=1/2时, 即为线性加速度法。借助式 (12) 、式 (13) 及式 (8) 可求出结构的速度、加速度及位移反应。

3.3 Rong-kut a法计算原理[3]

随着精度要求的不同, Rong-kuta法的计算公式也不相同。本文采用的算式, 取误差与Δt 5同阶, 此时的Rong-kuta法为:

式中,

4 计算与分析

某3层钢筋混凝土结构[4], 结构各层的参数为:第1~3层的质量分别为2762, 2760, 2300kg, 第1~3层的刚度分别为2.485×104, 1.921×104, 1.522×104N/m。地震波选用El-Centro波, 时间步长为0.02s, 且峰值调整为200cm/s2。

分别采用3种方法进行程序设计, 程序流程图如图1所示。 (1) 用Duhamel积分[5]方法进行分析计算, 在求解式 (9) 的积分时, 利用Matlab的内部函数conv进行求解。对两个长度都为m的向量, 利用conv函数求卷积可得到一个长度为2m-1的向量, 在求解卷积后, 取m前项作为计算结果。 (2) 用Newmark-Beta方法进行求解。 (3) 用Rong-kuta方法求解。3种计算方法的计算结果是一致的, 各层的位移、加速度时程曲线分别如图2、图3所示。

对采用Rong-kuta方法求解与文献[3]通过直接积分的Wilson-theta法进行对比。

图4为结构第三层的位移和加速度时程对比。由图4可知, 文中所用的Rong-kuta法与文献[3]所用的直接积分的Wilson-theta方法完全一致, 说明本文采用的Rong-kuta法是正确有效的。对于Duhamel积分分析方法和Newmark-Beta方法, 笔者已与文献[3]中的Wilson-theta方法进行了对比, 证明是正确有效的, 由于篇幅所限, 此处不再一一列出。不同方法计算时间对比如表1所示。

注:所用计算机的CPU型号:Intel奔腾双核T4200, 计算所用El-centro波的步数为400, 每步时间均为0.02s, 表中所列时间均为30次计算所费时间的平均值。

由表1可知, Wilson-theta方法计算所用时间是最长的, 主要是因为它是采用直接积分的计算方法, 每一时间步都要对矩阵进行乘与除的运算, 计算量较大。Duhamel积分方法也是效率比较低的, 且由于本文所用的conv函数本身的特点, 需要在计算完成后取其一半的有效计算结果, 因此, 也较费时。

5 结论

文中采用3种方法分别编制的几个程序都是基于单自由度结构体系, 因此, 在用于求解多自由度结构体系时需对结构体系进行振型分解。不同的数值分析方法的求解效率是不同的, 直接积分法计算量最大, 因此耗时最长。Duhamel积分的效率也相对较低, 主要是由于在计算卷积时需进行大量的乘法计算, 且由于conv函数本身的特点计算完成后要砍掉后半部分内容, 相当于做了一部分无用功。

很多文献对于几种数值计算方法都有介绍, 但基本都是采用Fortran或C语言编写, 循环较多、复杂且比较难懂。本文通过Matlab语言把多种数值方法都编写成了程序, 由于Matlab语言的优良特性使得编写程序变得比较容易。结果表明, 计算结果是正确有效的。

摘要：通过Matlab语言编写了3个时程分析程序, 分别采用Duhamel积分法、Newmark-beta法及Rong-kuta法。在求解多自由度结构体系时, 采用振型分解法将多自由度结构体系问题转化为多个广义单自由度结构体系问题, 然后对所得计算结果按一定准则进行叠加, 求出结构地震反应。首先, 对某3层钢筋混凝土框架结构进行时程分析, 计算结果与已有文献吻合较好, 证明了所编程序的正确性;其次, 通过对计算时间的分析对比, 表明Rong-kuta法效率最高。论文可为相关工程及设计人员在选择计算方法时提供参考。

关键词：Matlab程序,时程分析,振型分解法

参考文献

[1]Clough R.W, Penzien J.Dynamics ofstructures[M].NewYork:McGraw-Hill, 1993.

[2]Bathe K.J.Finite element procedures[M].Englewood Cliffs:Prentice hall, 1996.

[3]朱伯龙, 张琨联.建筑结构抗震设计原理[M].上海:同济大学出版社, 1993.

[4]徐赵东, 郭迎庆.Matlab语言在建筑抗震工程中的应用[M].北京:科学出版社, 2004.

Matlab数值分析 篇3

本文对Matlab中的各种数值积分函数进行了说明和分析, 并且通过对不同函数的比较总结了各个数值积分函数的优缺点。

1. Matlab数值积分函数的分析

1.1 基于梯形算法的积分函数

在Matlab中trapz (x, y, dim) 表示梯形法数值积分, 通过已知参数x, y按dim维使用梯形公式进行积分。向量x和y有相同的长度, (xi, yi) 代表曲线上的一点。曲线上点的距离不一定相等, x值也不一定有序。

另外, Matlab中函数cumtrapz (x, y, dim) 表示累积的梯形积分, 即不断的从第一个累积到当前的结果, 它的返回值是与x维数相同的向量。

1.2 基于Simpson算法的积分函数

quad (fun, a, b, tol, trace) 用自适应Simpson公式[3]进行数值积分的, 适用于精度要求低、平滑性较差的被积函数。其参数中, 被积函数fun必须是函数句柄;积分限[a, b]必须是有限的, 因此不能是inf;参数tol表示绝对误差, 默认值为10-6;当trace是非零值时则在递归过程中显示[fcnt a b-a Q]的值。值得注意的是, 该算法当积分递归计算超过10000次或当函数有不可积的奇点时, 则返回inf。

另外, Matlab中函数quadv (fun, a, b, tol, trace) 是用矢量化自适应simpson法进行数值积分的, 该函数是将quad函数矢量化, 然后一次可以计算多个积分, 所有的要求完全与quad相同。

1.3 基于Lobatto算法的积分函数

函数quadl (fun, a, b, tol, trace) 是用自适应Lobatto法进行数值积分的, 适用于精度要求高, 被积函数曲线比较光滑的数值积分。其中参数说明及注意事项与quad函数相同。

1.4 基于Gauss-Kronrod算法的积分函数

函数quadgk (fun, a, b, param1, val1, param2, val2, ...) 是用自适应Gauss-Kronrod法进行数值积分, 适用于高精度和震荡数值积分, 支持无穷区间, 并且能够处理端点包含奇点的情况, 同时还支持沿着不连续函数积分, 复数域线性路径的围道积分法。其参数中, 被积函数fun必须是函数句柄, 积分限[a, b]可以是[-inf, inf], 但必须快速衰减, param, val为函数的其它控制参数。

2. 几种数值积分函数的比较

从上述六个函数中选取具有代表性的trapz、quad、quadl以及quadgk四个函数进行分析。

2.1 基于精度的分析

使用这四个函数分别计算积分, 积分区间取[0, 1], 在Matlab中运行, 得出Matlab命令与计算结果, 见表 (2-1) 。

由表可看出, 函数trapz是定步长计算积分的, 它可以通过减小步长来提高精度, 但无法由该方法直接知道当前步长下的精度。而后面三个函数是由精度来控制运算的, quad及quadl是利用递归算法, 当达到指定精度时停止计算。由Matlab中这三种算法本身的实现过程可知, quad及quadl算法的默认精度都是10-6, 而quadgk的默认精度为10-10。

另外, 由于函数quad是用simpson算法实现的, 即用抛物线来代替被积函数, 只有三次代数精度[4], 所以对于光滑性[5]较差的函数具有较好的运算结果;而函数quadl是用Lobatto算法实现的, 故更适用于光滑性较好的函数。

2.2 基于运算速度的分析

由上面的研究可知, 用quad、quadl及quadgk进行数值计算时, 对运算结果的精度都能有较好的控制。为了进一步研究这三个函数的优缺点, 下面我们研究它们的运算速度随精度的变化。改变这三个函数的参数设置, 让它们的误差从10-4变化到10-16, 记录耗时随精度的变化情况, 得出图 (2-1) 。

分析图2-1发现, 当对精度的要求较低时, 利用函数quad计算积分耗时最短, 而随着精度增加, 它的耗时增加很快。故当精度要求较高时, quad函数的效率较低, 不建议继续使用。

观察函数quadl及quadgk (2009) 的耗时随精度的变化曲线发现的, 这两个函数的运算速度基本稳定, 但quadl当精度高于10-13时, 耗时在逐步增加。故当精度要求在10-6与10-13之间时, 建议读者优先考虑使用quadl函数计算积分。

2.3 基于适用范围的分析

由于函数quadgk在Matlab2009版才推出, 所以普及程度不够高, 但无论从精度、速度或是适用范围来看, 它是目前Matlab中最强大的数值积分函数。

3. 结论

本文对Matlab中的几种数值积分函数进行了简单阐述及说明, 对函数trapz、quad、quadl及quadgk进行具体分析及比较, 得出如下结论:

总的来说, 在Matlab的几种数值积分函数里, 函数quadgk的功能最为强大。

摘要：Matlab软件被广泛应用于科学计算, 而数值积分又是很重要的一个分支。本文首先对Matlab中的几种数值积分函数进行阐述和分析, 然后通过举例对函数trapz、quad、quadl及quadgk进行具体比较, 总结了各个数值积分函数的优缺点及适用性。

关键词：Matlab,数值积分函数,适用范围

参考文献

[1]张志涌, 等编著.精通MATLAB R2011a[M].北京:北京航空航天大学出版社, 2011, 11.

[2]华东师范大学数学系编著.数学分析上册[M].北京:高等教育出版社, 2011.

[3] (美) 里德 (Leader, J.J.) 著.数值分析与科学计算[M].张威, 等译.北京:清华大学出版社, 2008, 5.

[4]袁东锦编著.计算方法——数值分析[M].南京:南京师范大学出版社, 2004, 7.

Matlab数值分析 篇4

1 关于MATLAB软件

MATLAB是矩阵实验室Matrix Laboratory的简写, 它是一款数学软件, 主要用于数值运算, 它可以实现矩阵的运算、绘制函数和数据, 同时, 它还具有实现算法、创建用户界面的功能, 它的语法简单, 编程易于学习, 对现在的学生操作起来较容易, MATLAB软件强大的符号运算功能, 几乎包括了高等数学所涉及到所有的运算, 是当代理工科大学生应熟练掌握的一款数学软件。

2 关于数值分析课程的特点

数值分析是大学理工科开设的一门专业课, 学习内容涉及的范围较广, 主要有误差分析, 线性方程组的解法, 插值法和数据拟合, 数值微分与积分, 非线性方程及方程组的解法, 常微分方程数值解。③通过该课程的学习, 学生应该能够掌握一些常见数值方法的构造原理并且能达到在实际问题中的正确使用, 对某些现象能利用所学的知识给出较为合理的解释, 比如说关于方法的误差、收敛性、所分析的问题的性态等, 同时还要求学生具有一定的计算机编程能力, 能将本课程中常见的数值方法编写成程序且在计算机上运行, 得出正确的结果。

3 数值分析课程教学中存在的问题

从数值分析的教科书上可以看出来, 这门课程中的公式推导和算法比较多, 计算量相对比较大, 若教师整堂课都在讲公式推导, 没有加丝毫的实际演示过程, 这样的教学过程很难调动学生主动学习的积极性, 繁琐枯燥的数学公式推导只会导致学生失去对该课程的学习兴趣, 严重影响到课堂教学中学生与教师的互动。例如, 我们在讲三次样条插值时, 关于三次样条插值函数的计算非常的繁琐和复杂, 但若结合MAT-LAB软件, 学生还知道, 直接调用软件中的库函数也是可以实现的, 当然他们也可以通过自编程序实现插值过程。

4 基于MATLAB软件的数值分析课程的教学

4.1 用MATLAB软件实现具体的实例来进行数值分析课程的教学

在数值分析课程的教学中, 由于本门课程较强的应用性和教学内容的复杂性, 如果对于课程中遇到的较难理解的基本概念和主要方法, 教师都能够用尽可能多的具体实例来给予诠释, 那么, 这便已经是对数值分析课程教学改革的一个很好的尝试。因此, 结合数值分析这门课程理论教学的特点和存在问题, 教师只有改进传统的单一的教学模式, 才能更好地提高这门课程的教学效果, 从而能够使学生对所学内容有一个完整的理解和掌握。

对于实例教学, 如果我们使用MATLAB则很容易地可以实现大量复杂的数值计算和图形, 但如果用传统教学方式则难以解决同样的问题。例如在讲曲线拟合时, 就可以用人口增长问题作为实例。关于人口增长有两个基本的模型:指数增长模型和阻滞增长模型, 按照以上两种模型再结合线性最小二乘法, 我们就可以分别拟合出美国1790至1990年的人口数据, 还可以让学生观察下我们拟合出的结果与真实数据之间的差异。④又如三次样条插值问题, 常微分方程数值解问题, 数值积分问题等等, 我们都可以很容易地使用MATLAB软件提供的现有的函数直接得出想要的结果。⑤

4.2 用MATLAB软件实现教学内容的直观化形象化来进行数值分析课程的教学

现在, 随着计算机技术的迅速发展以及各类教学演示软件的开发, 在我们的课堂教学中, 可视化教学必会成为数值分析课程教学改革的一个方向, 因此, 对于本课程中问题的引入, 算法的讲解以及最终的结果, 教师尽可能都将其转化成图形等可视的结果展示给学生, 从而更形象, 减轻学生的学习压力, 来更好地激发他们的学习兴趣, 而MATLAB这款数学软件, 正具备有良好的可视化功能, 若将其合理地应用到我们的教学中, 必定会产生良好的教学效果。

对本课程中许多问题的求解过程, 在课堂教学中, 教师若能应用MATLAB软件直接进行演示, 就可以将枯燥抽象的数学内容, 繁杂庞大的计算过程非常直观地呈现出来, 从而使学生有着非常鲜明的感性认识。例如, 在讲利用雅可比迭代、高斯—赛德尔迭代和超松弛迭代法求线性方程组的数值解时, 教师直观形象的展示后, 要求学生比较这三种不同的方法, 它们各自的收敛速度如何以及收敛和发散的区别, 这些问题的提出, 都能够大大地激发学生的学习兴趣, 同时, 还可以逐渐培养学生的科学计算能力。

4.3 用MATLAB软件来进行数值分析课程的实验教学

实验是数值分析课程教学的一个重要组成部分, 是学生从理论知识到实际应用的一个重要过程, 利用计算机进行数学实验成为学生加深对所学知识的认识的一个重要手段, 但如果用C语言来实现, 就比较困难, 而且如果算法选择不好, 还可能得出错误的结果, 况且绘制图形本身也是C语言的一个难点, 而如果我们用MATLAB软件来实现就会很容易, 因为ATLAB具有强大的计算功能、图形处理和良好的交互界面等功能, 因此它是进行数值实验理想的工具, 这样, 对于C语言很复杂的算法, 利用MATLAB, 只需要直接调用其现有的函数, 加上几个简单的语句, 就可以得出其图形, 而且得出的图形美观、准确, 在用MATLAB软件进行数值实验的过程中, 通过数形结合学生不仅掌握了教学内容, 而且也会深刻体会到现代计算工具的魅力, 减轻学生学习压力的同时, 让学生自己更积极地来学习该门课程。

注释

11 韩旭里, 万中.数值分析与实验[M].北京:科学出版社, 2006.

224万中, 罗汉.加强开放性数学实验课程研究推动数学教育改革[J].大学教育科学, 2003.84 (4) :52-53.

33 丁丽娟, 程杞元.数值计算方法.北京:北京理工大学出版社, 2005.8.

Matlab数值分析 篇5

目前评估海浪载荷对浮式结构物的影响主要采用物理模型试验的方法, 然而这种方法费时费力, 而数值模拟具有参数设置灵活、计算结果精确等优点, 正逐步成为设计领域研究海浪载荷的重要手段。由于受到计算机硬件发展水平和波浪理论不成熟等的制约, 早期的海浪数值仿真主要以二维为主, 但是对于需要研究海浪和浮式结构物相互作用过程中产生的波浪折绕射、漩涡等现象时, 显然不能满足工程的需要。因此研究三维海浪数值模型, 实现对海浪现象更加真实准确的描述, 是海上复式结构物设计领域中研究海浪和结构物相互作用的必然发展方向。

本文利用谱分析的方法, 在MATLAB环境下对三维随机海浪模型进行了数值模拟仿真, 并给出了三维随机海浪波面图, 为浮式结构物设计中计算海浪载荷提供了参考。

1 Longuet-Higgins长峰波海浪模型

海浪现象十分复杂, 具有随机性和非线性的特征, 因此要建立精确的海浪模型是十分困难的。根据海浪谱和随机海浪理论, 可将实际海浪看成为不同频率、不同传播方向、不同波高和不同初始相位的正弦波叠加的结果[1]。Longuet-Higgins模型就是把长峰波海浪的波动用无数个随机余弦波的叠加来描述, 其波幅的表达式为:

其中ξai、ωi和εi分别表示第i个余弦波的波幅、频率和相位。

在实际应用中, Longuet-Higgins模型常用海浪的频谱来表达, 即:

ζ (t) =i=∑∞12SζΔωcos (ωi+εi)

式中Sζ (ωi) 为海浪的频谱。目前应用较多的海洋频谱有Pierson-Moscowitz谱 (PM谱) 、Neumann谱 (N谱) 、ITTC双参数波谱等, 其中PM谱应用最为广泛, 能很好地代表实际充分发育的随机海浪, 故本文采用PM谱。PM谱的表达式为:

其中v为离海面19.5m处的平均风速。

平均风速为10m/s、12m/s、15m/s的PM谱的仿真曲线如图1所示。从图1中可以看出, 海浪的PM谱为窄带谱, 能量集中于某一频段内, 随着风速的增加, 即高海情情况, PM谱峰值增大, 能量更为集中, 如当风速达到15m/s时 (六级海情) , 谱峰值为3.6m2∙s, 能量集中于0.3-1.7Hz频带内。因此对海浪进行数值仿真时, 可以选取能量集中的频段内的有限个谐波成分进行叠加, 可减少计算量, 从而加快仿真速度。

2 三维不规则短峰波随机海浪仿真

2.1 基于谱分析的三维不规则短峰波随机海浪模型

Longuet-Higgins模型假设海浪只沿一个固定方向传播, 且假设其波峰和波谷线相互平行且垂直于前进方向。实际上, 海面上的波浪不仅波高不同, 频率不同, 还会向各个方向传播。这些谐波除产生在主风向上的主浪向以外, 在主浪向两侧的±π2角度范围内都会有谐波的扩散, 这样的海浪称为三维不规则短峰波海浪。

与Longuet-Higgins模型不同, 三维不规则短峰波海浪是由振幅不等、频率不等、初始相位不同、传播方向不同的谐波多重求和得到。双叠加模型[2]为:

其中ζaij、ωi、ki、μj、εij分别为组成谐波的波幅、角频率、波数、方向角和随机初相位, εij为[0, π2]内的随机数, (ξ, η) 为波面上某点的坐标。

根据谐波的波幅和频谱的关系, 三维不规则短峰波海浪的波面方程可进一步写为:

其中Sζ (ωi, μj) 为短峰波海浪的方向波谱。

通常认为波能的方向和频率是无关的, 则短峰波海浪方向波谱可用两个独立函数的乘积来表示, 即:

其中Sζ (ω) 为海浪的PM谱, μ为海浪的方向角, ϕ (μ) 为波能的扩散函数。ITTC提出的建议性形式为:

波数ki表示谐波的传播方向上2π距离内波的个数, 于是有ki=2πλi, 其中λi为第i个谐波的波长。因为波长不易获得, 由线性波动理论可知, 在不考虑表面张力的情况下, 数值上的波数和角频率有如下关系, 即:

式中g为重力加速度, D为水深。

2.2 三维不规则短峰波随机海浪仿真

某浮式结构物的作业和生存海域的海洋条件为:风向30°, 风速13.8m/s, 五级浪, 为校验该结构物所能承受的波浪载荷, 需要对该海情下的海浪进行模拟。

选择100m×100m大小的海域, 利用上述理论进行三维随机波浪仿真。各参数设置为:仿真频段[0.3Hz, 1.7Hz], 海浪的方向角μ=30°, 风速v=13.8m/s, 时间t=10s。图2为该海域的模拟海浪图。

该仿真海浪中, 最大波高Hmax=3.32m, 根据标准浪级波高的参考值[3], 五级浪对应的波高范围为[) 2.5, 4.0 m, 最大波高Hmax位于允许的波高范围内, 说明利用海浪谱来模拟三维随机海浪能够得到比较精确的海浪波面图和波高值。进一步根据流体的势流理论就可以分析计算出该结构物受到的海浪载荷, 为校验结构物的结构强度提供了必要的基础。

3 结论

海上浮式结构物结构强度校验需要计算分析海浪载荷, 该文利用海浪谱分析的方法, 实现了在开阔海域主要由风力引起的海浪的模拟, 该仿真海浪的波面图和波高符合标准浪级波高的参考值。进一步利用流体的势流理论就可以分析计算出结构物受到的海浪载荷, 为进一步的结构强度校核提供基础。

参考文献

[1]俞聿修.随机波浪及其工程应用[J].大连理工大学出版社, 1992.

[2]Khatri S K.In the search of a coastal ocean wave model.OCEANS 97, MTS/IEEE Conference Proceedings, 1997, 213-218.

Matlab数值分析 篇6

关键词：土壤水分运移,RICHARDS模型,膜下滴灌,数值模拟

近年来,许多学者对滴灌条件下土壤中水分水平和垂直运移的变化过程进行研究表明,土壤中水分的运移过程一般都是抛物线型的二阶线性偏微分方程,其解法一般可分为解析解、半解析解及数值解法三种。实际中,水分在土壤中的运移过程是十分复杂的,只有对复杂的问题加以抽象、简化以及在很简单的边界条件下才能有解析解[1],而对复杂的问题和复杂的边界条件,只能借助于半解析解[2]或数值方法[3]求解。对于水分在非饱和土壤中的运移规律,只能采用数值拟合来进行计算。

文献[4]研究表明,由于空间和时间上的强烈变异性,采用直接大量的测量通常是不可行的,实验结果常存在明显的不确定性,土壤在水平方向和垂直方向的空间变异性也限制了直接测量的方法在实际工作中的应用。因此在更为便捷和精确的实验技术之前,需要寻求替代直接实验的其它方法来预测这一重要的水力特性。因此有学者利用MATLAB对土壤水分特征曲线进行数值模拟。文献[5]利用MATLAB对土壤水分特征曲线参数进行数值模拟,通过四参数模型计算,得到的模拟结果与实际结果拟合较好。可见MATLAB在土壤水分运移的数值模拟方面取得了比较满意的成果。本文利用MATLAB和RICHARDS数学模型,对滴灌土壤中水分的运移过程进行数值模拟,为研究非饱和土壤的水分运移问题提供了有效的手段。

1滴灌土壤水分运移的数学模型

Richards方程是描述水分运移的基本方程,是一个二阶非线性偏微分方程,其解析求解仅在极简单的边界条件和高度简化的土壤水分运移参数模型的情况下才有可能,绝大多数的条件下只能数值求解[6]。

对滴灌而言,滴灌条件下土壤中的水分运移为三维问题,且考虑边界条件是一个非常复杂的情况,通常采用RICHARDS的数值模型,其方程为:

(1)式中,θ为含水率;K(θ)为水平运移距离;D(θ)为垂直运移距离;X,Y,Z为滴灌土壤中水分运移的方向。

对滴灌土壤中水分水平运移和垂直运移过程进行数值模拟,可以采用二维的RICHARDS方程进行仿真。其方程为:

模型的上边界是一个动态边界,表现为土壤积水区半径的运移过程,因而,模型在上边界求解时需特别考虑土壤水分运移的上边界条件问题。

2土壤湿润区水分运移的Matlab数值模拟

2.1点源滴灌的数值拟合

利用Matlab中的PDE工具箱对点源滴灌土壤水分运移的过程进行模拟。

PDE中的标准抛物线方程为:

(3)式中,U为定义在Ψ上的未知函数;C、A、F为定义在Ψ上的函数;Ψ为求解域;F为源、汇项。

用到的边界条件是Dirichlet条件和Neumann条件:

式中,G、Q、H、R为定义在Ψ上的函数,可依赖于时间T,也可是空间坐标的函数;N为Ψ上的单位外法向矢量。

将PDE中的标准抛物型方程表达式与Richards方程相比较,得出抛物型方程中各参数与Richards方程中各参数之间的对应关系:

u=θ(x,z,t);c=D(θ);d=1;f=0;ZK(θ)在边界条件中反应。

对于滴灌土壤的水分运移过程,首先考虑上边界条件。对试验砂壤土的所有地表积水直径扩展过程进行拟合发现,积水区直径与时间变化的关系可用幂函数关系表达:

(6)式中,A、C>0;1>B>0。

计算域上边界条件是:原点(0,0)及0≤X≤0.5ρ(t)范围内为积水区,土壤含水率始终是饱和含水率θs。

可用Drichlet条件,其中,h=1,r=θs,θ=θS;0≤x≤0.5ρ(T),Z=0。

积水区以外是土壤湿润区,边界的法向上没有水分通量,采用Neumann条件,即式(7)

本模拟在积水区范围内采用Neumann条件。没有可靠的流量Q与积水区直径ρ(t)之间的关系,也没有关于积水区直径运移规律的成熟模型,所以,积水区直径ρ(t)用滴水点范围代替,即在滴水点范围内有:

底部边界条件(Z=Z)为Dirichlet条件;

右侧边界条件(X=X)均为Dirichlet条件;

底部和右侧边界条件对应式(4)中的H=1,R=0。

计算域左侧边界条件为Neumann条件;其对应式(5)中的P=0,G=0。

根据标准方程的要求,将K(θ)表达成时间T的函数。随着时间的增加,土壤含水率增加。而D(θ)同样也可表示成关于时间的函数。

根据试验结果对它们进行回归,得到砂壤土点源滴灌在Q=2L/H,V=20L条件下这两个参量的表达式:

式中,饱和含水率θs=46.7%,初始含水率θ0=3.18%(均为体积含水率)。

模拟结果如图1所示:

试验中观测到砂壤土点源滴灌在Q=2L/h,V=20L条件下的土壤最大水平运移距离为39.2cm,垂直运移距离为26.3cm。仿真模拟结果显示水平运移距离为37.5cm,垂直运移距离是25.2cm。所以,模拟结果基本与试验观测结果相吻合。

从模拟中可以知道滴水初期,土壤垂直运移速率比水平运移速率快,后期垂直运移速率慢下来,最终水平运移半径比垂直运移半径大。另外,滴灌土壤的最大水平运移半径不在地表,而在地表下(5-10)cm深度处,这与实测结果较为一致。但是,由于模拟中是用滴水点范围代替积水区直径ρ(t),简化了积水区运移对土壤水分运移的影响效果,模拟出的水平运移湿润半径比实测的数值略小。

2.2线源滴灌的数值模拟

用于线源滴灌土壤含水率分布数值模拟的模型与用于点源滴灌数值模拟的模型相同,都以Richards方程为基础。线源滴灌模拟中将滴水点设在上边界的两端点上。同样用Matlab中的PDE工具箱进行模拟。

线源滴水点在计算域上边界的两端,积水区范围内采用Neumann条件:

上边界条件(Z=0cm)为Neumann条件,其对应式(5)中的p=0,g=0。同时,在交汇区的上边界法向上水分梯度为零。

上边界中部也采用Neumann条件,即:

计算域两侧的边界条件(x=0cm和x=X=30cm)为Neumann条件,其对应式(5)中的p=0,g=0。

底部边界(Z=30cm)条件为Dirichlet条件,对应的式(4)中的h=1,r=0。

同样,根据标准方程的要求,将K(θ)表达成时间T的函数。随着时间的增加,土壤含水率增加,使土壤扩散率也增加。而D(θ)也可表示成关于时间的函数,其表达式的同上式(9)和式(10)。

对砂壤土线源滴灌在Q=2L/h,V=20L和滴头间距30cm条件下的土壤水分分布过程进行模拟,得到参量的表达式:

式中,饱和含水率θs=46.7%,初始含水率θ0=3.18%。

模拟结果如图2所示:

试验中观测到砂壤土线源滴灌在Q=2L/h,V=20L条件下模拟的土壤垂直运移距离为26.5cm,而实际的垂直运移距离为27.6cm,模拟结果与试验结果基本相吻合。

3结论

通过Matlab对滴灌土壤的水分运移进行数值模拟,得到以下结论:

(1)在点源滴灌水平和垂直运移过程中,水平和垂直的运移距离均随运移时间的增加而逐步增大,在滴水初期运移速率较大,随滴水时间的增加,水平运移速率逐步减小,垂直运移速率逐步增大。

(2)在线源滴灌水平和垂直运移过程中,在未形成交汇区时,土壤中水分的运移规律与点源中土壤水分运移的规律基本相同;当交汇区形成后,土壤水分的水平运移速率加快,含水率等值线变的宽浅。

参考文献

[1]雷廷武.滴灌湿润比的解析设计.水利学报,1994;25(1):1—9

[2]田富强,胡和平.基于常微分方程求解器的Richards方程数值模型.清华大学学报(自然科学版).2007;6(47):785—788

[3]池宝亮,黄学芳,张冬梅,等.点源地下滴灌土壤水分运动数值模拟及验证.农业工程学报,2005;21(3):56—59

[4] Van Genuchten MTh,Leij F J.On estimatiny the hydraulic propertiesof unsaturated soils.In:Van Genuchten MTh.et al.Proc Int Workshon Indirect Methods of Estimating the Hydraulic Properties of Unsatu-rated Soils,University of California,Riverside,California 1992;1—14

[5]彭建平,邵爱军.用Matlab确定土壤水分特征曲线参数.土壤,2007;39(3):433—438

Matlab数值分析 篇7

关键词：Matlab,地表滴灌,土壤,水分运动

0 引言

水分在土壤中的运动是影响作物生长的一个主要方面,对农业灌溉决策管理中实施土壤水分的动态监测、科学地指导灌溉、提高灌溉精度和农田水分利用效率及农作物产量有着重要的现实意义[1]。同时,地表滴灌具有显著的节水、增产以及提高作物品质的优势,因而成为农业工程学科的重要研究领域。近年来,随着非饱和多孔介质中水和溶质运移模拟的进一步开展,国内外的研究机构陆续开发了一系列数学模型软件,进行水、溶质、热运移等模型的模拟。例如,Li[2]等用HYDRUS模拟了地表滴灌条件下的水、氮在土壤中的运移规律;魏义长[3]、彭建平[4]等应用Matlab软件研究了土壤水分特征曲线参数;李晓斌[5]等人运用Matlab进行了土壤水分运移数值模拟。本文在前人研究的基础上,综合考虑作物根系吸水情况建立了地表滴灌条件下土壤水分运动方程,应用Matlab软件编程模拟了地表滴灌条件下土壤水分运动情况,并与田间实测数据进行直观比对。

1 土壤水分运动模型构建

1.1 土壤物理参数

土壤持水曲线是研究土壤水力学性质必不可少的,在已经建立的众多数学模型中,Van Genuchten模型是目前运用最广泛的模型,该模型适应性好,可以应用于不同耕作条件下的土壤水分分析[6]。Van Geneuchten[7]关于土壤水分特征曲线(θ-h曲线)的函数表达式为[7]

undefined

(1)

将θ(h)关系式代入Mualem[8]提出的K(θ,h)关系式[8],即

undefined (2)

经推导可得渗透系数K(h)的表达式为

K(h)=KsSundefined[1-(1-S1/me)m]2 (3)

其中,Se=(θ-θr)/(θs-θr),m =1-1/ n,n >1;h为土壤负压水头(cm);θ为土壤体积含水率(cm3/cm3);θs和θr分别为土壤饱和含水率和残余含水率(cm3/ cm3);Ks为土壤饱和导水率(cm3/cm3);α,m,n为拟合经验参数。上面方程中包含 5 个独立参数 θr,θs,α,n和 Ks,孔隙连通性参数 l 对大多数土壤来说可取为0.5。

1.2 土壤水分运动方程

地表滴灌条件下土壤水分运动为三维流动问题。假设各层土壤为均质、各向同性、骨架不变形的多孔介质,不考虑气相和温度对水分运动的影响,并假设地下滴灌点源条件下土壤水分运动为轴对称,则水分运动可简化为轴对称的二维问题来处理[9]。同时,综合考虑根系吸水因素,土壤水分运动方程可表示为[9]

undefined (4)

其中,θ为土壤体积含水率(cm3/ cm3);t为时间(h);x为水平坐标(cm);y为垂向坐标(cm);h为土壤负压水头(cm);K(h)为x,y方向的渗透系数(cm/h);S为根系吸水函数(h-1),代表单位时间内根系从单位体积土壤中吸收的水量。根据文献[10]用非线性回归分析方法可以得出根系吸水率S随含水量的变化经验公式为[10]

S(y,t)=E0(t)AB100θ(y,t) (5)

其中,E0(t)为一时段t内土壤入渗率或潜在蒸发强度;A=1.972 8×10-4;B=1.185 8;θ(y,t)为一时段t内沙壤土深度y处平均体积含水量。

2 定解条件

2.1 初始条件

假设土壤初始化含水率、负压水头及根系吸水项初始条件为

θi(y,0)=θ0i 0≤y≤80cm;i=1,…,5

hi(x,y,0)=h0i 0≤x≤40cm;

0≤y≤80cm;i=1,…,5 (6)

其中,θi,hi分别为第i层土壤的初始含水率(cm3/ cm3)和土壤负压水头(cm);θ0i和h0i分别为它们的初始值;i为土壤的层次。

2.2 边界条件

计算区域为一个长(垂向)80 cm,宽(径向)40 cm的圆柱形区域,假定整个模拟区域内的土壤质地是均一的,且各向同性。由于选取的模拟时段较长,因此忽略灌水过程中饱和区半径随时间的变化,假设饱和区半径为定值Rs。计算区域上边界饱和区在灌水过程中为随时间变化的通量边界,为第三类边界条件即

undefined (7)

其中,δ(t)为滴灌过程中进水边界的通量(cm/h),δ(t)=Q(t)/πRs2,Q(t)为滴头流量(cm3/h);Rs为灌水饱和区半径(cm);x为水平方向坐标(cm)。

模拟计算区域的左、右边界(x=0和x=40)处,假定为不透水边界,即为零通量边界,拟计算区域的下边界处(y=0)。考虑地下水埋深较大的情况,假定为自由出流边界条件,此时,水流边界条件为第一类边界条件,即

undefined (8)

undefined (9)

3 MATLAB数值计算与模拟

Matlab软件是一套可以实现数值分析、优化、统计、偏微分方程数值解等领域的计算和图形显示功能的软件包。其语言表达形式极其简单,不需像传统的算法语言那样进行编程,几乎与通常的数学表达形式相同,其函数的形式分类成库,使用时直接调用这些函数并赋予实际参数即可快速而准确地解决问题。本节应用Matlab对系统数学模型进行数值求解、绘图的过程如下。

3.1 Van Genuchten模型参数的求算

VG模型含有4个参数:a,n(其中n含有m),θr,θs,所以按4参数求解。以保定地区农田(粉砂壤土)实测土壤水吸力和相应含水率数据(如表1所示)为例进行数值求解。

先打开M-file创建窗口,输入:

function K1=fun1(n, ndata)

K1=n(1)+(n(2)-n(1))·/(1+(n(3)*ndata) ·∧n(4))·…(1-1·/n(4));

其中,Van Genuchten模型的4个参数θr,θs,α,n分别以n(1)、n(2)、n(3)、n(4)代表。将其保存为Fun1.m文件名。然后在命令窗口输入试验数据ndata(土壤吸力值数据)、ydata(对应的土壤体积含水率数据),并给出4个参数运算开始的初始值n0,执行lsqcurvefit函数命令,输入:

ndata= [27.2 53.04 62.56 69.36 81.6 95.2 108.8 126.48 159.12 197.2];

ydata= [0.4340.406 0.401 0.392 0.382 0.365

0.351 0.3350.312 0.287];

n0=[0.1,0.1,0.001,1];

[n, resnorm]=lsqcurvefit(@fun, n0, ndata, ydata)按回车后,即可得出4个参数的拟合结果及残差平方和为

n=0.1641 0.5652 0.0506 1.5003

resnorm=4·8701e-004

获得θr,θs,α,n参数后,即可建立粉砂壤土的土壤水分特征曲线如下

θ=0.164 1+(0.565 2-0.164)/(1+(0.050 6×h)^1.500 3)^(1-1/1.5003)

最后,还可运用Matlab软件的绘图函数plot绘制出土壤水吸力和土壤水分含量实测数据的散点图和拟合曲线。具体操作是在命令窗口输入:

n1=[0∶1∶197.2];

y1= 0.1641+(0.5652-0.1641) /(1+(0.0506*n1).^ 1.5003)·^(1-1/1.5003);

plot(ndata,ydata, ′k0′, n1, y1,′k-′)

回车后即可得到图1所示的结果。

这就是参数拟合并绘图显示出结果的整个过程。由此可以看出拟合效果非常好,拟合曲线几乎全部穿过了实测数据点。

3.2 渗透系数拟合经验公式确定

假定Ks土壤饱和导水率已经测得为0.073 9cm/min,将上面获得θr,θS,α,n参数,代入渗透系数K(h)的表达式(3)得

function K2=inline (h)

K2=0.073 9×(1+(0.05.6h)∧1.5003) ∧(1-1/1.5003)/2×(1-((1+(0.05.6h)∧1.5003) ∧(1-1/1.5003)∧1.5)2

经回归分析得砂壤土水分特征曲线表达式为

K2=0.073 9×h-3.5030

利用plot(h,K2)函数画图与“常水头法”测定结果比较,可以看出拟合效果非常好,如图2所示。

3.3 土壤水分运动方程求解

根据以上各方程确定地表滴灌土壤水分运动数学模型为

undefined

(10)

式(10)为二阶非线性偏微分方程,用数值法求解。取单位水平面积,高度为计算层厚度的沙柱体进行研究。将计算区域用两组相互正交的直线分成钜形网格,沿y方向将深度分为n(n=10)个单元,节点编号为i=0,1,2,…,n;空间步长为Δy(同前),将时间划分为m个时段,节点编号为j=0,1,2,…,m,时间步长Δt=1d。分别对各节点(含内节点和边界节点)用隐式差分格式列水量均衡方程式得

undefined

(11)

式(11)为非线性三对角方程组,迭代线性化后采用追赶法求解,求解过程含消元和回代两个过程。数值计算时,K(θ),K(h)由拟合经验公式给出其参数取值,参数拟合并绘图显示出结果如图3所示。

4 结论

结合Van Genuchten模型、土壤水分运动基本方程和根系吸水率经验公式,建立了地表滴灌条件下土壤水分运动二维数学模型,利用Matlab软件对模型进行模拟,与实测数据进行对比,主要结论如下:

1)运用Matlab非线性拟合函数求算土壤持水曲线Van Genuchten模型的4个参数,既快捷又准确,是切实可行的。

2)运用Matlab对地表滴灌条件下土壤水分运动数学模型进行模拟,拟合效果较好,该模型对研究土壤水分运动具有一定的参考意义。

3)Matlab数值分析软件在农业水土工程、土壤生态及相关学科的应用,必将有助于本学科的深入研究与发展。

参考文献

[1]王成武,戈振扬.滴灌条件下土壤渗流分形特征研究[J].昆明理工大学学报,2007,32(1):105-107.

[2]Li J,Zhang J,Ren L.Water and nitrogen distribution as af-fected by fertigation of ammonium nitrate from a point source[J].Irrig Sci,2003,22(1):19-30.

[3]魏义长,刘作新,康玲玲,等.土壤持水曲线van Genuchten模型求参的Matlab实现[J].土壤学报,2004,41(3):380-386.

[4]彭建平,邵爱军.基于Matlab方法确定VG模型参数[J].水文地质工程地质,2006(6):25-28.

[5]李晓斌,王海波,孙海燕,等.基于MATLAB的土壤水分运移数值模拟[J].科学技术与工程,2009,9(20):5 978-5 981.

[6]王小华,贾克力,刘景辉,等.Van Genuchten模型在土壤水分特征曲线拟合分析中的应用[J].干旱地区农业研究,2009,27(2):179-183.

[7]Vgenuchten M T H.A closed-form equation for predictingthe hydraulic conductivity of unsaturated soils[J].Soil Sci-ence Society of America Journal,1980,44:892–898.

[8]Mualem Y.Hydraulic conductivity of unsaturated porousmedia:generalized macroscopic approach[J].Water Re-sources Research,1978,14(2):325-334.

[9]Simunek J,Sejna M,Genuchten M T V.The hydrus-2D Soft-ware Package for Simulating the Two-Dimensional Move-ment of water,Heat,and Multiple Solutes in Variably-Satu-rated Media[M].[S.l.]:Salinity Laboratory,1999.