上海大学冶金工程专业本科生课程.ppt

上海大学冶金工程专业本科生课程 吴永全* 冶金过程数值模拟冶金过程数值模拟 Numerical Modelling of Numerical Modelling of Metallurgical ProcessingMetallurgical Processing 主讲：吴永全主讲：吴永全 上海大学，材料科学与工程学院，材料上海大学，材料科学与工程学院，材料 工程系，冶金工程教研室工程系，冶金工程教研室 本科生课程冶金过程数值模拟 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 授课内容 授课内容授课内容 绪论绪论数值求解方法数值求解方法冶金过程数值模拟冶金过程数值模拟数学描述数学描述绪论绪论数值求解方法数值求解方法冶金过程数值模拟冶金过程数值模拟数学描述数学描述 导热问题数值模拟导热问题数值模拟 对流与扩散数值模拟对流与扩散数值模拟 流场计算简介流场计算简介 导热问题数值模拟导热问题数值模拟 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 目录 一维稳态导热一维稳态导热 1 1 一维非稳态导热一维非稳态导热 2 2 二维导热二维导热 3 3 一维稳态导热一维稳态导热 1 1 三维导热三维导热 4 4 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维稳态导热 先考虑忽略对流换热的导热微分方程： 式中，为介质密度，kg/m3；cp为质 量横压热 容，J/(kg)；S为单 位体积内产 生的热量（内热源，产生为正，消耗为负 ），W/m3；为导热 系数， J/(m2ms)。 不同坐标系下的具体表达形式： 直角坐标直角坐标 柱坐标柱坐标 球坐标球坐标 其中，对于柱坐标系：x=rcos，y=rsin，z=z； 对于球坐标系：x=rsincos，y=rsincos，z=rcos。 控制方程控制方程 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维稳态导热 第一类边界条件：直接给定边界上的温度值或随时间变化的函数。 边界条件边界条件 第二类边界条件：直接给定边界上的温度的导数值或导数值随时间变化的函数。 第三类边界条件：直接给定边界上的温度的导数值与温度之间的函数。 对流换热对流换热 辐射换热辐射换热 对流辐射换热对流辐射换热 式中，TL, Tf, Tw分别为界面、流体和环境温度，K；h为表面换热系数，W/(m2K)；0为斯 芯藩玻尔兹曼常数，5.6710-8W/(m2K4)；s为表面黑度（发射率），其值取01。 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维稳态导热 一维导热问题的数值求解一维导热问题的数值求解第一步：求解区域离散化第一步：求解区域离散化 x=0 123i-1 i i+1nn+1 x=L x x/2x/2 WE P i-1 i i+1 we (x)w(x)e x 将求解区域n等分，得到(n+1)个节点，相 应地共得到(n+1)个控制容积，从左到右依 次编号。其中内节点共(n-1)个，从左到右 编号分别是2到n，相应的控制容积宽度为 x。边界节点共2个，均为半控制容积， 左边界节点编号为1，右边界节点编号为 n+1，相应控制容积宽 度为x/2。 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维稳态导热 一维导热问题的数值求解一维导热问题的数值求解第二步：控制方程离散化第二步：控制方程离散化 x=0 123i-1 i i+1nn+1 x=L x x/2x/2 WE P i-1 i i+1 we (x)w(x)e x 首先针对稳态导热问题。显然，稳态导热 中温度对时间的偏微分为零，于是 于是用差分格式代替微分格式 代入并整理 这就是内节点i的离散化方程，也称 内节点i的差分方程。 采用同样的方法对边界节点进行分 析，可以得到： 恒温 对流换热 绝热 辐射换热 非线性非线性 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维稳态导热 一维导热问题的数值求解一维导热问题的数值求解第三步：离散化方程的求解第三步：离散化方程的求解 x=0 123i-1 i i+1nn+1 x=L x x/2x/2 WE P i-1 i i+1 we (x)w(x)e x 根据区域离散化，得到(n+1)个节点，相应 地共有(n+1)个待求的未知温度，显然需要 建立(n+1)个方程才能得到确定的解。对于 (n-1)个内节点，我们构建的(n-1)个控制微 分方程对应。再加上左右两个边界节点的 差分方程就共有(n+1)个方程了。 联立这(n+1)个方程组成封闭的方程组，有 唯一解。 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维稳态导热 一维导热问题的数值求解一维导热问题的数值求解举例一：肋片稳态导热举例一：肋片稳态导热 图中一等截面直肋，处于温度为80的流体中。肋表面与 流体之间的表面传热系数为45W/(m2K)，肋基处温度300 。肋片端部绝热。肋片由铝合金制成，热导率为 110W/(m)。肋片的厚度0.01m，高度为H=0.1m。 肋片的一维导热模拟肋片的一维导热模拟 h, Tf h, Tf x H L=1 Tw x 分析：分析：由于肋片的横向毕渥数Bi=h/=0.0041=eps x=G*x0+f; %迭代公式 n=n+1; tol=norm(x-x0); x0=x; if(n=M) disp(Warning:迭代次数太多，可能不收敛！); return; end fprintf(%3in,n); fprintf(%3.4fn,x); fprintf(%3.5fnnn,tol); end %在前面具有gauseidel函数的基础上，在 Matlab中求解的指令集： A=1 -0.332 0 0 0 0 0 0 0 0;-0.498 1 -0.498 0 0 0 0 0 0 0;0 -0.498 1 -0.498 0 0 0 0 0 0;0 0 -0.498 1 -0.498 0 0 0 0 0;0 0 0 -0.498 1 - 0.498 0 0 0 0;0 0 0 0 -0.498 1 -0.498 0 0 0;0 0 0 0 0 -0.498 1 -0.498 0 0;0 0 0 0 0 0 - 0.498 1 -0.498 0;0 0 0 0 0 0 0 -0.498 1 - 0.498;0 0 0 0 0 0 0 0 -0.992 1; b=199.674;0.326;0.326;0.326;0.326;0.326; 0.326;0.326;0.326;0.649; x0=300;290;280;270;260;250;240;230;220 ;210; x,n=gauseidel(A,b,x0,1.0e-4) 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维稳态导热 线性方程组的求解线性方程组的求解TDMA TDMA法法(TriDiagonal Matrix Algorithm)(TriDiagonal Matrix Algorithm) 对于三对角线元素组成的线性方程组： 当i=1，由第一个方程得到： 当i=2，由第二个方程得到： 对任意i，有： 特殊地，对于TN有： 用上一用上一 步的步的P P i i 和和Q Q i i 代代 入到下入到下 一步求一步求 出下一出下一 步的步的 P P i i +1+1和和 Q Q i i +1+1 以上一以上一 步系数步系数 P P i i 和和Q Q i i 为前提为前提 ，因为，因为 T TN N = =Q Q N N ，直接，直接 往上代往上代 入即求入即求 出所有出所有 的的T T i i 追 赶 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维稳态导热 线性方程组的求解线性方程组的求解雅克比 雅克比(Jacobi)(Jacobi)法法 雅克比法是一种迭代法，所谓迭代就是按次序逐点地由老值（前次迭代的值，带 (k)上标的温度值）计算新值（后一次迭代的值，用(k+1)上标表示的温度值）， 当全部节点的新值都计算出来后，一起用新值替代老值，作为下一轮计算的老值 ，这是一种逐批更新的迭代方法。 迭代过程可以起步于整个温度场的初始试探值，在计入边界条件后，经过反复迭 代直到所有节点都满足 就得到问题的收敛解。 初始试探值对最终解的正确性没有影响。 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维稳态导热 线性方程组的求解线性方程组的求解高斯 高斯- -赛德尔赛德尔(Gauss-Seidel)(Gauss-Seidel)迭代法迭代法 该方法也是按次序逐点由老值计算新值，但与雅克比方法相比，它在任何节点只 要已经计算出来新值（方程右边带(k+1)上标的温度值）便立即取代老值（方程 右边带(k)上标的温度值）。所以逐点更新是高斯-赛德尔法的特点。 比如，如果计算次序是由右边界向左边界推进，则： 高斯-赛德尔法的优点是边界信息深入内部的速度比较快，因而收敛比较快。 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 目录 一维稳态导热一维稳态导热 1 1 一维非稳态导热一维非稳态导热 2 2 二维导热二维导热 3 3 一维稳态导热一维稳态导热 1 1 一维非稳态导热一维非稳态导热 2 2 三维导热三维导热 4 4 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维非稳态导热 常物性、无内热源的一维非稳态导热控制方程： 控制方程控制方程 其中，a=/(cp)为热扩散率。 区域离散化区域离散化 x i i+1 n-1n n+1 (i+1,n-1)(i+1,n)(i+1,n+1) (i,n) 对于一维非稳态导热问题，时间和空间坐标可用作平面 坐标系（如图）表示。任一节点(i,n)上的温度值Tni表示 空间第n个节点在i时刻的温度值。 控制方程离散化控制方程离散化 时间一阶微商采用向前差商： 空间二阶微商采用二阶中心差商： 代入并整理： 以以 x x为特征长度的傅里叶数，称为为特征长度的傅里叶数，称为 网格傅里叶数，用网格傅里叶数，用FoFo表示。表示。 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维非稳态导热 控制方程控制方程差分格式差分格式 差分格式图解差分方程 显示差分截断误差：O()+(x)2 稳定性条件：Fo=a/ (x)20 T=0oC t0 0 12345 6 因为采用显示格式，所以从稳定性考虑确定时间步长的极限值为 同时考虑计算精度，我们选取=2s。 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维非稳态导热 数值求解数值求解举例二：无限大薄板非稳态导热举例二：无限大薄板非稳态导热 一无限大薄板，初始状态为均匀温度200oC在某一时刻=0 ，板东侧面温度突然降到0oC，西侧面保持绝热。板厚 L=2cm，热导率为10W/(mK)，cp=10106J/(m3K)。采用 有限差分法的显示格式，选用一合理的时间步长，计算 =10s、20s等时刻板内的温度分布。 控制方程：初始条件：边界条件： Step 2: Step 2: 控制方程离散化控制方程离散化有限差分法有限差分法 对于内部节点2、3、4的显式差分格式为 L=2cm x=0. 4cm 绝热 t0 T=0oC t0 0 12345 6 对于节点1，采用元体平衡法时，因为西侧板面绝热，所以只有来自东面的传热 对于节点5，同样采用元体平衡法时，东西侧均有传热 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维非稳态导热 数值求解数值求解举例二：无限大薄板非稳态导热举例二：无限大薄板非稳态导热 一无限大薄板，初始状态为均匀温度200oC在某一时刻=0 ，板东侧面温度突然降到0oC，西侧面保持绝热。板厚 L=2cm，热导率为10W/(mK)，cp=10106J/(m3K)。采用 有限差分法的显示格式，选用一合理的时间步长，计算 =10s、20s等时刻板内的温度分布。 控制方程：初始条件：边界条件： Step 3: Step 3: 线性方程组的获得及其求解线性方程组的获得及其求解 将各种数据代入得到 L=2cm x=0. 4cm 绝热 t0 T=0oC t0 0 12345 6 同时代入T6=0，于是 显然，这种显示差分得到 的结果在计算过程中非常 方便，无需迭代，每次只 需要知道前一次的计算结 果，就能够代入直接计算 下一时刻的结果了。 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维非稳态导热 数值求解数值求解举例二：无限大薄板非稳态导热举例二：无限大薄板非稳态导热 步数时刻 s (0) x=0cm 1 x=0.2cm 2 x=0.6cm 3 x=1.0cm 4 x=1.4cm 5 x=1.8cm (6) x=2cm 00200200200200200200200 122002002002002001500 24200200200200193.75118.750 36200200200199.22185.1698.440 48200200199.9197.56176.0784.670 510199.99199.99199.62195.17167.3374.930 612199.94199.94199.11192.24159.2667.750 714199.84199.84198.36188.98151.9562.250 816199.65199.65197.37185.52145.3657.90 918199.37199.37196.17181.98139.4554.360 1020198.97198.97194.8178.44134.1351.40 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维非稳态导热 数值求解数值求解举例二：无限大薄板非稳态导热举例二：无限大薄板非稳态导热 一无限大薄板，初始状态为均匀温度200oC在某一时刻=0 ，板东侧面温度突然降到0oC，西侧面保持绝热。板厚 L=2cm，热导率为10W/(mK)，cp=10106J/(m3K)。采用 有限差分法的显示格式，选用一合理的时间步长，计算 =10s、20s等时刻板内的温度分布。 控制方程：初始条件：边界条件： Step 2: Step 2: 控制方程离散化控制方程离散化隐式差分法隐式差分法 对于内部节点2、3、4的隐式差分格式为 L=2cm x=0. 4cm 绝热 t0 T=0oC t0 0 12345 6 对于节点1，采用元体平衡法时，因为西侧板面绝热，所以只有来自东面的传热 对于节点5，同样采用元体平衡法时，东西侧均有传热 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维非稳态导热 数值求解数值求解举例二：无限大薄板非稳态导热举例二：无限大薄板非稳态导热 一无限大薄板，初始状态为均匀温度200oC在某一时刻=0 ，板东侧面温度突然降到0oC，西侧面保持绝热。板厚 L=2cm，热导率为10W/(mK)，cp=10106J/(m3K)。采用 有限差分法的显示格式，选用一合理的时间步长，计算 =10s、20s等时刻板内的温度分布。 控制方程：初始条件：边界条件： Step 3: Step 3: 线性方程组的获得及其求解线性方程组的获得及其求解 将各种数据代入得到 L=2cm x=0. 4cm 绝热 t0 T=0oC t0 0 12345 6 同时代入T6=0，于是 对于隐式差分，上一 时间层的结果直接作 为下一时间层的高斯 -赛德尔迭代的常数 向量b，进行计算。 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 一维非稳态导热 数值求解数值求解举例二：无限大薄板非稳态导热举例二：无限大薄板非稳态导热 步数时刻 s (0) x=0cm 1 x=0.2cm 2 x=0.6cm 3 x=1.0cm 4 x=1.4cm 5 x=1.8cm (6) x=2cm 00200200200200200200200 12199.9958199.9958199.9621199.6255196.2925163.29930 24199.9786199.9786199.8412198.7361190.5159136.08280 36199.9364199.9364199.5986197.3199183.7117115.67030 48199.8554199.8554199.2073195.4296176.5295100.1720 510199.7217199.7217198.6526193.1452169.363188.2490 612199.5226199.5226197.9296190.5532162.440678.94840 714199.247199.247197.0419187.7353155.884971.58840 816198.886198.886195.9985184.7633149.75265.67810 918198.4334198.4334194.8119181.6975144.057660.8620 1020197.8848197.8848193.4967178.587138.792856.88080 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 目录 一维稳态导热一维稳态导热 1 1 一维非稳态导热一维非稳态导热 2 2 二维导热二维导热 3 3 一维稳态导热一维稳态导热 1 1 一维非稳态导热一维非稳态导热 2 2 三维导热三维导热 4 4 二维导热二维导热 3 3 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 二维导热 1). 1). 区域离散化区域离散化 P N S EW n s we (x)w+(x)e- x (x)n-(x)s+ y 二维问题的网格划分如图所示。x方向的步长为x，y方向的步 长为y。控制容积的形状可以是正方形（均匀网格）、矩形、 三角形或者六边形等。网格的步长可以是均分，也可以不均分 ，是具体需求而定。 如果是非稳态问题 ，除了空间坐标的离散，还包括时间坐标 的离散，于是在原有的二维中需要增加一维，从而变为三维的 离散空间。 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 二维导热 2). 2). 控制方程的离散控制方程的离散有限差分法有限差分法 对于二维常物性控制方程（无热源）： 对时间的一阶微商用一阶向前差商代替，对坐标的二阶微商用二阶中心差商代替： 代入控制方程得到内部节点的显式差分方程： 如果是均匀离散x=y，则有： 上两式就是显式差分的一般性表达，很明显，每个时间层k+1的格点(i,j)温度都表达成了上 一时间层k的以该格点(i,j)为中心的5个格点温度的显函数。 该格式的稳定性条件是：Fo1/4。 改成全隐式格式该怎么办？改成全隐式格式该怎么办？ 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 二维导热 3). 3). 边界条件的离散边界条件的离散 节点差 分 方 程稳定性条件 内节点 显式Fo1/4 隐式无条件 平直绝 热边界 节点 显式Fo1/4 隐式无条件 平直对 流边界 节点 显式Fo1/(4+2Bi) 隐式无条件 平直辐 射边界 节点 显式Fo1/2 隐式无条件 二维非稳态导热问题节点差分方程（无内热源，x=y） 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 二维导热 4). 4). 举例三举例三 如图所示一长400mm，宽200mm的矩形板，导热率是常数，无内热源。设板的周围各边界 的温度给定，且不随时间变化，厚度方向的温度变化忽略。试求该矩形板内的温度分布。 L=400mm H=200mm x y Step 1: Step 1: 生成计算网格生成计算网格 将求解区域（矩形板）的x方向8等分（m=8），y方向4等分 （n=4），x=y=50mm，因各边界温度均已知或可求，故只 需要计算内节点的温度即可。 Step 2: Step 2: 控制方程离散化控制方程离散化 内节点的控制方程离散后可直接写成： 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 二维导热 4). 4). 举例三举例三 如图所示一长400mm，宽200mm的矩形板，导热率是常数，无内热源。设板的周围各边界 的温度给定，且不随时间变化，厚度方向的温度变化忽略。试求该矩形板内的温度分布。 L=400mm H=200mm x y Step 3: Step 3: 建立线性方程组求解建立线性方程组求解 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 二维导热 x y 050100150200250300350400 200-100100100100100100100- 150100114.486130.1506148.8958174.5715215.6815292.2935456.2718860 100100127.7935157.2205190.8611233.7085295.8611397.2206572.7936860 50100139.4676180.0769223.6196273.5405336.8339427.9341577.6819860 0-150200250300350400450- 采用高斯-斯德尔迭代法，初值温度均为100oC，收敛精度为1.010-4oC，经过39次迭代得到 的收敛结果如下： 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 目录 一维稳态导热一维稳态导热 1 1 一维非稳态导热一维非稳态导热 2 2 二维导热二维导热 3 3 一维稳态导热一维稳态导热 1 1 一维非稳态导热一维非稳态导热 2 2 三维导热三维导热 4 4 二维导热二维导热 3 3 三维导热三维导热 4 4 上海大学冶金工程专业本科生课程 吴永全* 冶金数值 导热问题数值模拟 三维导热 针对三维带热源的非稳态导热控制微分方程 该式在时间间隔到+，空间上沿x方向自w至e，y方向自s至n，z方向自b至t，对控制容 积P进行积分，并整理成通用形式的离散方程 式中 上海大学冶金工程专业本科生课程 吴永全* 下次课再见！下次课再见！ 冶金过程数值模拟 上海大学冶金工程专业本科生课程 吴永全 8KcNfQiUlXp#s%v)y0B3F6IaLdOgSjVmYq!t*w-z1D4G7JbMeQhTkWoZr$u(x+A2E5H9KcNfRiUlXp#s&v)y0C3F6IaLdPgSjVnYq!t*w-A1D4G8JbMeQhTlWoZr%u(x+B2E6H9KcOfRiUmXp!s&v)z0C3F7IaMdPgSkVnYq$t*w-A1D5G8JbNeQhTlWo#r%u(y+B2E6H9LcOfRjUmXp!s&w)z0C4F7IaMdPhSkVnZq$t*x-A2D5G8KbNeQiTlXo#r%v(y+B3E6H9LcOgRjUmYp!s&w)z1C4F7JaMdPhSkWnZq$u*x-A2D5H8KbNfQiTlXo#s%v(y0B3E6I9LdOgRjVmYp!t&w-z1C4G7JaMePhSkWnZr$u*x+A2D5H8KcNfQiUlXo#s%v)y0B3F6I9LdOgSjVmYq!t&w- z1D4G7JbMePhTkWoZr$u(x+A2E5H9KcNfRiUlXp#s&v)y0C3F6IaLdOgSjVnYq!t*w-z1D4G8JbMeQhTkWoZr%u(x+B2E5H9KcOfRiUmXp#s&v)z0C3F7IaLdPgSkVnYq$t*w-A1D5G8JbNeQhTlWoZr%u(y+B2E6H9KcOfRjUmXp!s&v)z0C4F7IaMdPgSkVnZq$t*x-A1D5G8KbNeQiTlWo#r%v(y+B3E6H9LcOgRjUmYp!s&w)z1C4F7JaMdPhSkVnZq$u*x-A2D5G8KbNfQiTlXo#r%v(y0B3E6I9LcOgRjVmYp!t&w)z1C4G7JaMePhSkWnZr$u*x+A2D5H8KcNfQiUlXo#s%v(y0B3F6I9LdOgRjVmYq!t&w- z1C4G7JbMePhTkWnZr$u(x+A2E5H8KcNfRiUlXp#s%v)y0C3F6IaLdOgSjVnYq!t*w-z1D4G7JbMeQhTkWoZr$u(x+B2E5H9KcNfRiUmXp#s&v)y0C3F7IaLdPgSjVnYq$t*w-A1D4G8JbNeQhTlWoZr%u(y+B2E6H9KcOfRjUmXp!w-z1D4G7JbMeQhTkWoZr$u(x+B2E5H9KcNfRiUmXp#s&v)y0C3F7IaLdPgSjVnYq$t*w-A1D4G8JbNeQhTlWoZr%u(x+B2E6H9KcOfRiUmXp!s&v)z0C3F7IaMdPgSkVnYq$t*x-A1D5G8JbNeQiTlWo#r%u(y+B3E6H9LcOfRjUmYp!s&w)z0C4F7JaMdPhSkVnZq$t*x- A2D5G8KbNeQiTlXo#r%v(y+B3E6I9LcOgRjUmYp!t&w)z1C4F7JaMePhSkWnZq$u*x+A2D5H8KbNfQiUlXo#s%v(y0B3E6I9LdOgRjVmYp!t&w-z1C4G7JaMePhTkWnZr$u*x+A2E5H8KcNfQiUlXp#s%v)y0B3F6IaLdOgSjVmYq!t*w-z1D4G7JbMePhTkWoZr$u(x+A2E5H9KcNfRiUlXp#s&v)y0C3F6IaLdPgSjVnYq!t*w-A1D4G8JbMeQhTlWoZr%u(x+B2E6H9KcOfRiUmXp!s&v)z0C3F7IaLdPgSkVnYq$t*w-A1D5G8JbNeQhTlWo#r%u(y+B2E6H9LcOfRjUmXp!s&w)z0C4F7IaMdPhSkVnZq$t*x- A2D5G8KbNeQiTlWo#r%v(y+B3E6H9LcOgRjUmYp!s&w)z1C4F7JaMdPhSkWnZq$u*x-A2D5H8KbNfQiTlXo#s%v(y0B3E6I9LdOgRjVmYp!t&w-z1C4G7JaMePhSkWnZr$u*x+A2D5H8KcNfQiUlXo#s%v)y0B3F6I9LdOgSjVmYq!t&w-z1D4G7JbMePhTkWoZr$u(x+A2E5H9KcNfRiUlXp#s%v)y0C3F6IaLdOgSjVnYq!t*w-z1D4G8JbMeQhTkWoZr%u(x+B2E5H9KcOfRiUmXp#s&v)z0C3F7IaLdPgSkVnYq$t*w-A1D4G8JbNeQhTlWoZr%u(y+B2E6H9KcOfRjUmXp!s&v)z0C4F7IaMdPgSkVnZq$t*x-A1D5GOfRiUmXp#s&v)z0C3F7IaLdPgSjVnYq$t*w- A1D4G8JbNeQhTlWoZr%u(y+B2E6H9KcOfRjUmXp!s&v)z0C4F7IaMdPgSkVnZq$t*x-A1D5G8KbNeQiTlWo#r%u(y+B3E6H9LcOfRjUmYp!s&w)z0C4F7JaMdPhSkVnZq$u*x-A2D5G8KbNfQiTlXo#r%v(y0B3E6I9LcOgRjVmYp!t&w)z1C4G7JaMePhSkWnZq$u*x+A2D5H8KbNfQiUlXo#s%v(y0B3F6I9LdOgRjVmYq!t&w-z1C4G7JbMePhTkWnZr$u(x+A2E5H8KcNfRiUlXp#s%v)y0B3F6IaLdOgSjVmYq!t*w-z1D4G7JbMeQhTkWoZr$u(x+B2E5H9KcNfRiUmXp#s&v)y0C3F7IaLdPgSjVnYq$t*w- A1D4G8JbMeQhTlWoZr%u(x+B2E6H9KcOfRiUmXp!s&v)z0C3F7IaMdPgSoZr$u(x+B2E5H9KcNfRiUmXp#s&v)y0C3F7IaLdPgSjVnYq!t*w-A1D4G8JbMeQhTlWoZr%u(x+B2E6H9KcOfRiUmXp!s&v)z0C3F7IaMdPgSkVnYq$t*x-A1D5G8JbNeQiTlWo#r%u(y+B2E6H9LcOfRjUmXp!s&w)z0C4F7IaMdPhSkVnZq$t*x-A2D5G8KbNeQiTlXo#r%v(y+B3E6I9LcOgRjUmYp!t&w)z1C4F7JaMePhSkWnZq$u*x-A2D5H8KbNfQiTlXo#s%v(y0B3E6I9LdOgRjVmYp!t&w-z1C4G7JaMePhTkWnZr$u*x+A2E5H8KcNfQiUlXp#s%v)y0B3F6I9LdOgSjVmYq!t&w- z1D4G7JbMePhTkWoZr$u(x+A2E5H9KcNfRiUlXp#s&v)y0C3F6IaLdPgSjVnYq!t*w-A1D4G8JbMeQhTkWoZr%u(x+B2E5H9KcOfRiUmXp#s&v)z0C3F7IaLdPgSkVnYq$t*w-A1D5G8JbNeQhTlWo#r%u(y+B2E6H9LcOfRjUmXp!s&w)z0C4F7IaMdPgSkVnZq$t*x-A1D5G8KbNeQiTlWo#r%v(y+B3E6H9LcOgRjUmYp!s&w)z1C4F7JaMdPhSkWnZq$u*x-A2D5H8KbNfQiTlXo#r%v(y0B3E6I9LcOgRjVmYp!t&w)z1C4G7JaMePhSkWnZr$u*x+A2D5H8KcNfQiUlXo#s%z1C4F7JaMdPhSkWnZq$u*x- A2D5G8KbNfQiTlXo#r%v(y0B3E6I9LcOgRjVmYp!t&w)z1C4G7JaMePhSkWnZr$u*x+A2D5H8KcNfQiUlXo#s%v)y0B3F6I9LdOgRjVmYq!t&w-z1C4G7JbMePhTkWnZr$u(x+A2E5H8KcNfRiUlXp#s%v)y0C3F6IaLdOgSjVnYq!t*w-z1D4G8JbMeQhTkWoZr$u(x+B2E5H9KcNfRiUmXp#s&v)y0C3F7IaLdPgSjVnYq$t*w-A1D4G8JbNeQhTlWoZr%u(y+B2E6H9KcOfRjUmXp!s&v)z0C4F7IaMdPgSkVnYq$t*x-A1D5G8JbNeQiTlWo#r%u(y+B3E6H9LcOfRjUmYp!s&w)z0C4F7JaMdPhSkVnZq$u*x- A2D5G8KbNfQiTlXo#r%v(y+B3E6I9LcOgRjU