基于fluent的阻力计算.docx

上传人:b****6 文档编号:5888348 上传时间:2023-01-01 格式:DOCX 页数:34 大小:1.14MB
下载 相关 举报
基于fluent的阻力计算.docx_第1页
第1页 / 共34页
基于fluent的阻力计算.docx_第2页
第2页 / 共34页
基于fluent的阻力计算.docx_第3页
第3页 / 共34页
基于fluent的阻力计算.docx_第4页
第4页 / 共34页
基于fluent的阻力计算.docx_第5页
第5页 / 共34页
点击查看更多>>
下载资源
资源描述

基于fluent的阻力计算.docx

《基于fluent的阻力计算.docx》由会员分享,可在线阅读,更多相关《基于fluent的阻力计算.docx(34页珍藏版)》请在冰豆网上搜索。

基于fluent的阻力计算.docx

基于fluent的阻力计算

基于fluent的兴波阻力计算

本文主要研究内容

本文的工作主要涉及小型航行器在近水面航行时的绕流场及兴波模拟和阻力的数值模拟两个方面。

在阅读大量文献资料的基础上,通过分析、比较上述领域所采用的理论和方法,针对目前需要解决的问题,选择合理的方法加以有机地综合运用。

具体工作体现在以下几个方面:

1.本人利用FLUENT软件的前处理软件GAMBIT自主建立简单回转体潜器模型,利用FLUENT求解器进行计算,得出在不同潜深下潜器直线航行的绕流场、自由面形状及阻力系数的变化情况。

2.通过对比潜器在不同潜深情况下的阻力系数,论证了增加近水面小型航行器的深度可以有效降低阻力。

通过对模型型线的改动,为近水面小型航行器的型线设计提供了一定的参考。

通过改变附体形状和位置计算了附体对阻力的影响程度,为附体的优化设计提供了一定的依据。

计算模型

航行器粘性流场的数值计算理论

水动力计算数学模型的建立

根据流体运动时所遵循的物理定律,基于合理假设(连续介质假设)用定量的数学关系式表达其运动规律,这些表达式成为流体运动的数学模型,它们是对流体运动的一种定量模型化,称为流体运动控制方程组。

根据控制方程组,结合预先给定的初始条件和边界条件,就可以求解反映流体运动的变量值,从而实现对流体运动的数值模拟预报,形成分析报告。

基于连续介质假设的流体力学中流体运动必须满足要遵循的物理定律:

1)质量守恒定律

2)动量守恒定律

3)能量守恒定律

4)组分质量守恒方程

针对具体研究的问题,有选择的满足上述四个定律。

船体的粘性不可压缩绕流运动,如果不考虑水温对水物理性质的影响,水的密度和分子粘性系数都是常数,同时没有能量的转换,就仅仅需要满足质量守恒定律、动量守恒定律。

在满足这些定律下所建立的数学模型称为Navier-Stokes方程。

另外,自由液面的存在也需要建立合适的数学模型。

本文是利用FLUENT进行数值模拟,而软件里面关于自由液面模拟是用界面追踪方法的一种-流体体积法(VOF),基于该方法所建立的数学模型称为流体体积分数方程。

另外,高雷诺数下的水动力问题还需要考虑粘性不可压缩流体的湍流运动。

对于湍流运动的数值模拟一直是流体力学数值计算的一个难点。

直接数值模拟(DNS)目前还仅仅在院校中研究,而且也仅限于二维流体问题。

大涡模拟(LES)向工程应用的过渡似乎还没有完成,并且就高雷诺数问题而言,对计算机硬件要求很苛刻。

目前,从算法的可行性、硬件要求的可实现性、完成任务所消耗时间和人力等方面看,基于湍流模型的数值计算更为工程实际所接受。

本章将会对各种湍流模型加以介绍。

粘性不可压缩流体流动数学模型

连续方程

任何流动问题都必须满足质量守恒方程即:

连续方程。

根据连续介质假设,单位时间内流体微团的质量变化等于同时间间隔内进入微团的总净质量。

按照这一定律,连续方程数学表达式写为:

(2.1)

以上是在笛卡尔直角坐标系下表示,上面给出的是瞬态可压流体连续方程。

由于对于潜艇粘性流场介质的不可压缩,密度ρ为常数,引入散度算子,则方程(2.1)变成为:

(2.2)

式中:

速度矢量V={u,v,w}。

上式为粘性不可压缩流体运动的连续方程。

动量方程

动量守恒方程也是任何流动系统都必须满足的定律。

根据牛顿第二定律,流体微团中流体的动量对时间的变化等于微团所受外力之和,即:

(2.3)

(2.4)

(2.5)

式中,p代表流体微团所受的压力;τxx、τxy、τxz等是因分子粘性作用而产生的作用在流体微团表面上的粘性应力τ的分量;Fx、Fy、Fz表示直角坐标系下三个方向上流体微团的体积力分量,如果体积力只有重力,且Z竖直向上,则Fx=Fy=0,Fz=-ρg。

式(2.3)~(2.5)是对任何类型的流体(包括非牛顿流体)均成立的动量守恒方程。

本文研究的范围属于牛顿流体,故粘性应力τ与流体的变形率成比例,有:

(2.6)

式中,µ是动力粘度系数,λ是第二粘度,一般可取λ=−2/3,将(2.6)代入式(2.3)~(2.5)得到张量形式的动量守恒方程:

(2.7)

式(2.7)就是动量守恒方程。

方程(2.1)和(2.7)组成了控制粘性不可压缩流体运动的基本数学模型。

对于低雷诺数的层流运动,上述方程组已经可以确切描述流体运动。

但湍流流动以脉动的速度场为基本特征,各速度在时间和空间上变化很快,给流场的数值模拟带来很大困难。

再则,湍流是一种极度复杂的物理现象,包含无规律性,扩散性,三维涡旋波动及耗散。

在实际工程计算中要对湍流进行数值模拟代价十分高昂。

然而研究表明,大尺度涡在流体运动中起主要作用。

由此可见,若采用时间平均、集合平均或者其他人工处理方法略去小尺度运动,将小尺度运动模型化后代入大尺度中,从而替代求解原有瞬时控制方程,就会花费较小的计算代价获得较高精度的数值解。

以此为出发点,提出了将速度分解成平均值和脉动值,则瞬时速度分量u可以表达为:

(2.8)

将式(2.8)代入(2.1)和(2.7)再对时间积分就会得到下面的平均流方程。

(2.9)

(2.10)

(2.11)

方程(2.9)是时均形式的连续方程,方程(2.10)是时均形式的Navier-Stokes方程。

方程(2.11)为Reynolds应力。

由于式(2.7)采用的是Reynolds平均法,因此方程(2.10)被成为Reynolds平均Navier-Stokes方程(Reynolds-AveragedNavier-Stokes,简称RANS方程)。

有式(2.9)和(2.10)组成的方程组共有五个方程(RANS方程实际是3个)现在新增了6个Reynolds应力,再加上原来4个时均未知量,总共9个未知量,因此,方程组不封闭,必须引入新的湍流模型(方程)才能使方程组(2.9)和(2.10)封闭。

湍流模型

为了使雷诺平均N-S方程(RANS方程)封闭可解,要根据湍流的运动规律来寻求附加的条件和关系式,这就形成了不同的湍流模型。

在FLUENT计算软件中可以供选择的湍流模型有:

一方程模型Spalart-Allmaras(S-A)、两方程模型k−ε(包括Sk−ε、RNGk−ε)和k-ω(包括Sk-ω和SSTk-ω)以及雷诺应力模型(RSM),下面将本文所用到的四种湍流模型加以介绍。

标准k−ε模型(Sk−ε)

标准k−ε模型是典型的两方程模型,该模型是目前应用最广泛的湍流模型。

k和ε是两个基本未知量,与之相对应的输运方程为:

湍流动能k方程为:

(2.12)

湍流耗散率ε方程为:

(2.13)

式中的湍流涡粘度µt可表示为:

(2.14)

(其中=0.09,为一常数)

式中:

Gk式由于平均速度梯度引起的湍动能k的产生项,Gb是由于浮力引起湍动能k的产生项,YM代表可压缩湍流中脉动扩张的贡献,C1ε、C2ε和C3ε为经验常数,σk和σε分别是与湍动能k和耗散率ε对应的湍流普朗特数,Sk和Sε是用户定义的源项。

模型常数C1ε、C2ε、Cµ、σk、σε的取值为:

C1ε=1.44,C2ε=1.92,Cµ=0.09,σk=1.0,σε=1.3

RNGk−ε模型

RNGk−ε湍流模型是由Yakhot及Orzag提出的,该模型中的RNG是英文“renormalizationgroup”的缩写。

在RNGk−ε湍流模型中,通过在大尺度运动和修正后的粘度项体现小尺度的影响,而使这些尺度运动有系统地从控制方程中去除。

所得到的k方程和ε方程,与标准k−ε模型非常相似:

湍流动能k方程:

(2.15)

湍流耗散率ε方程:

(2.16)

与标准k−ε模型相比较发现,RNGk−ε模型的主要变化:

通过修正湍动粘度,考虑了平均流动中的旋转流动情况;

在ε方程中增加了一项,从而反映了主流的时均应变率Eij,这样,RNGk−ε模型中产生项不仅与流动情况有关,而且在同一问题中也还是空间坐标的函数。

模型常数C1ε,C2ε由RNG理论:

C1ε=1.42,C2ε=1.68,

其他常数:

Cµ=0.0845,σk=1.0,σε=1.3

Sk−ω模型

本文采用的Sk−ω湍流模型是基于湍流动能k和特殊湍流动能耗散率ω的输运方程建立起来的经验公式。

是由Wilcox在1998年提出的对原k−ω模型的改进模型。

湍流动能k方程为:

(2.17)

特殊耗散率ω方程为:

(2.18)

Γk、Γω表示k、ω的有效扩散率,表示为:

(2.19)

(2.20)

σk、σω分别为湍流动能k和湍流耗散率ω的普朗特数,湍流涡粘度µt可表示为:

(2.21)

α∗为低湍流雷诺数修正系数:

(2.22)

上式中α∗=βt/3,Ret为雷诺数:

(2.23)

以上各式中的常数取值为:

剪切应力输运k−ω模型(SSTk−ω)

SSTk−ω湍流模型由Menter提出,该模式的湍流动能方程和湍流耗散率方程与标准Sk−ω模型的形式相似:

湍流动能k方程为:

(2.24)

特殊耗散率ω方程为:

(2.25)

Γk、Γω和µt见式(2.18)~(2.20),常数σk、σω表示为:

(2.26)

(2.27)

式中F1是混合函数:

(2.28)

(2.29)

(2.30)

式中:

Gk式由于平均速度梯度引起的湍动能k的产生项,Gω是由于特殊湍流动能耗散率ω的产生,Dω为横向扩散项,Yk、Yω表示湍流k、ω的消耗,Sk和Sε是用户定义项。

边界条件

边界条件类型简介

流体在运动的过程中会受到边界的限制,反映到物理模型上,就是要给控制方程加一些关于变量Ui、P、k、ε相应的边界条件。

最常见的线性边界条件有两大类:

第一类边界条件(Dirichlet条件)和第二类边界条件(Neumann条件)。

前者描述的是计算区域的边界或部分边界上变量的值,后者则描述边界

上变量梯度的法向分量值,即:

Dirichlet条件:

φ=φb在边界上

Neumann条件:

ф=φn在边界上

式中φ为任意的物理量,

表示物体表面的单位外法线矢量,φb为给定的边界上的数值,φn为给定的

ф在边界上的法向分量。

对于潜器粘性绕流,入流边界是一种人工边界,它不由物体的性质决定,因而不是固定不变的,它需要取得离潜器表面足够远才能尽量地反映真实情况。

入口处边界条件属于Dirichlet条件:

其速度是预先给定的,一般是均匀来流条件,湍动能k和耗散率s也是预先给定的。

出流边界条件则是虚拟的,出流边界到艇尾的距离也要合理确定以消除对流场计算的影响。

对于粘性流动,在固壁边界(如艇体表面)须满足对速度和湍动能k的无

滑移边界条件,即:

u=v=w=0,k=0

然而在靠近壁面的区域,由于湍动能被强烈地耗减,耗散率达到最大值,在固壁上不易给出s的边界条件,因为它在壁面上不等于零。

在与壁面相邻的粘性子层中,由于粘性的影响,局部雷诺数变得很小。

由于前述k-ε模型是一种高雷诺数模型,因而对粘性子层不再适用,一般采取Launder和Spalding提出的壁函数方法来处理。

使用边界条件的注意事项

1.边界条件的组合

在CFD计算域内的流动是由边界条件驱动的。

从某种意义上说,求解实际问题的过程,就是将边界线或边界面上的数据,外推扩展到计算域内部的过程。

因此,提供符合物理实际且适合的边界条件是极其重要的,否则,求解过程将很难进行。

CFD模拟过程中迅速发散的一个最常见的原因就是边界条件选择的不合理。

例如,只给定进口边界和壁面边界,而没有给定出口边界,那么,将不可能得到计算域的稳定解,CFD将越计算越发散。

这样的边界条件组合显然是不合理的。

在使用出口边界时需要特别注意,该边界只在进入计算域的流动是以进口边界条件给定(如在进口给定速度和标量)时才使用,而且仅推荐在只有一个出口的计算域中使用。

物理上,出口压力控制着流体在多出口间的分流情况,因此,在出口给定压力值要比给定出口条件(梯度为零)合理。

将出口条件和一个或多个恒压边界结合使用是不允许的,因为零梯度的出口条件不能指定出口的流量,也不能指定出口的压力,这样将使问题不可解。

2.流动出口边界的位置选取

如果流动出口边界太靠近固体障碍物,流动可能尚未达到充分发展的状态(在流动方向上梯度为零),这将导致相当大的误差。

一般来讲,为了得到准确的结果,出口边界必须位于最后一个障碍物后10倍于障碍高度或更远的位置。

对于更高的精度要求,还要研究模拟结果对出口位于不同距离时的影响的敏感程度,以保证内部模拟不受出口位置选取的影响。

3.近壁面网格

在CFD模拟时,为了获得较高的精度,常需要加密计算网格,而另一方面,在近壁面处为快速得到解,就必须将k-ε模型与结合了准确经验数据的壁面函数法一起使用。

要保证壁面函数法有效,就需要使离壁面最近的以内节点位于湍流的对数律层之中,即Y+必须大于11.63(最好是在30~500之间)。

这就相当于给最靠近壁面的网格到壁面的距离设定了一个下限。

但是,在流动的任意位置都使上述要求得到保证常常不太可能,典型的例子就是包含回流的流动。

4.随时间变化的边界条件

这类边界条件是针对非稳定问题而言的。

就是说边界上的有关流动变量并不是一成不变,而是随着时间变化的。

对于这类边界条件,需要将边界条件离散成与时间步长相应的离散结果,然后存储起来,供计算到相应的时间步时调用。

这类边界条件一般是与初始条件一同给定的。

k和ε的计算公式

在入口、出口或远场边界流入流域的流动,FLUENT需要指定输运标量的值。

在某些情况下流动流入开始时,将边界处的所有湍流量指定为统一值是适当的。

在大多数湍流流动中,湍流的更高层次产生于边界层而不是流动边界进入流域的地方,因此这就导致了计算结果对流入边界值相对来说不敏感。

然而必须注意的是要保证边界值不是非物理边界。

非物理边界会导致你的解不准确或者不收敛。

对于外部流来说这一特点尤其突出,如果自由流的有效粘性系数具有非物理性的大值,边界层就会找不到了。

湍流强度I定义为相对于平均速度u_avg的脉动速度的均方根。

小于或等于1%的湍流强度通常被认为低强度湍流,大于10%被认为是高强度湍流。

从外界测量数据的入口边界,你可以很好的估计湍流强度。

例如:

如果你模拟风洞试验,自由流湍流强度通常可以从风洞指标中得到。

在现代低湍流风洞中自由流湍流强度通常低到0.05。

对于内部流动,入口的湍流强度完全依赖于上游流动的历史,如果上游流动没有完全发展者没有被扰动,你就可以使用低湍流强度。

如果流动完全发展,湍流强度可能就达到了百之几。

完全发展的管流的核心的湍流强度可以用下面的经验公式计算:

I=0.16×(雷诺数Re)-1/8

雷诺数Re=速度×当量直径×(密度÷粘度)

湍流尺度l是和携带湍流能量的大涡的尺度有关的物理量。

例如在完全发展的管流中,l被管道的尺寸所限制,因为大涡不能大于管道的尺寸。

L和管的物理尺寸之间的计算关系如下:

L=0.07×l

其中L为管道的相关尺寸。

因子0.07是基于完全发展湍流流动混合长度的最大值的。

在速度入流,压力边界条件中,需要输入如下两个物理量,其计算公式如下:

其中0.09是湍流模型中指定的经验常数(近似为0.09)。

为最大速度,I为湍流强度

,L为相当尺寸。

数值计算方法

网格生成

用CFD方法进行流场计算时,首先要将计算区域离散化,即划分网格。

网格是CFD模型的几何表达形式,也是模拟与分析的载体。

计算网格的好坏直接影响到数值计算的可行性、收敛性以及计算精度。

对于复杂的CFD问题,网格生成极为耗时,且极易出错,生成网格所需时间常常大于实际CFD计算的时间。

因此,有必要对网格生成方式给以足够的关注。

网格(grid)分为结构网格和非结构网格两大类。

把节点看成是控制体积的代表。

在离散过程中,将一个控制体积上的物理量定义并存储在该节点处。

若节点排列有序,即当给出了一个节点的编号后,立即可以得出其相邻节点的编号。

这种网格称之为结构网格(structuredgrid)。

结构网格是一种传统的网格形式,网格自身利用了几何体的规则形状。

近几年来,还出现了非结构网格(unstructuredgrid)。

非结构网格的节点以一种不规则的方式布置在流场中。

这种网格虽然生成过程比较复杂,但却有着极大的适应性,尤其对具有复杂边界的流场计算问题特别有效。

无论是结构网格还是非结构网格,都需要按下列过程生成网格:

1.建立几何模型。

几何模型是网格和边界的载体。

对于二维问题,几何模型是二维面;对于三维问题,几何模型是三维实体。

2.划分网格。

在所生成的几何模型上应用特定的网格类型、网格单元和网格密度对面或体进行划分,获得网格。

3.指定边界区域。

为模型的每个区域指定名称和类型,为后续给定模型的物理属性、边界条件和初始条件做好准备。

网格生成是一个“漫长而枯燥”的工作过程,经常需要进行大量的试验才能取得成功。

因此,出现了许多商品化的专业网格生成软件。

如GAMBIT、TGrid、GeoMesh、preBFC和ICEMCFD等。

此外,一些CFD或有限元结构分析软件,如ANSYS、I—DEAS、NASTRAN、PATRAN和ARIES等,也提供了专业化的网格生成工具。

这些软件或工具的使用方法大同小异,且各软件之间往往能够共享所生成的网格文件,例如FLUENT就可读取上述各软件所生成的网格。

有一点需要说明,由于网格生成涉及几何造型,特别是3D实体造型,因此,许多网格生成软件除自己提供几何建模功能外,还允许用户利用CAD软件(AutoCAD、Pro/ENGINEER)先生成几何模型,然后再导入到网格软件中进行网格划分。

因此,使用前处理软件,往往需要涉及CAD软件的造型功能。

方程离散

数值计算是将描述物理现象的偏微分方程在一定的网格系统内离散,用网格节点处的场变量值近似描述微分方程中各项所表示的数学关系,按一定的物理定律或数学原理构造与微分方程相关的离散代数方程组。

引入边界条件后求解离散代数方程组,得到各网格节点的场变量分布,用这一离散的场变量分布近似代替原微分方程的解析解。

当前求解流体流动和传热方程的数值计算方法比较多,如有限差分法(FiniteDifferenceMethod),有限元法(FiniteElementMethod)、有限体积法(FiniteVolumeMethod)、边界元法、特征线法、谱方法、有限分析法和格子类方法等。

每种数值计算方法有各自的特点和各自的适用范围,其中通用性比较好、应用比较广泛的是前

3种。

1有限差分法

有限差分法(FiniteDifferenceMethod,简称FDM)是数值解法中最经典的方法。

它是将求解域划分成差分格式,用有限个网格节点代替连续的求解域,然后将偏微分方程(控制方程)的导数用差商代替,推导出含有离散点上有限个未知数的差分方程组。

求差分方程组(代数方程组)的解,就是微分方程定解问题的数值近似解,这是一种直接将微分问题变成代数问题的近似数值解法。

这种方法发展较早,比较成熟,较多的用于求解双曲型和抛物型问题。

用它求解边界条件复杂、尤其是椭圆问题不如有限元或有限体积法方便。

2有限元法

有限元法(FiniteElementMethod,简称FEM)与有限差分法都是广泛应用的流体动力学数值计算方法。

有限元是将一个连续的求解域任意分成适当形状的许多微小单元,并与各小单元分片构造插值函数,然后根据极值原理(变分或加权余量法),将问题的控制方程转化为所有单元上的有限元方程,把总体的极值作为各单元极值之和,即将局部单元总体合成,形成嵌入了指定边界条件的代数方程组,求解该方程就得到各节点上待求的函数值。

有限元法的基础是极值原理和划分插值,它吸收了有限差分法中离散处理的内核,又采用了变分计算中选择逼近函数并对区域进行积分的合理方法,是这两类方法互相结合、取长补短发展的结果。

它具有很广泛的适应性,特别适用于几何及物理条件比较复杂的问题,而且便于程序的标准化。

对椭圆型方程问题又更好的适应性。

有限元法因求解速度较有限差分法和有限体积法慢,因此,在商用CFD软件中应用不普遍。

3有限体积法

有限体积法(FiniteVolumeMethod)又称控制体积法(ControlVolumeMethod,CVM)。

其基本思想是:

将计算区域划分为网格,并使每个网格点周围有一个互不重复的控制体积;将待解微分方程(控制方程)对每个控制体积积分,从而得出一组离散方程。

其中的未知量是网格点上因变量φ。

为了求出控制体积的积分,必须假定因变量φ在网格点之间的变化规律。

从积分区域的选取方式来看,有限体积法属于加权余量法中的子域法,从未知解的近似方法看来,有限体积法属于采用局部近似的离散方法。

简而之,子域法加离散,就是有限体积法的基本方法。

就离散方法而言,有限体积法可视作有限元法和有限差分法的中间物。

有限元法必须假定φ值在网格节点之间的变化规律(即插值函数),并将其作为近似解。

有限差分法只考虑网格点上φ的数值而不考虑φ值在网格点之间如何变化。

有限体积法只寻求φ值节点值,这与有限差分法相类似;但有限体积法在寻求控制体积的积分时,必须假定φ值在网格点之间的分布,这又与有限单元法相类似。

在有限体积法中,插值函数只用于计算控制体积的积分,得出离散方程之后,便可忘掉插值函数;如果需要的话,可以对微分方程中不同的项采取不同的插值函数。

综上所述,有限体积法是目前在流体流动和传热问题求解中最有效的数值计算方法。

有限体积法也称为控制容积积分法,是20世纪六七十年代逐步发展起来的一种主要用于求解流体流动和传热问题的数值计算方法。

有限体积法是在有限差分法的基础上发展起来的,同时它又吸收了有限元法的一些优点。

有限体积法与有限元法和有限差分法一样,要对求解域进行离散,将其分割成有限大小的离散格式。

在有限体积法中每一网格点按一定的方式形成一个包围该节点的控制容积矿,有限体积法的关键步骤是将控制微分方程式在控制容积内进行积分。

有限体积法获得的离散方程,物理上表示的是控制容积的通量平衡,方程中各项有明确的物理意义。

有限体积法区域离散的节点网格与进行积分的控制容积分立。

由于有限体积法的诸多优点,当今大多数CFD软件都采用它来离散求解,如PHOENICS、FLUENT、STAR-CD、NUMECA等。

SIMPLE算法

流场计算方法的本质是对离散后的控制方程组的求解。

目前各种商用CFD软件普遍采纳的算法是压力耦合方程组的半隐式方法(SIMPLE算法),它属于压力修正法的一种。

SIMPLE是英文Semi—ImplicitMethodforPressure-LinkedEquations的缩写。

SIMPLE方法由Patankar与Spalding于1972年提出,是~种主要用于求解不可压流场的数值方法,也可用于求解可压流动。

它的核心是采用“猜测一修正”的过程,在交错网格的基础上来计算压力场,从而达到求解动量方程(Navier—Stokes方程)的目的。

SIMPLE方法的基本思想是对于给定的压力场(它可以是假定的值或者是上一次迭代计算所得到的结果),求解离散形式的动量方程,得出速度场。

因为压力场是假定的或不精确的,由此得到的速度场一般不满足连续方程,所以,必须对给定的压力场加以修正。

修正的原则是与修正后的压力场相对

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 自然科学

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1