




版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
空气动力学方程:动量方程的物理意义详解1空气动力学方程:动量方程的物理意义1.1基础概念1.1.1动量与动量守恒定律动量是物体运动状态的一种度量,定义为质量与速度的乘积。在流体动力学中,动量的概念被扩展到流体的每一小部分,即流体质点。动量守恒定律是物理学中的一个基本原理,它指出在一个系统中,如果没有外力作用,系统的总动量保持不变。在流体动力学中,这一原理被用于描述流体在管道、大气或任何其他流体系统中的运动。1.1.2流体动力学基本原理流体动力学是研究流体(液体和气体)的运动和静止状态的学科。其基本原理包括连续性方程、动量方程和能量方程。连续性方程描述了流体质量的守恒,动量方程描述了流体动量的守恒,而能量方程描述了流体能量的守恒。这些方程共同构成了流体动力学的基础,用于分析和预测流体的行为。1.2动量方程的物理意义动量方程,也称为纳维-斯托克斯方程(Navier-Stokesequations),是流体动力学中描述流体运动的关键方程。它基于牛顿第二定律,即力等于质量乘以加速度,来描述流体内部的力和流体的加速度之间的关系。动量方程可以写作:ρ其中,ρ是流体的密度,u是流体的速度向量,p是流体的压力,μ是流体的动力粘度,f是作用在流体上的外力向量。方程左边描述了流体的加速度,右边描述了作用在流体上的力。1.2.1动量方程的组成部分惯性项:ρ∂u∂t+u⋅压力梯度项:−∇粘性力项:μ∇外力项:f描述了作用在流体上的外力,如重力、电磁力等。1.2.2动量方程的应用动量方程在空气动力学中有着广泛的应用,例如:飞机翼型分析:通过求解动量方程,可以分析飞机翼型周围的气流分布,预测升力和阻力。风洞实验:在风洞实验中,动量方程用于模拟和分析不同风速下物体的空气动力学特性。大气运动:动量方程也是气象学中预测天气和大气运动的基础。1.2.3动量方程的数值求解动量方程通常是非线性的,解析求解非常困难,因此在实际应用中,常采用数值方法求解。以下是一个使用Python和NumPy库求解一维动量方程的简单示例:importnumpyasnp
importmatplotlib.pyplotasplt
#参数设置
rho=1.225#空气密度,单位:kg/m^3
mu=1.7894e-5#空气动力粘度,单位:Pa*s
dx=0.1#空间步长,单位:m
dt=0.01#时间步长,单位:s
L=1.0#计算域长度,单位:m
T=1.0#计算时间,单位:s
u0=1.0#初始速度,单位:m/s
#创建网格
x=np.arange(0,L,dx)
t=np.arange(0,T,dt)
u=np.zeros((len(t),len(x)))
u[0,:]=u0
#数值求解
forninrange(len(t)-1):
foriinrange(1,len(x)-1):
u[n+1,i]=u[n,i]-(u[n,i]*dt/dx)*(u[n,i]-u[n,i-1])+(mu*dt/dx**2)*(u[n,i+1]-2*u[n,i]+u[n,i-1])
#绘制结果
plt.plot(x,u[-1,:],label='FinalVelocity')
plt.legend()
plt.show()代码解释参数设置:定义了空气的密度、动力粘度、空间步长、时间步长、计算域长度、计算时间以及初始速度。网格创建:使用NumPy创建空间和时间网格。数值求解:通过迭代更新速度场,使用显式欧拉方法求解一维动量方程。结果可视化:使用Matplotlib库绘制最终的速度分布。这个示例仅用于说明动量方程的数值求解过程,实际应用中需要考虑更多的边界条件和物理效应。1.3结论动量方程是空气动力学和流体动力学的核心,它描述了流体动量的守恒,是理解和预测流体行为的关键。通过数值方法求解动量方程,可以解决复杂流体系统中的问题,为工程设计和科学研究提供重要支持。2动量方程的推导2.1牛顿第二定律在流体中的应用牛顿第二定律,即力等于质量乘以加速度,是动量方程推导的基础。在流体力学中,我们考虑的是流体微元的运动,而非刚体。对于流体微元,牛顿第二定律可以表述为:作用在流体微元上的总力等于流体微元质量与加速度的乘积。假设我们有一个流体微元,其体积为dV,密度为ρ,速度为v2.1.1表面力对于表面力,我们考虑流体微元的六个面。在x方向上,流体微元受到的力为:F其中,p是压力,τ是剪切应力,dAF2.1.2体积力体积力fx直接作用在流体微元的每单位体积上,因此在xF2.1.3质量与加速度流体微元的质量为ρdV,其在x方向上的加速度为∂vρ2.1.4动量方程将上述方程两边的dV约去,得到xρ类似地,可以得到y和z方向的动量方程,从而构成完整的动量方程组。2.2连续性方程与动量方程的关系连续性方程描述了流体质量的守恒,即流体在任意体积内的质量不会随时间改变,除非有流体流入或流出。对于不可压缩流体,连续性方程可以简化为:∂其中,v是流体的速度向量。对于可压缩流体,连续性方程则为:∂其中,J是质量流密度。连续性方程与动量方程的关系在于,它们都是基于流体微元的平衡条件推导出来的。连续性方程描述了质量的守恒,而动量方程描述了动量的守恒。在实际应用中,连续性方程和动量方程通常需要同时求解,以获得流体的完整运动状态。2.2.1示例:使用Python求解一维不可压缩流体的动量方程importnumpyasnp
fromegrateimportsolve_ivp
#定义参数
rho=1.225#空气密度,单位:kg/m^3
L=1.0#空间域长度,单位:m
N=100#空间离散点数
dx=L/(N-1)#空间步长
x=np.linspace(0,L,N)#空间坐标
t_span=(0,1)#时间跨度
t_eval=np.linspace(t_span[0],t_span[1],100)#时间坐标
#定义动量方程
defmomentum_eq(t,v):
dvdt=np.zeros_like(v)
dvdt[1:-1]=-1/rho*(np.gradient(v,dx)[1:-1]+np.gradient(v,dx,edge_order=2)[1:-1])
returndvdt
#定义初始条件
v0=np.zeros(N)
v0[1:-1]=1.0#初始速度为1m/s
#使用solve_ivp求解动量方程
sol=solve_ivp(momentum_eq,t_span,v0,t_eval=t_eval)
#输出结果
print("速度随时间变化:")
print(sol.y)在这个例子中,我们使用了Python的numpy和scipy库来求解一维不可压缩流体的动量方程。numpy用于数值计算,egrate.solve_ivp用于求解常微分方程。我们定义了流体的密度、空间域长度、空间离散点数等参数,并使用solve_ivp求解了动量方程。最后,我们输出了速度随时间的变化情况。2.2.2结论动量方程和连续性方程是流体力学中描述流体运动状态的基本方程。它们的推导基于牛顿第二定律和质量守恒原理,通过考虑流体微元的力平衡和质量平衡,可以得到描述流体运动的数学模型。在实际应用中,这些方程通常需要数值求解,如上述Python示例所示。3动量方程的物理意义3.1动量方程的各分量解释动量方程是流体力学中的基本方程之一,描述了流体在运动过程中动量守恒的原理。在空气动力学中,动量方程通常以Navier-Stokes方程的形式出现,它包含了流体的惯性力、压力梯度力、粘性力和外力等项。动量方程可以写为:ρ其中:-ρ是流体的密度。-u是流体的速度矢量。-p是流体的压力。-τ是应力张量,描述了流体内部的粘性力。-f是作用在流体上的外力,如重力。3.1.1惯性力项ρ这一项描述了流体的惯性,即流体保持其运动状态的倾向。第一部分∂u∂t3.1.2压力梯度力项−这一项表示了压力梯度力,即由于压力在空间上的变化而产生的力。在空气动力学中,这一项对于描述流体的压缩性和流体在物体表面的流动至关重要。3.1.3粘性力项∇粘性力项描述了流体内部由于粘性而产生的力。在空气动力学中,这一项对于理解边界层的形成和流体的湍流行为非常重要。3.1.4外力项f外力项包括所有作用在流体上的外部力,如重力、电磁力等。在空气动力学中,重力是常见的外力,尤其是在大气流动和飞行器设计中。3.2动量方程在空气动力学中的作用动量方程在空气动力学中扮演着核心角色,它帮助我们理解和预测流体在物体表面的流动行为,包括边界层的形成、分离点的确定、以及流体的湍流特性。通过求解动量方程,工程师和科学家能够设计更高效的飞行器,优化飞机的气动性能,减少阻力,提高飞行效率。3.2.1边界层分析边界层是指流体紧贴物体表面的一层薄薄的流体,其速度从物体表面的零逐渐增加到自由流的速度。动量方程在边界层分析中用于描述这一层流体的速度分布和厚度变化。3.2.2分离点预测分离点是指流体从物体表面开始分离的点,这通常发生在物体的后缘或流体速度减小到零的地方。分离点的预测对于理解物体的气动阻力至关重要,动量方程能够帮助我们确定这一关键点。3.2.3湍流模拟湍流是流体流动的一种复杂状态,其中流体的速度和压力在时间和空间上随机波动。动量方程在湍流模拟中用于描述流体的平均行为,以及湍流脉动对流体动量的影响。3.2.4示例:使用Python求解二维边界层方程下面是一个使用Python和SciPy库求解二维边界层方程的简单示例。我们将使用有限差分法来近似方程的导数。importnumpyasnp
fromscipy.sparseimportdiags
fromscipy.sparse.linalgimportspsolve
#定义网格参数
nx=100#网格点数
L=1.0#物体长度
x=np.linspace(0,L,nx)
dx=x[1]-x[0]
#定义流体参数
rho=1.225#空气密度
u_inf=1.0#自由流速度
nu=1.5e-5#空气动力粘度
#初始化速度和压力
u=np.zeros(nx)
v=np.zeros(nx)
p=np.zeros(nx)
#定义边界条件
u[0]=0.0#物体表面速度为零
v[0]=0.0#物体表面垂直速度为零
#构建系数矩阵
main_diag=np.ones(nx)*(1+nu/(rho*dx**2))
off_diag=-nu/(rho*dx**2)
A=diags([main_diag,off_diag,off_diag],[0,-1,1],shape=(nx,nx)).toarray()
#构建右侧向量
b=np.zeros(nx)
b[0]=-nu*u_inf/(rho*dx)
#求解速度分布
foriinrange(1,nx):
b[i]=-nu*(u[i-1]-u_inf)/(rho*dx**2)
u=spsolve(A,b)
#输出速度分布
print("速度分布:")
print(u)在这个示例中,我们使用了有限差分法来近似边界层方程中的导数。我们首先定义了网格参数和流体参数,然后初始化了速度和压力。接着,我们构建了系数矩阵和右侧向量,用于求解速度分布。最后,我们使用SciPy的spsolve函数求解了线性方程组,并输出了速度分布。3.2.5结论动量方程在空气动力学中是理解和预测流体流动行为的关键工具。通过分析动量方程的各分量,我们可以深入理解流体的惯性、压力梯度、粘性和外力如何影响流体的运动。此外,动量方程的求解对于边界层分析、分离点预测和湍流模拟等空气动力学问题的解决至关重要。4动量方程的应用4.1动量方程在翼型设计中的应用4.1.1理解动量方程动量方程是流体力学中的一个基本方程,它描述了流体在流动过程中动量守恒的原理。在空气动力学中,动量方程特别重要,因为它帮助我们理解翼型(机翼的横截面形状)如何通过改变流过其表面的空气动量来产生升力。4.1.2翼型设计中的动量方程在设计翼型时,工程师们利用动量方程来预测不同翼型形状在特定飞行条件下的性能。例如,通过调整翼型的弯度和厚度,可以改变流过翼型上表面和下表面的空气速度,从而影响升力的产生。示例:计算翼型上的升力假设我们有一个翼型,其上表面和下表面的空气速度分别为Vupper和Vlower,空气密度为动量方程简化后的形式可以表示为:C_L=2*(V_upper^2-V_lower^2)/(V^2)其中V是自由流速度,可以通过翼型的攻角和来流速度计算得到。代码示例#翼型升力系数计算示例
defcalculate_lift_coefficient(V_upper,V_lower,V):
"""
根据动量方程计算翼型的升力系数。
参数:
V_upper:float
翼型上表面的空气速度。
V_lower:float
翼型下表面的空气速度。
V:float
自由流速度。
返回:
C_L:float
翼型的升力系数。
"""
C_L=2*(V_upper**2-V_lower**2)/(V**2)
returnC_L
#示例数据
V_upper=150#翼型上表面速度,单位:m/s
V_lower=140#翼型下表面速度,单位:m/s
V=145#自由流速度,单位:m/s
#计算升力系数
C_L=calculate_lift_coefficient(V_upper,V_lower,V)
print(f"升力系数C_L:{C_L}")4.1.3解释在上述代码示例中,我们定义了一个函数calculate_lift_coefficient,它接受翼型上表面和下表面的空气速度以及自由流速度作为输入,根据动量方程计算并返回翼型的升力系数。通过调整这些输入值,工程师可以评估不同翼型设计在特定飞行条件下的升力性能。4.2动量方程在飞行器控制中的作用4.2.1动量方程与飞行器控制动量方程在飞行器控制中扮演着关键角色,尤其是在理解飞行器如何通过改变其周围空气的动量来实现转向、爬升或下降。例如,飞机的副翼、升降舵和方向舵都是通过改变局部空气流的动量来控制飞机姿态的。4.2.2示例:副翼操作对飞机转向的影响当飞机的副翼操作时,一侧的翼面会向下偏转,而另一侧的翼面会向上偏转。这会导致一侧的翼型产生更多的升力,而另一侧产生较少的升力,从而产生侧向力,使飞机转向。代码示例#副翼操作对飞机转向影响的计算示例
defcalculate_roll_moment(V_left,V_right,rho,S,c,delta_aileron):
"""
根据动量方程计算副翼操作对飞机产生的滚转力矩。
参数:
V_left:float
左侧翼型的空气速度。
V_right:float
右侧翼型的空气速度。
rho:float
空气密度。
S:float
翼型的参考面积。
c:float
翼型的弦长。
delta_aileron:float
副翼偏转角度。
返回:
M_roll:float
飞机的滚转力矩。
"""
#假设升力系数与速度的平方成正比,与副翼偏转角度成正比
C_L_left=2*(V_left**2)*delta_aileron
C_L_right=2*(V_right**2)*delta_aileron
#计算两侧翼型的升力
L_left=0.5*rho*V_left**2*S*C_L_left
L_right=0.5*rho*V_right**2*S*C_L_right
#计算滚转力矩
M_roll=(L_left-L_right)*c/2
returnM_roll
#示例数据
V_left=150#左侧翼型的空气速度,单位:m/s
V_right=145#右侧翼型的空气速度,单位:m/s
rho=1.225#空气密度,单位:kg/m^3
S=10#翼型的参考面积,单位:m^2
c=1#翼型的弦长,单位:m
delta_aileron=0.1#副翼偏转角度,单位:弧度
#计算滚转力矩
M_roll=calculate_roll_moment(V_left,V_right,rho,S,c,delta_aileron)
print(f"滚转力矩M_roll:{M_roll}N*m")4.2.3解释在飞行器控制的代码示例中,我们定义了一个函数calculate_roll_moment,它计算了副翼操作对飞机产生的滚转力矩。通过改变副翼的偏转角度和两侧翼型的空气速度,我们可以评估副翼操作对飞机转向效果的影响。这个计算是基于动量方程的原理,通过比较两侧翼型的升力差异来实现的。通过这两个应用示例,我们可以看到动量方程在空气动力学中的重要性,它不仅帮助我们设计更有效的翼型,还使我们能够精确控制飞行器的飞行姿态。5动量方程的数值解法5.1有限差分法求解动量方程5.1.1原理有限差分法是求解偏微分方程的一种数值方法,它通过将连续的偏微分方程离散化为一系列离散的代数方程来近似求解。在空气动力学中,动量方程描述了流体在空间和时间上的动量变化,通过有限差分法,我们可以将动量方程在时间和空间上进行离散,从而得到可以数值求解的方程组。5.1.2内容动量方程的一般形式为:ρ其中,ρ是流体密度,u是流体速度向量,p是压力,τ是应力张量,f是体积力向量。在有限差分法中,我们首先将空间和时间离散化。假设我们有一个一维的流体流动问题,空间上离散为一系列网格点,时间上离散为一系列时间步。对于空间上的一个网格点i和时间步n,动量方程可以被近似为:ρ这里,Δt和Δx分别是时间步长和空间步长,上标n和n+1表示当前和下一个时间步,下标i,i5.1.3示例假设我们有一个简单的流体流动问题,流体在x方向上流动,初始速度为u0=1,流体密度ρ=1,压力p=0,应力τimportnumpyasnp
#参数设置
rho=1.0#流体密度
dx=0.1#空间步长
dt=0.01#时间步长
L=1.0#空间长度
N=int(L/dx)#网格点数量
u=np.ones(N)#初始速度分布
#边界条件
u[0]=0.0#左边界速度为0
u[-1]=2.0#右边界速度为2
#时间迭代
forninrange(100):
un=u.copy()#保存上一步的速度分布
foriinrange(1,N-1):
u[i]=un[i]-un[i]*dt/dx*(un[i]-un[i-1])
#输出最终速度分布
print(u)这段代码使用了有限差分法中的显式欧拉法来求解一维动量方程。在每个时间步,它通过上一步的速度分布来更新当前步的速度分布。5.2有限体积法在动量方程中的应用5.2.1原理有限体积法是另一种求解偏微分方程的数值方法,它基于守恒定律,将计算域划分为一系列控制体积,然后在每个控制体积上应用守恒定律,得到控制体积的守恒方程。在空气动力学中,动量方程的有限体积形式可以更准确地处理流体的守恒性质,特别是在处理复杂的流体动力学问题时。5.2.2内容动量方程的有限体积形式可以写作:d这里,V和S分别是控制体积和控制体积的表面,uu是动量通量,u5.2.3示例考虑一个二维的流体流动问题,流体在x和y方向上流动,我们使用有限体积法求解x和y方向上的动量方程。importnumpyasnp
#参数设置
rho=1.0#流体密度
dx=0.1#x方向空间步长
dy=0.1#y方向空间步长
dt=0.01#时间步长
Lx=1.0#x方向空间长度
Ly=1.0#y方向空间长度
Nx=int(Lx/dx)#x方向网格点数量
Ny=int(Ly/dy)#y方向网格点数量
u=np.zeros((Ny,Nx))#x方向速度分布
v=np.zeros((Ny,Nx))#y方向速度分布
#边界条件
u[:,0]=0.0#左边界速度为0
u[:,-1]=2.0#右边界速度为2
v[0,:]=0.0#下边界速度为0
v[-1,:]=1.0#上边界速度为1
#时间迭代
forninrange(100):
un=u.copy()#保存上一步的x方向速度分布
vn=v.copy()#保存上一步的y方向速度分布
foriinrange(1,Nx-1):
forjinrange(1,Ny-1):
u[j,i]=un[j,i]-un[j,i]*dt/dx*(un[j,i]-un[j,i-1])-vn[j,i]*dt/dy*(un[j,i]-un[j-1,i])
v[j,i]=vn[j,i]-un[j,i]*dt/dx*(vn[j,i]-vn[j,i-1])-vn[j,i]*dt/dy*(vn[j,i]-vn[j-1,i])
#输出最终速度分布
print(u)
print(v)这段代码使用了有限体积法中的显式欧拉法来求解二维动量方程。在每个时间步,它通过上一步的速度分布来更新当前步的速度分布,同时处理了x和y方向的速度分量。以上就是关于动量方程的有限差分法和有限体积法的原理和内容,以及具体的代码示例。通过这些方法,我们可以有效地求解空气动力学中的动量方程,为流体动力学问题提供数值解。6高级主题6.1动量方程与能量方程的耦合6.1.1耦合原理在空气动力学中,动量方程和能量方程的耦合是理解流体动力学行为的关键。动量方程描述了流体运动的速度分布,而能量方程则关注流体的温度和压力变化。两者通过流体的密度、速度、温度和压力之间的关系紧密相连,共同决定了流体的动态特性。密度与速度的关系动量方程基于牛顿第二定律,表达为:ρ其中,ρ是流体密度,u是流体速度,p是压力,τ是应力张量,f是体积力。温度与压力的关系能量方程描述了流体的内能变化,可以表示为:ρ其中,Cp是定压比热,T是温度,k是热导率,ϕ6.1.2耦合示例考虑一个简单的二维不可压缩流体流动问题,使用Python和NumPy库进行数值模拟。我们将展示如何在计算动量方程的同时,考虑能量方程的影响。importnumpyasnp
importmatplotlib.pyplotasplt
#定义网格和时间步长
nx,ny=100,100
dx,dy=1.0,1.0
dt=0.01
#初始化速度、温度和压力
u=np.zeros((nx,ny))
v=np.zeros((nx,ny))
T=np.zeros((nx,ny))
p=np.zeros((nx,ny))
#定义物理参数
rho=1.0#密度
mu=0.1#动力粘度
k=0.1#热导率
Cp=1.0#定压比热
f=np.zeros((nx,ny))#体积力
#边界条件
u[:,0]=1.0#左边界速度为1
u[:,-1]=0.0#右边界速度为0
T[0,:]=100.0#下边界温度为100
T[-1,:]=200.0#上边界温度为200
#主循环
forninrange(1000):
#更新速度
u[1:-1,1:-1]=u[1:-1,1:-1]-dt*(1/dx*(p[1:-1,2:]-p[1:-1,:-2])-mu/dx**2*(u[2:,1:-1]-2*u[1:-1,1:-1]+u[:-2,1:-1]))
v[1:-1,1:-1]=v[1:-1,1:-1]-dt*(1/dy*(p[2:,1:-1]-p[:-2,1:-1])-mu/dy**2*(v[1:-1,2:]-2*v[1:-1,1:-1]+v[1:-1,:-2]))
#更新温度
T[1:-1,1:-1]=T[1:-1,1:-1]+dt*(k/dx**2*(T[2:,1:-1]-2*T[1:-1,1:-1]+T[:-2,1:-1])+k/dy**2*(T[1:-1,2:]-2*T[1:-1,1:-1]+T[1:-1,:-2])+rho*Cp*f[1:-1,1:-1])
#更新压力
#这里省略了压力更新的复杂步骤,通常需要求解泊松方程
#可视化结果
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.imshow(u,cmap='viridis')
plt.colorbar()
plt.title('速度分布')
plt.subplot(1,2,2)
plt.imshow(T,cmap='viridis')
plt.colorbar()
plt.title('温度分布')
p
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- PET-cGMP-生命科学试剂-MCE
- Lobatoside-C-生命科学试剂-MCE
- Cy3-PEG-FA-MW-2000-生命科学试剂-MCE
- 跑道维修合同范本
- 试桩合同范本
- 品牌木门合同范本
- 2025年运动场馆灯具项目发展计划
- 诉讼保全担保合同(2篇)
- 教学设备整改报告范文
- 建筑学术交流报告范文
- 企业名称预先核准通知书
- 统筹管理方案
- 建筑工程安全文明施工标准化图集(附图丰富)
- 人教版 美术二年级上册 第9课 蜻蜓飞飞 教案
- Unit 1 Travel教案-2023-2024学年高一下学期 中职英语高教版(2023修订版)基础模块2
- DB3206T 1083-2024机关会议服务人员操作技术规范
- 习作《我的家人 》教案-2024-2025学年统编版语文四年级上册
- 眼镜学智慧树知到答案2024年温州医科大学
- 24年注安-其他安全-必背考点-王培山
- 《积极心理学(第3版)》 课件 第9章 积极的爱
- 中医淋巴排毒
评论
0/150
提交评论