Fluent中用户自定义函数应用举例.doc

第10章 应用举例本章包含了FLUENT中UDFs的应用例子。10.1 边界条件10.2源项10.3物理属性10.4反应速率（Reacting Rates）10.5 用户定义标量（User_Defined Scalars）10.1边界条件这部分包含了边界条件UDFs的两个应用。两个在FLUENT中都是作为解释式UDFs被执行的。10.1.1涡轮叶片的抛物线速度入口分布要考虑的涡轮叶片显示在Figure 10.1.1中。非结构化网格用于模拟叶片周围的流场。区域从底部周期性边界延伸到顶部周期性边界，左边是速度入口，右边是压力出口。Figure 10.1.1: The Grid for the Turbine Vane Example常数速度应用于入口的流场与抛物线速度应用于入口的流场作了比较。当采用分段线性分布的型线的应用是有效的对边界型线选择，多项式的详细说明只能通过用户定义函数来完成。常数速度应用于流场入口的结果显示在Figure 10.1.2和Figure 10.1.3中。当流动移动到涡轮叶片周围时初始常速度场被扭曲。Figure 10.1.2: Velocity Magnitude Contours for a Constant Inlet x VelocityFigure 10.1.3: Velocity Vectors for a Constant Inlet x Velocity现在入口速度将用以下型线描述：这里变量在人口中心是0.0，在顶部和底部其值分别延伸到。这样速度在入口中心为20m/sec，在边缘为0。UDF用于传入入口上的这个抛物线分布。C源代码（vprofile.c）显示如下。函数使用了Section 5.3中描述的Fluent提供的求解器函数。 /*/* vprofile.c */* UDF for specifying steady-state velocity profile boundary condition */*/#include udf.hDEFINE_PROFILE(inlet_x_velocity, thread, position) real xND_ND;/* this will hold the position vector */ real y; face_t f; begin_f_loop(f, thread) F_CENTROID(x,f,thread); y = x1; F_PROFILE(f, thread, position) = 20. - y*y/(.0745*.0745)*20.; end_f_loop(f, thread)函数，被命名为inlet_x_velocity,使用了DEFINE_PROFILE定义并且有两个自变量：thread 和position。Thread是一个指向面的thread的指针，position是一个整数，它是每个循环（loop）内为变量设置的数值标签。函数通过声名变量f作为face_t的数据类型。一维数组x和变量y被定义为real数据类型。循环宏用于循环区域内每个面来创建型线，或数据数组。每个循环内，F_CENTROIDS为带指标f的面输出面质心（数组x）的值，指标f在被thread指向的线（thread）上。存储在x1中的y坐标用于为变量y赋值，然后用于计算x速度。这个值接着被分配到F_PROFILE，使用整数position（在你从Velocity Inlet面板中选择的UDF作为x速度的边界条件的基础上通过求解器传递给它）来设置内存中x速度面值。为了在FLUENT中使用这个解释式UDF，你必须首先编译它。Define User-Defined Functions Interpreted.在Interpreted UDFs面板中，在Source File Name区域命名你的函数。如果必要，在CPP Command Name区域输入你的C预处理程序类型和在Stack Size下输入堆栈大小。打开Display Assembly Listing在你的控制台窗口中能看到汇编列表，当函数编译时。点击Compile然后Close面板。为了选择这个用户定义函数作为所选区域的速度边界条件，打开Velocity Inlet面板。在X-Velocity下拉列表下，选择udf inlet_x_velocity，这个名字在上面已经命名给函数。这个函数将被使用，而不是出现在X-Velocity区域的0值（在这个例子中）。点OK接受这个新边界条件并关闭面板。在求解运行到收敛之后，获得了一个修正的速度场如Figure 10.1.4 and 10.1.5所示。速度场在入口中心显示了最大值，到边缘上降为0。Figure 10.1.4: Velocity Magnitude Contours for a Parabolic Inlet x VelocityFigure 10.1.5: Velocity Vectors for a Parabolic Inlet x Velocity10.1.2管内流动的瞬态速度入口分布 这个例子中，瞬时的周期性的速度边界条件使用UDF应用于管子的入口。速度形式为： （10.1.1） 管子长1m，半径为0.2m。假设充满密度为1kg/m3和粘度为kg/m-s的空气。空气速度关于平衡值的波动振幅为5m/s，频率为10rad/s。速度分布UDF（unsteady.c）的源代码显示如下。 /*/* unsteady.c */* UDF for specifying a transient velocity profile boundary condition */*/#include udf.hDEFINE_PROFILE(unsteady_velocity, thread, position) face_t f; begin_f_loop(f, thread) real t = RP_Get_Real(flow-time); F_PROFILE(f, thread, position) = 20. + 5.0*sin(10.*t); end_f_loop(f, thread)函数，被命名为unsteady_velocity,使用DEFINE_PROFILE macro定义。The utility RP_Get_Real(flow-time)用于访问real流动时间，它被分配给变量t。（RP_Get_Real的详细内容见Section 6.9）。在你编译这个UDF之前，你必须在Solver面板中指定为非稳态流动计算。Define Models Solver.接下来你将打开Interpreted UDFs面板。在Source File Name下面的文本框内输入函数的名字，如果必要，还有你的C预处理程序的名字。打开Display Assembly Listing。点击Compile然后关闭面板。Define User-Defined Functions Interpreted.由UDF定义的正弦速度边界条件现在可选为入口区域的X-Velocity。在X-Velocity面板中，在X-Velocity区域右边的下拉列表中选择udf unsteady_velocity，并点击OK。时间步参数在Iterate面板中设置。Solve Iterate. 在这个例子中，使用了0.0314的Time Step Size，因此在入口速度一个完整的震荡周期内将完成20个时间步。在每个时间步内执行的迭代次数被限制为20。在求解的早期阶段这个迭代次数不足以获得收敛的流场，但是在通过几个时间步之后，很可能导致收敛的时间步。UDF Profile Update Interval被设置为1因此每次迭代速度都被更新。在完成了60个时间步之后（或3个周期），你可以检查跨越压力出口由振荡入口条件引起速度数量级。为了在计算阶段收集这个信息，在开始迭代之前打开Surface Monitors面板。Solve Monitors Surface.增加Surface Monitors的指标到1。这将使你定义monitor-1（如果需要，你可重命名，在Name下的文本输入框内）的参数。选择Plot，所以所选的数量将作为计算的收益被画出来。选择Print在控制台窗口来看所选数的变化值。选择Write，因而信息将被写入文件，它将被给一个名字为monitor-1.out。（如果你改变了monitor的名字，这个名字将用作输出文件的前缀）。在Every下面，你可选择Iteration, Time Step, or Flow Time。为了监视每个时间步的结果，你应选择Time Step选项。通过点击Define.你可在Define Surface Monitor面板中指定被监视的数量。在这个例子中，在Report Of下的下拉列表中选择Velocity. and Velocity Magnitude。报告的区域是pressure-outlet-5,，在Surfaces列表中选择它。Area-Weighted Average在Report Type的下拉列表中选择，Flow Time在X Axis的下拉列表中选择。一旦完成了你要求的迭代次数，监视器应出现在所选绘图窗口中。作为选择，你可通过打开File XY Plot面板读入文件。Plot File.你可通过在Files下面的文本框中输入它的名字来读取输出文件和选择这个文件并点击Add.，然后点击Plot，你可获得Figure10.1.6的图形显示。Figure 10.1.6: Average Velocity Magnitude at the Pressure Outlet这个图精细地显示了速度以期望的振幅值5m/s围绕平衡值20m/s上下振荡。10.2 源项这部分包含了源项UDF的应用。在FLUENT中它是作为解释式UDF被执行的。10.2.1给管道流动添加动量源项（Adding a Momentum Source to a Duct Flow）当使用UDF模拟源项时，理解函数调用的内容是很重要的。当你增加了一个源项，它在单元上执行全局的循环时，FLUENT将调用你的函数。你的函数应只计算源项并返回它到求解器。在这个例子中，动量源项将被添加到2D笛卡儿管道流动中。管道长4m, 宽2m,使用通过中心线的对称边界条件来模拟。液体金属（其属性列举在表10.2.1中）在290K温度下从左边以1mm/s的速度进入管道。在金属沿着管道前进了0.5m之后，它受到了冷壁面的冷却，壁面温度保持在280K。为了模拟金属的凝固，只要它的温度一降到288K以下，动量源项就应用于金属。动量源项是与速度的分量成比例的，并有相反的符号： S x = - C v x （10.2.1） 这里C是常数。当液体凝固时，它的运动将减小到0，模拟的是固体的形成。（在这个例子中，能量方程不能用来解决凝固潜热。速度场只能作为凝固区域指示来使用）。为了提高求解的稳定性和收敛性，求解器应线性化源项。为了能使求解器做到这点，你必须在你的UDF中以导数的形式指定源项和求解变量之间的依赖关系。源项 只依赖于求解变量。它关于的导数为： （10.2.2）下面的UDF指定了源项和它的导数： Table 10.2.1: Properties of the Liquid MetalPropertyValueDensity8000 kg/m 3Viscosity5.5 kg/m-sSpecific Heat680 J/kg-KThermal Conductivity30 W/m-K/*/* UDF that adds momentum source term and derivative to duct flow */ /*/#include udf.h#define CON20.0DEFINE_SOURCE(cell_x_source, cell, thread, dS, eqn) real source; if (C_T(cell,thread) 288K）流体的的粘度对液体有一个分子值（），而更冷区域（）的粘度有更大的值（1.0kg/m-s）。在中等温度范围内（），粘度在上面给出的两个值之间按线性分布： （10.3.1）这个模型的基础是假设液体冷却时很快地变为高粘性，它的速度降低，所以模拟的是凝固。这里，没有为能量场做修正以包含凝结潜热。UDF的C源代码显示如下。/*/* UDF for specifying a temperature-dependent viscosity property */*/#include udf.hDEFINE_PROPERTY(cell_viscosity, cell, thread) real mu_lam; real temp = C_T(cell, thread); if (temp 288.) mu_lam = 5.5e-3; else if (temp 286.) mu_lam = 143.2135 - 0.49725 * temp; else mu_lam = 1.; return mu_lam;函数，被命名为cell_viscosity，使用了DEFINE_PROPERT定义在单元上。引入了两个real变量：temp,C_T(cell,thread)的值，和mu_lam, 由函数计算的层流粘度。温度值被检测了，在它的下降范围的基础上计算了适当的mu_lam值。在函数的结尾，mu_lam的计算值返回到求解器。为了使用用户定义的属性，你将使用Materials面板。在Viscosity的下拉列表下，选择user-defined选项。一旦你选了这个选项，User-Defined Functions面板将打开，在其中你可选择合适的函数名字。这个例子中，只有一个是有效的，但在别的例子中，你可从中选择几个函数。（如果你需要编译多于一个的解释式UDF，这些函数应在编译前连接。详细内容见Section 7.2.1）。这个模拟的结果相似于Section 10.2.1中得到的。Figure 10.3.1显示由于应用用户定义函数而导致的粘度场。粘度在很窄的范围内很快地从常数值0.0055变化到1.0kg/m-s。速度场（Figure 10.3.2）证明响应于粘度的增加，液体流动减慢，这在预料之中。在这个模拟中，有一个大的“mushy”区域，在其中液体的运动逐渐减小。这是与第一个模拟的对比，在第一个模拟中应用了动量源项，在流体运动中观察到一个急剧的变化。Figure 10.3.1: Laminar Viscosity Generated by a User-Defined FunctionFigure 10.3.2: Contours of Velocity Magnitude Resulting from a User-Defined ViscosityFigure 10.3.3: Stream Function Contours Suggesting Solidification10.4反应速率（Reaction Rates）这部分包含了一个反应速率UDF的应用。它在FLUENT中以编译式的方式执行。10.4.1自定义体积反应速率（A Custom Volume Recation Rate）作为一个编译式UDF的例子，考虑为一个两种气态物质的简单系统自定义体积反应速率。两种物质分别被命名为species-a和spexies-b。反应速率是一个由下面的表达式给出的从pecies-a转变species-b的速率： （10.4.1）这里是species-a的质量分数，和是常数。2D（平面）区域由90度的弯管组成。管子宽16 inches，约114 inches长。一个6 inches厚的多孔区域覆盖在底部和右壁面，反应仅发生在多孔区域。管子内的物质有相同的属性。密度是1.0kg/m3, 粘度为。区域的轮廓显示在Figure 10.4.1中。多孔介质是从左边入口延伸到区域顶部压力出口的线的下方和右方的区域。Figure 10.4.1: The Outline of the 2D Duct通过左边入口，气态纯物质species-a以0.1m/s的速度进入。气体进入多孔介质顶部的开区域和多孔介质本身，这里在两个坐标方向中的每一个上都有5m-1的惯性阻力。层流流场（Figure 10.4.2）显示了大部分气体是从多孔区域转向开区域。Figure 10.4.2: Streamlines for the 2D Duct with a Porous Region流动方式通过显示在Figure 10.4.3中的矢量图得到进一步证实。多孔区域的流动比开区域的流动慢的多。Figure 10.4.3: Velocity Vectors for the 2D Duct with a Porous Region用于模拟发生在多孔区域反应的包含UDF的C源代码（rate.c）显示如下：/*/* rate.c */* UDF for specifying a reaction rate in a porous medium */*/#include udf.h#define K1 2.0e-2#define K2 5.DEFINE_VR_RATE(user_rate, c, t, r, mole_weight, species_mf, rate, rr_t) real s1 = species_mf0; real mw1 = mole_weight0; if (FLUID_THREAD_P(t) & THREAD_VAR(t).fluid.porous) *rate = K1*s1/pow(1.+K2*s1),2.0)/mw1; else *rate = 0.;函数，被命名为user_rate, 用DEFINE_VR_RATE为给定物质的质量分数定义在单元上。UDF执行测试来检查多孔区域，仅应用反应速率方程于多孔区域。宏FLUID_THREAD_P(t)用于确定一个单元thread是否是流体（而不是固体）thread。变量THREAD_VAR(t).fluid.porous用于检测流体单元thread是否是多孔区域。为了使用这个UDF，设置名字为librate的目录来控制库目录结构。然后，建立一个如Section 7.3.3描述的共享库。接着，启动FLUENT，或者在你的case文件读入或者set it up。打开在Makefile过程建立的库，通过在Compiled UDFs面板中指定librate为Library Name和点击Open连接这个库到FLUENT执行。Define User-Defined Functions Compiled.现在你可通过在User-Defined Function Hooks面板的Volume Reaction Rate Function的下拉列表中选择它来指定user_rate作为反应速率。Define User-Defined Function Hooks.初始化和运行计算。Species-a的质量分数收敛解显示在Figure 10.4.4中。穿过多孔区域的运动的气体在管子的水平部分逐渐转化成species-b。尽管多孔区域之外species-b的某些扩散被100% and 0% species-a区域之间的宽阔的过度层所暗示，但在流体区域没有反应发生。Figure 10.4.4: Mass Fraction for species-a Governed by a Reaction in a Porous Region10.5 用户定义标量（User-Defined Scalars）这部分包含了用于联合用户定义标量输运方程的UDF的例子。关于在FLUENT中模化用户定义标量输运方程的信息见第9章。10.5.1使用用户定义标量的后处理（Postprocessing Using User-Defined Scalars）下面是一个编译式UDF例子，该UDF计算温度四次幂的梯度，并存储其数量于拥护定义标量中。例如，被计算的温度梯度可随后用于绘制等值线图。尽管这个UDF的实际应用是可疑的，这里它的目的在于显示用于后处理的任意数量梯度的计算方法。/*/* UDF for computing the magnitude of the gradient of T4 */*/#include udf.h/* Define which user-defined scalars to use. */enum T4, MAG_GRAD_T4, N_REQUIRED_UDS;DEFINE_ADJUST(adjust_fcn, domain) Thread *t; cell_t c; face_t f; /* Make sure there are enough user-defined scalars. */ if (n_uds t0,SV_T) T = C_T(F_C0(f,t),t-t0); F_UDSI(f,t,T4) = pow(T,4.); end_f_loop (f,t) /* Fill second UDS with magnitude of gradient. */ thread_loop_c (t,domain) if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4) & NULL != T_STORAGE_R_NV(t,SV_UDSI_G(T4) begin_c_loop (c,t) C_UDSI(c,t,MAG_GRAD_T4) = NV_MAG(C_UDSI_G(c,t,T4); end_c_loop (c,t) thread_loop_f (t,domain) if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4) & NULL != T_STORAGE_R_NV(t-t0,SV_UDSI_G(T4) begin_f_loop (f,t) F_UDSI(f,t,MAG_GRAD_T4)=C_UDSI(F_C0(f,t),t-t0,MAG_GRAD_T4); end_f_loop (f,t) 条件语言if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4)用于检测用户定义标量T4的存储是否被分配，而NULL != T_STORAGE_R_NV(t,SV_UDSI_G(T4)检测用户定义标量T4的梯度的存储是否被分配。除了编译和连接这个UDF，如在第7章描述的，你还必须在FLUENT中求解用户定义标量输运方程。Define User-Defined Scalars.更详细的内容见Section 9.3。10.5.2实现FLUENT的P-1辐射模型（Implementing FLUENTs P-1 Radiation Model）这部分给出了一个例子示范了P1辐射模型是如何使用UDFs和用户定义标量输运方程来实现的。在P1模型中，入射辐射变量G，在区域内通过由扩散（diffusion）和源项组成的方程描述。入射变量G的输运方程由方程10.5.1给出。扩散系数，由方程10.5.2给出，源项由方程10.5.3给出。这些方程更详细的讨论在用户指南的Section 11.3.4中。 （10.5.1） （10.5.2） （10.5.3）如用户指南中所显示的，壁面上G的边界条件等于辐射壁面热流的负数（方程10.5.4），在这个方程中，是外法向量。辐射壁面热流由方程10.5.5给出。 （10.5.4） （10.5.5）根据壁面入射辐射指定的边界条件的这种形式是不适宜的。这种混合边界条件可通过先使用方程10.5.4和10.5.5求解，导致方程10.5.6来避免。那么，的表达式由方程10.5.7给出的辐射壁面热流来取代方程10.5.5给出的。 （10.5.6） （10.5.7）出现在方程10.5.6和10.5.7中的附加项和是估计方程10.5.4中入射辐射梯度的结果。在FLUENT中，被控制标量的梯度的分量垂直于单元边界（面），是主要和次要分量总和的估计。主要分量代表由单元质心定义方向上的梯度，次要分量是沿着分隔两个单元的面的方向。从这个信息中，面法向量可确定。梯度的次要分量可使用FLUENT宏BOUNDARY_SECONDARY_GRADIENT_SOURCE来建立。这个宏的使用首先要求定义单元几何信息，它可使用第二个宏BOUNDARY_FACE_GEOMETRY（见Section 5.4.5）容易的获得。你将在定义G的壁面边界条件的 UDF中看到这些宏的调用。为了完成P1模型的实现，辐射能量方程必须与热能方程耦合。这可通过修改源项和能量方程的壁面边界条件来完成。首先应考虑如何修改能量方程的源项。入射辐射梯度与辐射热流是成比例的。辐射热流局部的增加（或减小）可归结于能量方程通过吸收和散射机理的局部减小（或增加）。所以，辐射热流梯度是能量方程的（负的）源项。如用户指南中显示的，入射辐射能量方程 10.5.3的源项等于辐射热流的梯度，因此，它的负值指定了需要修改的能量方程的源项。现在，考虑如何修改壁面上能量的边界条件。局部地，能量从壁面传递默认的计算流体的唯一模式是传导。在包含辐射影响时，必须计算流体与壁面之间的辐射传热。（如果你使用了FLUENT的内建的P1模型，这项将被自动地计算）。宏DEFINE_HEAT_FLUX允许通过指定Section4.3.3中讨论的方程的系数来修改壁面边界条件以适应第二种传热模式。对壁面的净辐射热流已经由方程10.5.5给出。比较这个方程与Section 4.3.3中方程将为得出合适的系数。在这个例子中，P1模型的实现通过六个分离的UDFs来完成。它们都被包含在一个源文件中，它作为编译式UDF来执行。入射辐射G的单独的用户定义标量输运方程，使用DEFINE_DIFFUSIVITY UDF来定义方程10.5.2中的和DEFINE_SOURECE UDF来定义方程10.5.3中的源项。壁面上G的边界条件通过在DEFINE_PROFILE中指定方程 10.5.7的负值作为指定热流。 DEFINE_ADJUST UDF用于通知FLUENT检查已定义的（在求解器中）合适的用户定义标量数。最后，能量方程必须指定源项等于入射辐射方程中使用的源项的负值，而DEFING_HEAT_FLUX UDF用于改变能量方程壁面上的边界条件。在求解器中，必须至少激活一个用户定义标量输运方程。标量扩散率在Materials面板中为标量方程指定。标量源项和能量源项在边界条件面板中为流体区域指定。壁面上标量方程的边界条件在边界条件面板中为壁面区域指定。DEFINE_ADJUST and DEFINE_HEAT_FLUX函数在User-Defined Function Hooks面板中指定。/*/* Implementation of the P1 model using user-defined scalars */*/#include udf.h/* Define which user-defined scalars to use. */enum P1, N_REQUIRED_UDS;static real abs_coeff = 1.0;/* absorption coefficient */static real scat_coeff = 0.0;/* scattering coefficient */static real las_coeff = 0.0;/* linear-anisotropic */ /* scattering coefficient */static real epsilon_w = 1.0;/* wall emissivity */DEFINE_ADJUST(p1_adjust, domain) /* Make sure there are enough user defined-scalars. */ if (n_uds t0; /* Do nothing if areas arent computed yet or not next to fluid. */ if (!Data_Valid_P() | !FLUID_THREAD_P(t0) return; begin_f_loop (f,thread) cell_t c0 = F_C0(f,thread); BOUNDARY_FACE_GEOMETRY(f,thread,A,ds,es,A_by_es,dr0); At = NV_MAG(A); if (NULLP(T_STORAGE_R_NV(t0,SV_UDSI_G(P1)Gsource = 0.;/* if gradient not stored yet */ elseBOUNDARY_SECONDARY_GRADIENT_SOURCE(Gsource,SV_UDSI_G(P1), dG,es,A_by_es,1.); gamma0 = C_UDSI_DIFF(c0,t0,P1); alpha0 = A_by_es/ds; beta0 = Gsource/alpha0; aterm = alpha0*gamma0/At; Ibw = SIGMA_SBC*pow(WALL_TEMP_OUTER(f,thread),4.)/M_PI; /* Specify the radiative heat flux. */ F_PROFILE(f,thread,position) =aterm*Ew/(Ew + aterm)*(4.*M_PI*Ibw - C_UDSI(c0,t0,P1) + beta0); end_f_loop (f,thread)DEFINE_HEAT_FLUX(heat_f