高等数学笔记-乐经良老师-第九章-重积分

高等数学笔记-乐经良老师

第九章 重积分

第一节 二重积分的概念和性质

一、典型例子

01 平面薄板的质量

平面薄板位于 $xy$ 平面区域 $D$，其面密度为 $\mu (x,y)$ 如何求其质量？

类似一元的处理方法，采用：

(1) 分割

将 $D$ 任意划分成 $n$ 个小区域 $\Delta D_1,\Delta D_2,\cdots ,\Delta D_n,\Delta D_i$ 的面积记为 $\Delta {\sigma }_i,(i=1,2,\cdots ,n)$

(2) 作和

在小区域分得很小时，近似认为质量均匀，任取 $\left({\xi }_i,{\eta }_i\right)\in \Delta D_i$，薄板的质量近似地表达为
$$m=\sum _{i=1}^n\Delta m_i\approx \sum _{i=1}^n\mu \left({\xi }_i,{\eta }_i\right)\Delta {\sigma }_i$$
(3) 取极限

记 λ = max ⁡ 1 ≤ 1 ≤ n { d i } \lambda=\max \limits_{1 \leq 1 \leq n}\left\{d_{i}\right\} λ=1≤1≤nmax​{di​}，( $d_i$ 是小区域 $\Delta D_i$ 的直径 ) 那么若
$$m=\underset{\lambda \to 0}{\lim }\sum _{i=1}^n\mu \left({\xi }_i,{\eta }_i\right)\Delta {\sigma }_i$$
存在，就给出了薄板的质量。

02 曲顶柱体的体积

柱体的侧面是母线垂直 $xy$ 平面的柱面，顶面为曲面 $S:z=f(x,y)$，

底面是 $xy$ 平面上区域 $D$，如何求此曲顶柱体的体积？

(1) 分割

用曲线将 $D$ 分成小区域 $\Delta D_1,\Delta D_2,\cdots ,\Delta D_n$，而 $\Delta D_i$ 的面积记为 $\Delta {\sigma }_i$

(2) 求和

区域分得很小时，用柱体来近似小曲顶柱体的体积，任取 $\left({\xi }_i,{\eta }_i\right)\in \Delta D_i$，则总体积近似为
$$\sum _{i=1}^{n}f\left({\xi }_i,{\eta }_i\right)\Delta {\sigma }_i$$
(3) 取极限

记 λ = max ⁡ 1 ≤ 1 ≤ n { d i } \lambda=\max \limits_{1 \leq 1 \leq n}\left\{d_{i}\right\} λ=1≤1≤nmax​{di​}，( $d_i$ 是小区域 $\Delta D_i$ 的直径 ) 则体积 $V$ 由如下极限给出
$$V=\underset{\lambda \to 0}{\lim }\sum _{i=1}^{n}f\left({\xi }_i,{\eta }_i\right)\Delta {\sigma }_i$$
从以上例子抽象出来就得到二重积分的概念，这类问题要计算在一个平面区域上分布率不均匀的量的总量。

二、二重积分定义

01 二重积分的定义

设 $D$ 是 $xy$ 平面的有界闭区域，函数 $f(x,y)$ 在 $D$ 定义，$I$ 为实数，

若将 $D$ 任意划分成个小区域 $\Delta D_1,\Delta D_2,\cdots ,\Delta D_n$，

任取 $\left({\xi }_i,{\eta }_i\right)\in \Delta D_i,(i=1,2,\cdots ,n)$，作和
$$\sum _{i=1}^{n}f\left({\xi }_i,{\eta }_i\right)\Delta {\sigma }_i(\Delta {\sigma }_i表示\Delta D_i的面积)$$
总有
$$\underset{\lambda \to 0}{\lim }\sum _{i=1}^{n}f\left({\xi }_i,{\eta }_i\right)\Delta {\sigma }_i=I$$
( 其中 λ = max ⁡ 1 ≤ i ≤ n { d i } \lambda=\max \limits_{1 \leq i \leq n}\left\{d_{i}\right\} λ=1≤i≤nmax​{di​}，$d_i$ 是小区域 $\Delta D_i$ 的直径 )

则称函数 $f(x,y)$ 在 $D$ 上可积，$I$ 称为$f(x,y)$ 在 $D$ 的二重积分，记为 $\underset{D}{\iint }f(x,y)d\sigma$ .

其中，$\iint -$ 积分号，$D-$ 积分区域，$f(x,y)-$被积函数， $d\sigma -$面积元素。

02 二重积分的几何意义

以区域 $D$ 为底，以曲面 $S:z=f(x,y)$ 为顶的曲顶柱体的体积。

03 可积的充分条件

若函数 $f(x,y)$ 在有界区域 $D$ 上分片连续有界，则 $f(x,y)$ 在 $D$ 可积。

三、二重积分的性质

设以下性质中出现的积分均存在

• 性质1 (线性) ：若 $\alpha ,\beta$ 是常数，
  $$\underset{D}{\iint }(\alpha f(x,y)+\beta g(x,y))d\sigma =\alpha \underset{D}{\iint }f(x,y)d\sigma +\beta \underset{D}{\iint }g(x,y)d\sigma$$

• 性质2 (可加性) ：若积分区域 $D$ 分成 $D_1,D_2$ 两个子区域，
  $$\underset{D}{\iint }f(x,y)d\sigma =\underset{D_1}{\iint }f(x,y)d\sigma +\underset{D_2}{\iint }f(x,y)d\sigma$$

• 性质3：
  $$\underset{D}{\iint }1d\sigma =A_D(D的面积)$$

• 性质4 (单调性) ：若 $f(x,y)\le g(x,y)$，则
  $$\underset{D}{\iint }f(x,y)d\sigma \le \underset{D}{\iint }g(x,y)d\sigma$$
  性质4的推论：
  $$\begin{array}{rl}& (1)若f(x,y)\ge 0,则\underset{D}{\iint }f(x,y)d\sigma \ge 0\\ & (2)\left|\underset{D}{\iint }f(x,y)d\sigma \right|\le \underset{D}{\iint }|f(x,y)|d\sigma (三角不等式的推广)\\ & (3)若m\le f(x,y)\le M,则m{A}_D\le \underset{D}{\iint }f(x,y)d\sigma \le M{A}_D\end{array}$$

• 性质5 (中值定理) ：若 $D$ 是有界闭区域，$f(x,y)\in C(D)$，则存在 $(\xi ,\eta )\in D$，
  $$\underset{D}{\iint }f(x,y)d\sigma =f(\xi ,\eta )A_D$$
  对应一元函数定积分中的平均值定理（积分第一中值定理）

第二节 二重积分的计算

一、直角坐标系下的计算

01 $x$ 型正则区域

设区域 D = { ( x , y ) ∣ a ≤ x ≤ b , φ 1 ( x ) ≤ y ≤ φ 2 ( x ) } D=\left\{(x, y) \mid a \leq x \leq b, \varphi_{1}(x) \leq y \leq \varphi_{2}(x)\right\} D={(x,y)∣a≤x≤b,φ1​(x)≤y≤φ2​(x)}

二重积分 $\underset{D}{\iint }f(x,y)dxdy$ $(d\sigma =dxdy)$ 的值

等于以 $D$ 为底，以曲面 $S$ : $z=f(x,y)$ 为顶的曲顶柱体的体积。

利用定积分来求体积考虑垂直 $x$ 轴过 $x$ 处的平面截曲顶柱体所得截面积 $A(x)$

截面曲边梯形的面积 $A(x)$
$$\mathrm{A}(x)={\int }_{{\phi }_1(x)}^{{\phi }_2(x)}f(x,y)dy$$
可得曲顶柱体的体积
$$V={\int }_a^b\mathrm{A}(x)dx={\int }_a^b\left({\int }_{{\phi }_1(x)}^{{\phi }_2(x)}f(x,y)dy\right)dx$$
导出
$$\underset{D}{\iint }f(x,y)dxdy={\int }_a^b\left({\int }_{{\phi }_1(x)}^{{\phi }_2(x)}f(x,y)dy\right)dx$$
写成
$$\underset{D}{\iint }f(x,y)dxdy={\int }_a^{b}dx{\int }_{{\phi }_1(x)}^{{\phi }_2(x)}f(x,y)dy$$

02 $y$ 型正则区域

若积分区域
D = { ( x , y ) ∣ ψ 1 ( y ) ≤ x ≤ ψ 2 ( y ) , c ≤ y ≤ d } D=\left\{(x, y) \mid \psi_{1}(y) \leq x \leq \psi_{2}(y), c \leq y \leq d\right\} D={(x,y)∣ψ1​(y)≤x≤ψ2​(y),c≤y≤d}
则有
$$\underset{D}{\iint }f(x,y)dxdy={\int }_c^{d}dy{\int }_{{\psi }_1(y)}^{{\psi }_2(y)}f(x,y)dx$$
小总结

对于一般区域的二重积分可将其分成若干个正则子区域，
利用积分的可加性，分别在各子区域积分后求和。

当积分区域关于 $x$ 轴或 $y$ 轴对称时，注意被积函数是否有奇偶性，从而使积分简化。

对称性非常重要！

二、极坐标系下的计算公式

当积分区域的边界曲线或被积函数用极坐标表示较为简单时，二重积分有时可用极坐标来计算

我们来考虑面积元素 $\Delta \sigma$ 在极坐标下的形式

用 $r$ 为常数所表示的圆周族和 $\theta$ 为常数所表示的射线族分割区域 $D$，那么小区域面积
$$\begin{array}{rl}& \Delta \sigma =\frac{1}{2}\left[(r+\Delta r{)}^2\Delta \theta -r^2\Delta \theta \right]=\frac{1}{2}\left[2r\Delta r+(\Delta r{)}^2\right]\Delta \theta \\ & ⟹d\sigma =rdrd\theta \end{array}$$
从直角坐标变换为极坐标时的二重积分的变换公式
$$\underset{D}{\iint }f(x,y)dxdy=\underset{D}{\iint }f(r\cos \theta ,r\sin \theta )rdrd\theta$$
若区域 D = { ( r , θ ) ∣ r 1 ( θ ) ≤ r ≤ r 2 ( θ ) , α ≤ θ ≤ β } D =\left\{(r, \theta) \mid r_{1}(\theta) \leq r \leq r_{2}(\theta), \alpha \leq \theta \leq \beta\right\} D={(r,θ)∣r1​(θ)≤r≤r2​(θ),α≤θ≤β}

二重积分化为累次积分
$$\underset{D}{\iint }f(r\cos \theta ,r\sin \theta )rdrd\theta ={\int }_{\alpha }^{\beta }d\theta {\int }_{r_1(\theta )}^{r_2(\theta )}f(r\cos \theta ,r\sin \theta )rd\theta$$

三、二重积分的变量代换

设变换 $\{\begin{array}{l}x=x(u,v)\\ y=y(u,v)\end{array}$ 有连续偏导数，且满足
$$J=\frac{\partial (x,y)}{\partial (u,v)}=\left|\begin{array}{ll}{x}_u& y_u\\ x_v& y_v\end{array}\right|\ne 0$$
而 $f(x,y)\in C(D)$，那么
$$\underset{D}{\iint }f(x,y)dxdy=\underset{D^{\prime }}{\iint }f(x(u,v),y(u,v))|J|dudv$$
$uv$ 平面小矩形 $A^{\prime }{B}^{\prime }{C}^{\prime }{D}^{\prime }⟶xy$ 平面曲边四边形 $ABCD$
$$\begin{array}{rl}& A^{\prime }(u,v)⟶A(x(u,v),y(u,v))\\ & B^{\prime }(u+\Delta u,v)⟶B(x(u+\Delta u,v),y(u+\Delta u,v))\\ & C^{\prime }(u+\Delta u,v+\Delta v)⟶C(x(u+\Delta u,v+\Delta v),y(u+\Delta u,v+\Delta v))\\ & D^{\prime }(u,v+\Delta v)⟶D(x(u,v+\Delta v),y(u,v+\Delta v))\end{array}$$
（此处应该有图）

$ABCD$ 近似平行四边形，只需求出一组邻边的向量表示
$$\begin{array}{rl}& \overrightarrow{AB}\approx \left(x_u(u,v)\Delta u,y_u(u,v)\Delta u\right)\\ & \overrightarrow{AD}\approx \left(x_v(u,v)\Delta v,y_v(u,v)\Delta v\right)\\ & \Rightarrow \\ & \Delta \sigma \approx \left|\begin{array}{ll}{x}_u(u,v)\Delta u& x_v(u,v)\Delta u\\ y_u(u,v)\Delta v& y_v(u,v)\Delta v\end{array}\right|=|J|\Delta u\Delta v\\ & \\ & d\sigma =|J|dudv\end{array}$$

第三节 三重积分

一、三重积分的定义

00 引入

${\iiint }_{\Omega }f(x,y,z)dV$ 你能说出它的含义吗？

概念联系的问题

一个占据三维空间中区域 $Q$ 的几何体，其密度为 $f(x,y,z)$，那么其质量为多少?

回顾定积分和二重积分的概念

求在三维区域上分布率非均匀的某种物理量 (或其它量) 的总量

分割—求和—求极限

01 定义

设 $\Omega$ 是 $R^3$ 中的有界闭区域，函数 $f(x,y,z)$ 在 $\Omega$ 上定义，$I$ 为实数，

若将区域 $\Delta {\Omega }_1,\Delta {\Omega }_2,\cdots ,\Delta {\Omega }_n$，任取 $\left({\xi }_i,{\eta }_i,{\varsigma }_i\right)\in \Delta {\Omega }_i$，作和
$$\sum _{i=1}^{n}f\left({\xi }_i,{\eta }_i,{\varsigma }_i\right)\Delta V_i(\Delta V_i是\Delta {\Omega }_i的体积)$$
总有
$$\underset{i\to 0}{\lim }\sum _{i=1}^{n}f\left({\xi }_i,{\eta }_i,{\varsigma }_i\right)\Delta V_i=I$$
( 其中 λ = max ⁡ 1 ≤ i ≤ n { d i } , d i \lambda=\max _{1 \leq i \leq n}\left\{d_{i}\right\}, d_{i} λ=max1≤i≤n​{di​},di​ 是小区域 $\Delta {\Omega }_i$ 的直径 )

则称函数 $f(x,y,z)$ 在 $\Omega$ 可积，$I$ 称为 $f$ 在 $\Omega$ 的三重积分，记为
$${\iiint }_{\Omega }f(x,y,z)dV(dV-体积元素)$$

02 意义

一种物理意义（三维物体的质量）

若 $f(x,y,z)$ 表示占有三维空间区域 $Q$ 的物体的质量密度函数，则

$\underset{\Omega }{\iiint }f(x,y,z)dV$ 给出了物体的质量

02 性质

(1) 类似二重积分，有线性、可加性、单调性和中值定理，还有

(2) $\underset{\Omega }{\iiint }1dV=V_{\Omega }(\Omega 的体积)$

二、在直角坐标系下的计算公式

直角坐标系下
$$\underset{\Omega }{\iiint }f(x,y,z)dV=\underset{\Omega }{\iiint }f(x,y,z)dxdydz$$

01 柱线法

设 $\Omega$ 是以曲面 $\mathrm{z}={\mathrm{z}}_1(x,y)$ 为底，曲面 $\mathrm{z}={\mathrm{z}}_2(x,y)$ 为顶，

而侧面是母线平行 $z$ 轴的柱面所围成的区域。

设 $\Omega$ 在 $xy$ 平面上的投影区域为 $D$ ，则 $\Omega$ 可表示为（$xy$ 型正则区域）
{ ( x , y , z ) ∣ z 1 ( x , y ) ≤ z ≤ z 2 ( x , y ) , ( x , y ) ∈ D } \left\{(x, y, z) \mid z_{1}(x, y) \leq z \leq z_{2}(x, y),(x, y) \in D\right\} {(x,y,z)∣z1​(x,y)≤z≤z2​(x,y),(x,y)∈D}
从质量角度求三重积分，则 $f(x,y,z)$ 为密度，对 $(x,y)\in D$，
$$\mu (x,y)={\int }_{z_1(x,y)}^{z_2(x,y)}f(x,y,z)dz$$
给出了 $\Omega$ 内由 $z_1(x,y)$ 到 $z_2(x,y)$ 的线段上所分布的质量密度。

物体的总质量就是
$$\underset{D}{\iint }\left({\int }_{z_1(x,y)}^{z_2(x,y)}f(x,y.z)dz\right)dxdy$$
从而
$$\underset{\Omega }{\iiint }f(x,y,z)dxdydz=\underset{D}{\iint }dxdy{\int }_{z_1(x,y)}^{z_2(x,y)}f(x,y\cdot z)dz$$

02 截面法

设区域 $\Omega$ 在 $z$ 轴上投影区间为 $[h_1,h_2]$，即 $\Omega$ 介于平面 $\mathrm{z}=h_1$ 与 $\mathrm{z}=h_2$ 之间，

垂直 $z$ 轴过 $z$ 处的平面截 $\Omega$ 所得截面为区域 $D_z$，则（$z$ 型空间区域）
Ω = { ( x , y , z ) ∣ ( x , y ) ∈ D z , h 1 ≤ z ≤ h 2 } \Omega=\left\{(x, y, z) \mid(x, y) \in D_{z}, h_{1} \leq z \leq h_{2}\right\} Ω={(x,y,z)∣(x,y)∈Dz​,h1​≤z≤h2​}
仍从质量角度考虑，对 $z\in [h_1,h_2]$，二重积分
$$F(z)=\underset{D_z}{\iint }f(x,y,z)dxdy$$
给出了物体在截面 $D_z$ 上所分布的质量。物体的总质量为：
$${\int }_{z_1}^{z_2}\left(\underset{D_z}{\iint }f(x,y,z)dxdy\right)dz$$
从而有：
$$\underset{\Omega }{\iiint }f(x,y,z)dV={\int }_{h_1}^{h_2}dz\underset{D_x}{\iint }f(x,y,z)dxdy$$

三、三重积分变量代换

与二重积分的变量代换类似，

设变换 $\{\begin{array}{l}x=x(u,v,w)\\ y=y(u,v,w)\\ z=z(u,v,w)\end{array}$ 有连续偏导数，且满足
$$J=\frac{\partial (x,y,z)}{\partial (u,v,w)}=\left|\begin{array}{lll}{x}_u& y_u& z_u\\ x_v& y_v& z_v\\ x_w& y_w& z_w\end{array}\right|\ne 0$$
而 $f(x,y,x)\in C(\Omega )$，那么
$$\underset{\Omega }{\iiint }f(x,y,z)dxdydz=\underset{\Omega }{\iiint }f(x(u,v,w),y(u,v,w),z(u,v,w))|J|dudvdw$$

01 柱面坐标系下的三重积分

这个坐标系实际上就是 $xy$ 坐标转变为极坐标，即变换公式为 $\{\begin{array}{c}x=r\cos \theta \\ y=r\sin \theta \\ z=z\end{array}$

由于
$$\frac{\partial (x,y,z)}{\partial (r,\theta ,z)}=\left|\begin{array}{ccc}\cos \theta & \sin \theta & 0\\ -r\sin \theta & r\cos \theta & 0\\ 0& 0& 1\end{array}\right|=r$$
得到柱坐标积分公式
$$\underset{\Omega }{\iiint }f(x,y,z)dV=\underset{\Omega }{\iiint }f(r\cos \theta ,r\sin \theta ,z)rdrd\theta dz$$
注意，事实上，在具体计算时，可以用柱线法或截面法得到 $D$ ( 或 $D_z$ ) 的二重积分，再转化为极坐标。

02 球面坐标系下的三重积分

设点 $M(x,y,z)$ 是空间一点，引进球坐标 $(\rho ,\phi ,\theta )$

$\rho =\Vert \overrightarrow{OM}\Vert$，$\phi :\overrightarrow{OM}$ 与z轴正向的夹角，$\theta :\overrightarrow{OP}$ 与 $x$ 轴正向的夹角，

且 $\rho$，$\phi$，$\theta$ 满足
$$0\le \rho \le +\infty ,0\le \phi \le \pi ,0\le \theta \le 2\pi 或-\pi \le \theta \le \pi$$
坐标变换关系式$⟹\{\begin{array}{c}x=\rho \sin \phi \cos \theta \\ y=\rho \sin \phi \sin \theta \\ z=\rho \cos \phi \end{array}$

由于雅可比行列式
$$\begin{array}{rl}\frac{\partial (x,y,z)}{\partial (\rho ,\phi ,\theta )}& =\left|\begin{array}{ccc}\sin \phi \cos \theta & \sin \phi \sin \theta & \cos \phi \\ \rho \cos \phi \cos \theta & \rho \cos \phi \sin \theta & -\rho \sin \phi \\ -\rho \sin \phi \sin \theta & \rho \sin \phi \cos \theta & 0\end{array}\right|={\rho }^2\sin \phi \end{array}$$
导出
$$\begin{array}{rl}& \underset{\Omega }{\iiint }f(x,y,z)dV=\underset{{\Omega }^{\ast }}{\iiint }f(\rho \sin \phi \cos \theta ,\rho \sin \phi \sin \theta ,\rho \cos \phi ){\rho }^2\sin \phi d\rho d\phi d\theta \end{array}$$
使用球坐标时，

$\rho =$ 常数:球面，$\phi =$ 常数: 锥面，$\theta =$ 常数: 平面，

且球面和锥面的中心在原点，平面过 $z$ 轴。

注意，围成区域的部分曲面有上述特点，或被积函数含 $x^2+y^2+z^2$，可考虑用球坐标

第四节 重积分的应用

一、曲面面积

设空间曲面 $S$ 为 $z=f(x,y)$，$(x,y)$ 定义于 $D$，即曲面 $S$ 在 $xy$ 平面的投影区域为 $D$，如何求 $S$ 的面积？

用分割求和（微元法）的思想将 $D$ 分割成小区域，对应 $D$ 上小区域面积 $\Delta \sigma =\Delta x\Delta y$ ，$S$ 上的小曲面面积为 $\Delta S$，当小区域微小时，$\Delta S|\cos \phi |\approx \Delta \sigma$，其中 $\phi$ 是 $\Delta S$ 所在平面与 $xy$ 平面夹角。

$\phi$ 是 $\Delta S$ 上的法向量与 $z$ 方向的夹角，这两个方向向量分别为 { z x , z y , − 1 } , { 0 , 0 , 1 } \left\{z_{x}, z_{y},-1\right\},\{0,0,1\} {zx​,zy​,−1},{0,0,1}

导出
$$|\cos \phi |=\frac{1}{\sqrt{1+z_x^2+z_y^2}}$$
因此得到曲面面积的有关公式
$$\begin{array}{rl}& dS=\sqrt{1+z_x^2+z_y^2}dxdy(曲面面积元素)\\ & S=\underset{D}{\iint }\sqrt{1+z_x^2+z_y^2}dxdy\end{array}$$

二、重积分的物理应用举例

01 质心

物体的质心 ( 或重心 ) 与它的质量和静力矩有关。

设面密度为 $\mu (x,y)$ 薄板占据平面区域 $D$

考虑 $D$ 上面积元素 $dxdy$ 其质量为 $\mu (x,y)dxdy$，

从而对 $y$ 轴的静力矩为 $d{M}_y=x\mu (x,y)dxdy$

质量 $m$ 和静力矩 $M_y$ 为：$m=\underset{D}{\iint }\mu (x,y)d\sigma$，$M_y=\underset{D}{\iint }x\mu (x,y)d\sigma$

同样，对 $x$ 轴电静力矩 $M_x=\underset{D}{\iint }y\mu (x,y)d\sigma$

于是薄片的质心位置 $(\bar{x},\bar{y})$
$$\bar{x}=\frac{M_y}{m}=\frac{\underset{D}{\iint }x\mu (x,y)d\sigma }{\underset{D}{\iint }\mu (x,y)d\sigma },\bar{y}=\frac{M_x}{m}=\frac{\underset{D}{\iint }y\mu (x,y)d\sigma }{\underset{D}{\iint }\mu (x,y)d\sigma }$$
注意，当 $\mu =1$ 时得到平面图形的形心

思考：对三维物体如何求质心。

02 转动惯量

转动惯量也是一种矩 ( 二次矩 )，设平面区域 $D$ 上薄板的面密度为 $\mu (x,y)$ ，

那么面积元素 $d\sigma$ 处微量物体对 $y$ 轴的转动惯量
$$\begin{array}{rl}& d{I}_y=x^2\mu (x,y)d\sigma ⟹I_y=\underset{D}{\iint }{x}^2\rho (x,y)d\sigma \end{array}$$
问题：对 $x$ 轴和对原点 $O$ 的转动惯量?

最后

😊为防止河蟹，链接已经通过“与熊论道/熊曰加密”加密处理，将下面的文字复制到“与熊论道/熊曰加密”页面的第二个输入框，点击“领悟熊所言的真谛”即可查看链接啦：
😊熊曰：呋食食和嚁非象嗚家吃呱山萌萌笨有哞魚既魚性蜜覺呆食哮性洞哮山噗眠嗥嚄萌洞擊嗄襲呱物人你
😊如果嫌麻烦的话请私信咨询博主，谢谢！
😊PS：繁星依月/惟欢/一舟均为博主的马甲，本篇作品完全原创，再次感谢！