基于超声能量扩散的CFST 界面剥离损伤检测方法与验证

2023-01-04 07:19张海月高树灵李忠献
工程力学 2023年1期
关键词:扩散系数钢管尺寸

李 宁,张海月,高树灵,李忠献

(1. 天津大学建筑工程学院,天津 300350;2. 滨海土木工程结构与安全教育部重点实验室(天津大学),天津 300350;3. 中国地震局地震工程综合模拟与城乡抗震韧性重点实验室(天津大学),天津 300350)

近年来,由于钢管混凝土(CFST)结构在强震作用下延性好、刚度大、承载能力高等优点得益于其施工过程快捷高效,正被广泛应用于各类型工程结构部件中[1]。CFST 结构在施工过程中,核心混凝土的浇筑方式、施工工艺及温度变化导致的混凝土徐变收缩、CFST 构件受压初期产生的负围压等不利因素,都可能引起核心混凝土与钢管界面的剥离损伤[2]。一方面,界面剥离损伤削弱了钢管对核心混凝土的约束作用,导致CFST 极限承载能力、延性及刚度降低[3],因此钢管与核心混凝土间的粘结情况一直是一个备受关注的问题;另一方面,由于CFST 构件界面剥离损伤的隐蔽性,传统意义上的无损检测手段无法直接检测到,因而开展界面粘结状态的检测具有重要的工程应用需求和实际意义[4]。

针对上述问题和需求,许多学者基于超声检测对CFST 剥离损伤检测开展了试验及仿真研究。LIANG 等[5]基于超声阻抗法研究了钢板与混凝土间的粘结滑移情况,结合均方根差识别滑移发展过程。ZHANG 等[6]利用压电智能骨料监测低频循环载荷作用下L 型CFST 柱内部损伤状况,采用小波包分析损伤指数进行数据分析。许斌等[7]提出基于PZT 的CFST 柱主动界面剥离损伤检测方法,并采用小波包能量谱的均方根偏差作为评价指标。栾乐乐等[8]建立CFST 界面剥离谱元法模型,探明了应力波在其内部的传播机理。CHEN等[9]利用随机聚合方法模拟核心混凝土,数值分析应力波通过CFST 构件不同剥离损伤及微观骨料时的传播过程。

以上研究表明,对CFST 界面剥离损伤定量识别迄今为止研究较少,且大部分注重定性地判别CFST 损伤程度,判别损伤范围十分有限。为了更加全面提取信号损伤信息、定量的对CFST 剥离损伤程度进行描述和检测,一种对微小损伤敏感、识别准确且能直接表征损伤的方法—超声能量扩散方法逐渐被广泛应用[10]。

已有学者基于超声能量扩散方法,对不同构件中存在的损伤进行定量识别。QUIVIGER 等[11]利用超声能量扩散有效识别了混凝土中的裂纹,并判别了裂纹数量及深度。IN 等[12]通过扩散超声波原位监测开裂混凝土在模拟海洋环境下的自愈合过程,表明扩散率可以有效预测自愈趋势。FARID 等[13]基于超声能能量扩散法表征钢筋混凝土板在弯曲试验中的机械损伤,结论表明该方法能够区分混凝土构件的损伤阶段。AHN 等[14]利用能量扩散系数D及耗散系数σ 表征混凝土试样的微裂纹,并探讨了波频、压电换能器的间距对超声扩散的影响。LU 等[15]研究了结构载荷对超声能量扩散的影响,得出能量扩散系数与拉伸力指数的关系曲线,在此基础上可进行螺栓松动预测。WANG 等[16]提取超声扩散信号损伤信息来识别列车轨道的缺陷,并对钢轨健康状态进行连续监测,证明该方法在实际铁路应用中的有效性。以上研究表明,超声能量扩散方法中提取的能量扩散系数D及耗散系数 σ 可直接表征损伤,根据两参数数值可确定损伤尺寸,可提取应力波损伤信息,对损伤分析更全面,因此本文利用超声能量扩散方法对CFST 剥离损伤进行定量检测。

本文研究表明,扩散及耗散系数随界面剥离损伤变化规律明显,在此基础上所提出的判定模型可对钢管与核心混凝土间的界面剥离损伤进行准确、合理的定量评价,为工程实际中的超声检测应用提供了重要的理论与方法支撑。

1 超声能量扩散分析理论

应力波在CFST 中的传播极其复杂,能量衰减主要包括2 方面:① 由于核心混凝土[17]的内部组成物几何尺寸的差异,从极微小的未水化水泥颗粒到5 mm~ 40 mm 的粗骨料,颗粒的固有频率各异,当某种频率成分通过该介质时,相应自振频率的颗粒便产生共振,形成新的球面波源使应力波能量不断减弱;② 钢管与混凝土组合界面以及粗、细骨料的边界的反射、折射,使应力波发生不规则的散射、衰减和模态转换现象[14]。可见,提取应力波在界面剥离处传播携带的损伤信息是分析剥离损伤的有效方法。

应力波在CFST 试件(有剥离损伤)中的传播过程如图1 所示,PZT 驱动器在输入电压信号的激励下产生应力波,主要由钢管外壁处产生的Lamb 波[18]和在核心混凝土内部产生的体波组成,经过带有剥离损伤的钢管混凝土后由试件另一侧的压电接收器接收,本文着重于对接收波进行整体能量分析。Lamb 波主要在钢管薄壁上传播,体波在钢管及混凝土中传播,其中Lamb 波和体波皆存在剥离损伤反射波信息,但是由于Lamb 波为面波能量消散较快且只在钢板中产生[18],因此Lamb 中损伤成分较少,在混凝土中传播的体波所携带的损伤信息更多。当存在大面积剥离损伤时,应力波在损伤边界的反射、折射导致能量衰减及波形的变化,因此本文对接收应力波整体分析其能量特征是识别损伤的有效途径。WEAVER等[19]基于超声能量扩散理论,对应力波在混凝土等复合材料中的能量衰减机制进行了分析,结果表明应力波在传播介质中的能量衰减主要包括2 个方面:① 在传输介质中的细小颗粒处产生的能量扩散,其能量损失主要与微颗粒的几何参数及结构参数有关;② 由于介质振动而产生的能量耗散,主要与基体及介质的材料属性相关。随后,WEAVER[21]分析了应力波在泡沫铝板中的扩散形式,应用能量扩散及耗散系数进行表征,表明上述参数有一定的研究意义。

图1 超声波在钢管混凝土中的扩散及传播Fig. 1 Diffusion and propagation of ultrasonic wave in CFST

假设应力波的扩散与热传导过程存在相似的机理,有学者[19−22]利用修正的热传导控制方程来描述应力波能量的变化规律,观测点P 在t时刻的能量变化可用三维超声能量扩散方程表示:

式中:D和σ 分别代表超声能量扩散系数和耗散系数;P(x,y,z,t,f)代表超声波激励信号能量谱密度;E(x,y,z,t,f)代表超声波能量谱密度函数;x,y,z是三维笛卡尔坐标系的坐标。

式(1)中,D主要受介质内部结构的影响较大,对内部微粒及边界反射较敏感,包括随机分布的晶粒及介质中产生的损伤,其值大小与最大能量密切相关。σ 主要与传播介质的粘弹特性及介质本身的能量吸收相关。损伤的存在改变了应力波的传播路径、波能及能量耗散速度,理论上,参数D和σ 皆可有效辨识CFST 剥离损伤尺寸。

式(1)是扩散控制方程,本文对其简化考虑,三维结构中的超声能量谱密度可表示为[10]:

式中:E0为初始超声能量密度;r为信号激励源与接收器之间的距离。为了求解超声能量扩散拟合曲线,式(2)可用对数形式表示:

其中,C0代表与信号初始能量有关的常数。根据此方程,通过数值拟合的方法即可求解参数D和σ。

2 钢管混凝土剥离缺陷检测试验

基于《钢管混凝土结构技术规范》(GB 50936−2014)[23]中对实心钢管凝土柱的设计规定及试验设备的限制,本节选取试件尺寸长×宽×高=100 mm×100 mm×300 mm,外钢管采用Q235 钢材,核心混凝土采用C30。共设计4 组CFST 试件,每组2 个,各试件4 面皆由不同的损伤工况组成,均布置无损伤面作对比。使用带凹槽的亚克力板模拟剥离损伤,在浇筑混凝土之前贴于钢管内壁表面,用有机硅密封胶粘贴好亚克力板四周,防止水分进入。剥离损伤厚度为4 mm。试件各面均贴有对称的压电传感器用来发射与接收信号,具体剥离损伤布置见图2,损伤工况见表1,剥离损伤尺寸面积占比为0%~28%,涵盖了工程中出现剥离损伤尺寸面积的大多数情况。其中损伤面积占比Q=SD/STube(剥离损伤面积/相应钢管面面积)。

图2 钢管混凝土剥离损伤布置图Fig. 2 Debonding defects layout in CFSTs

表1 剥离损伤工况尺寸表Table 1 Debonding defects condition size

图3 为试验装置示意图。图4 为试验试件及设备。函数发生器生成峰值为2V 的五峰波(不同频率)超声脉冲信号(发射频率为50 kHz、100 kHz、150 kHz、200 kHz、250 kHz);宽带功率放大器将脉冲信号放大20 倍,由粘贴在钢管外壁的压电驱动器传至CFST 内部,由贴在另一侧的压电接收器接收,以25 MHz/s 的采样率进行采集。

图3 试验装置示意图 /mm Fig. 3 Test device schematic diagram

图4 试验构件及设备Fig. 4 Test specimens and equipment

3 数据处理

为了确定式(3)中C0、D和σ 的值,对信号进行频域分析。首先,将接收的应力波时域信号(图5)进行矩形窗划分,窗长∆t取0.14 ms,为了保持信号信息的完整性,窗重叠率取90%。其次,对分窗后的时域信号进行快速傅里叶变换求其频谱,计算各频窗的总能量,该能量为其对应的时域信号的扩散能量。最后,对计算得到的扩散能量取对数,按照式(3)对D和σ 进行拟合。

图5 试件在150 K 激发频率下的应力波时域图Fig. 5 Time domain diagram of stress wave of the specimen at 150 K excitation frequency

有学者[11]提出最大能量到达时间ATME对剥离损伤十分敏感,因此本文基于D、σ 及ATME三个参数对CFST 剥离损伤进行研究。数据处理流程如图6 所示。如图7 所示,对照时域中的超声波能量密度数据,以上D和σ 两个参数可通过式(3)回归分析来确定,其中的附加参数ATME可通过求解信号能量得到。拟合过程中使用了信赖域算法和加权平方和,最终得到的扩散系数D与耗散系数σ 用于表征本研究中试样的剥离损伤。

图6 超声能量扩散系数数据拟合流程图Fig. 6 Flow chart of ultrasonic energy diffusion coefficient data fitting operation

图7 不同激发频率下的超声能量密度曲线Fig. 7 Ultrasonic energy density curves at different excitation frequencies

4 结果与讨论

本文选取能量扩散系数D、耗散系数σ 及ATME定量研究上述参数与激发频率及剥离损伤尺寸的关系。由于混凝土是非均质多相凝聚体且弹性模量相对较小、信号能量衰减较快,因此对混凝土的超声检测常用频率范围一般在20 kHz ~100 kHz;而外部钢管是高密度材料且信号能量衰减较小,检测频率一般较大。对于两者组成的CFST 而言,由于CFST 中核心混凝土内部较为复杂[8],剥离损伤处于混凝土与钢管两者界面之间,使用的超声波频率较低[17],本文初选频率范围50 kHz ~ 250 kHz。

4.1 超声能量扩散系数D、耗散系数σ 及ATME与波频及剥离损伤尺寸的关系

超声能量扩散主要取决于材料特性以及测试配置(例如激发频率、环境及温度),为了消除本研究可能存在环境条件的影响,整个试验在无噪声,温度、湿度恒定的环境中进行。CFST 剥离损伤和激励频率对超声扩散形为的影响是研究重点。图7 是无剥离损伤的CFST(Q= 0)在5 种不同激励频率下测量的超声能量密度数据。通过回归分析得到相应频率超声密度数据的最佳拟合曲线。由式(3)分析可知,测量数据达到峰值(图8 中约为1 ms)的上升部分(即上升至峰值的时间及其幅度)由D决定。相比,图8 中2.5 ms~8 ms 的实测数据的下降部分(即斜率)由σ 决定[22]。

图8 扩散方程参数拟合示意图Fig. 8 Fitting plot of diffusion equation parameters

由图9、图10、图11 可知,随着波频的增加,应力波能量衰减加快,D减小,σ 增加,ATME减小。由此可知,激励频率对超声扩散影响显著,因此,激励频率是D及σ 的重要影响要素。

图9 不同激发频率超声能量扩散系数Fig. 9 Ultrasonic energy diffusion coefficient at different excitation frequencies

图10 不同激发频率超声能量耗散系数Fig. 10 Ultrasonic energy dissipation coefficient at different excitation frequencies

图11 不同波频下最大能量到达时间Fig. 11 Maximum energy arrival time at different wave frequencies

较小的扩散系数意味着从激励源到接收器的超声波能量不断衰减,因此,D值随剥离损伤尺寸的增加而减小,如图9 所示。界面剥离损伤的存在增加了应力波的反射、折射及入射波束向周围的散射,因此到达压电接收器的应力波能量不断衰减。可见D值是一个衡量CFST 剥离损伤尺寸的有效参数。

超声能量耗散系数σ 与传播介质的材料属性相关,由图10 可知,σ 随着剥离损伤尺寸的增加而增加,整个选定的频率范围内,健康工况的耗散值最低,由此可知,当CFST 存在剥离损伤时,信号能量耗散远大于健康工况。原因是剥离损伤的存在,传播介质的内部边界增加,应力波的散射增大,能量耗散增大,并且CFST 剥离损伤的传播介质主要为空气,应力波在气体中传播相比固体在固体中传播衰减更快,随着剥离损伤尺寸的增加,空气占比更大,因此能量衰减加快,σ 增大。

如图11 所示,一方面,在不同频率下ATME呈现下降的趋势,频率越高,ATME越小,原因是频率越高,应力波能量衰减越快;另一方面,随着损伤尺寸增大,ATME呈现下降的趋势。原因是,剥离损伤阻碍了应力波向其内部的传播,更多传播至钢管外壁,而波在钢管的传播速度远大于混凝土[3],因此ATME随着剥离损伤尺寸的增加而减小。由于ATME不呈现严格的递减趋势,有一定的数据波动,因此对于CFST 而言,此参数无法充分表征剥离损伤尺寸,故在以下研究中舍弃ATME,研究重点为D及σ。

4.2 最佳激励频率的选择

50 kHz~250 kHz 频率中,得出适合本试验的最佳激励频率并在此基础上提出损伤判定模型。将不同波频及相应的D及σ 进行归一化,如图12、图13 所示。首先归一化扩散及耗散系数曲线,求解拟合直线斜率,斜率越大,意味着D和σ 随剥离损伤尺寸波动越大,对损伤尺寸判别灵敏度较高,D和σ 线性拟合后的斜率值如表2 所示。

表2 归一化D 和σ 线性拟合斜率绝对值Table 2 Unified absolute value of slope of normalized Dand σ with linear fitting

图12 超声能量扩散系数归一化D 值Fig. 12 Normalized Dvalue of ultrasonic energy diffusion coefficient

图13 超声能量耗散系数归一化σ 值Fig. 13 Normalized σ value of ultrasonic energy dissipation coefficient

由表2 可知,在激励频率150 kHz 下,基于D值的斜率最大且线性拟合系数R2= 97%。对σ,150 kHz 和200 kHz 斜率皆较大,150 kHz 的线性拟合系数R2= 97.28%>95.84%(200 kHz)。因此,150 kHz 为最佳激励频率,即此频率下D和σ 对损伤更为敏感。

5 基于超声能量扩散定量评估CFST剥离损伤尺寸

基于上述数据分析结论,ATME 对损伤灵敏度较低,其变化不呈现严格的递减趋势,D和σ 均随剥离损伤尺寸增加,呈现线性变化的趋势,因此基于D和σ 提出CFST 剥离损伤判定模型。通过线性拟合最佳激励频率150 kHz 下D和σ 与剥离损伤尺寸的判定曲线,提出剥离损伤判定模型。基于判定模型可根据归一化D和σ 值计算出剥离损伤面积占比Q,判定损伤面积实现定量检测。判定模型如表3 所示。

表3 D 及σ 剥离损伤判定模型Table 3 Debonding defects criteria of Dand σ

如图14、图15 为线性拟合150 kHz 波频下D和σ 与剥离损伤尺寸的判定曲线,其线性拟合系数分别为97.40%、97.28%,对于在剥离损伤面积占比为0%~28%之间的D值,拟合判定模型与判定曲线最大差值为损伤面积占比6%对应下的8.40%,最小差值为损伤面积占比18%对应下的1.9%,整体平均误差为3.60%。对于σ 值,除了无损伤工况及损伤面积占比28%工况下,其他工况下两者值非常相近,在此期间最大差值为3.76%,最小差值为0.03%,整体平均误差2.53%,均小于5%且线性拟合系数均大于97%,因此,该损伤判定模型是准确的。

图14 归一化D 值线性拟合Fig. 14 Normalized Dvalue with linear fitting

图15 归一化σ 值线性拟合Fig. 15 Normalized σ value linear fitting

6 基于COMSOL 的CFST 剥离检测模型与验证

本节研究是对上述损伤判定模型的仿真验证。考虑试验构件大小、剥离尺寸条件,模拟了不同界面剥离损伤情况下CFST 构件的应力波传播过程,基于仿真对试验提出的判定模型进行了仿真验证。

6.1 基于时域信号的COMSOL 模型验证

超声能量扩散法主旨是分析时域信号信息,因此对比试验及仿真时域信号验证模型的准确性。图16 为试验及仿真时域信号对比。可见仿真信号幅值略大于试验信号,原因是试验过程中由于环境、温度的影响,存在应力波能量损耗,因此幅值略低。两者幅值均有下降趋势且波形相似,表明仿真信号较准确。

图16 试验及仿真信号对比Fig. 16 Comparison of experimental and simulated signals

6.2 基于波速的COMSOL 模型验证

根据同等条件下的首波波速大小进一步验证模型的准确性。图17 为试验条件下,首波到达的时间Δt1= 5×10−5s,以及两传感器之间的距离s=0.21 m,试验波速V1= 0.21/5 × 10−5= 4200 m/s。图18 为仿真条件下,首波到达时间Δt2= 5.16 ×10−5s,两传感器之间的距离仍为s,模拟波速V2=0.21/5.16 × 10−5= 4069.77 m/s,两者的误差百分比为3.1% < 5%,因此基于波速分析的有限元模拟结果与试验相比是准确的。

图17 试验条件下激励与接收信号时间差Fig. 17 Time difference between excitation and reception signal during tests

图18 仿真条件下激励与接收信号时间差Fig. 18 Time difference between excitation and reception signal in simulation

6.3 基于应力波传播理论对COMSOL 模型准确性的验证

本节模拟了CFST 存在界面剥离损伤时的应力波的传播过程,通过应力波传播理论对照仿真现象进行模型验证。为了简易计算,核心混凝土视为均匀材料。应力波穿过钢管和核心混凝土之间的界面时,存在反射、衍射及散射现象,且应力波形会发生改变。为了凸显应力波对界面剥离损伤的敏感性及传播模式,应用COMSOL 有限元软件仿真验证,更直观观察应力波传播过程,以损伤工况为A12(剥离损伤尺寸为60 mm×140 mm×4 mm)的CFST 为例。由于探讨界面剥离损伤对应力波传播过程的影响,因此以钢管面为主方向观察应力波传播过程。

如图19 所示,首先压电驱动器激发应力波,主要由Lamb 波和体波组成[18](图19(a))。Lamb 波主要在钢管薄壁上传播;体波向内部核心混凝土传播,经界面剥离损伤会在钢管表面发生反射弹射。在3.80 × 10−5s 时刻,应力波不断发散,经钢管边界发生反射并沿钢管向前传播(图19(b)、图19(c)),此时出现了清晰的波阵面,随后边界反射波逐渐耗散,应力波继续向前传播,在剥离损伤中心处发生反射现象(图19(d)、图19(e)),最后在损伤区域不断反射,携带损伤信息(图19(f))。

图19 应力波在带有剥离损伤的CFST 外壁表面传播过程Fig. 19 Propagation of stress waves on the surface of CFST with debonding defects

理论上,应力波在传播过程中不断反射、折射,能量逐渐耗散且波能伴随着衰减。界面剥离损伤的存在增加了钢板边界条件的复杂性,导致了波不同方向的散射使其波形发生变化。且界面剥离减弱了从钢管到核心混凝土的传播,导致应力波在剥离区域反射增加,同样增加应力波的绕射,这也是存在剥离损伤时传播时间对比于健康工况下更长的原因。研究表明,在此模型中应力波遇到剥离损伤的传播过程与理论上保持一致,证明了仿真的准确性。

6.4 基于模拟分析的剥离损伤判定模型验证

利用COMSOL 有限元软件对CFST 3 种剥离损伤工况H、A5、A12(Q= 0%、14%和28%)进行数值模拟,目的是得到有限元分析时域信号,求解归一化D及σ 值对比于判定模型计算求解的归一化D值与σ 值,验证判定模型的准确性。

由图20 可知,归一化仿真与判定模型D值在不同剥离损伤面积占比时、表现出不一致的误差,当剥离损伤面积占比分别为0%、14%和28%时,两者求解的D值误差分别为3.03%、3.08%和9.33%。其中误差最大时对应的是剥离损伤占比28%时的工况,其归一化模拟D值为25.44%,判定模型D值为16.11%,结果表明剥离损伤面积占比越大增加了应力波传播的不确定性因素,也使得模型仿真与试验效果差别越大。由于所有误差项均小于10%,表明本文拟合的扩散系数判定模型是准确的。在此基础上进一步对比了在不同损伤面积占比下的归一化仿真σ 值与判定模型σ 值之间的误差,如图21 所示。随剥离损伤面积占比的增加,归一化仿真σ 值与判定模型求解的σ 值误差表现出和D值类似的特征。当不存在剥离损伤面积时误差为9.06%,且随损伤面积占比的增大误差项增加,但是仿真及试验的σ 值呈现相似的变化趋势,均线性增加,这表明σ 值判定模型同样准确,能合理地表征CFST 构件剥离损伤的耗散特性。

图20 归一化仿真D 值与判定模型D 值对比Fig. 20 Comparison of normalized simulation Dvalues and the debonding criteria Dvalues

图21 归一化仿真σ 值与判定模型σ 值对比Fig. 21 Comparison of normalized simulation σ values and the debonding criteria σ values

7 结论

本文利用超声能量扩散法拟合超声扩散及耗散系数,定量表征CFST 剥离损伤特性,研究了在不同波频下、不同剥离损伤尺寸(剥离面积占比0%~28%)对扩散、耗散系数及ATME 的影响,并且基于归一化D及σ 与损伤尺寸的关系曲线,提出CFST 剥离损伤判定模型,经有限元仿真最后通过误差分析验证了判定模型的准确性。主要结论如下:

(1) 超声波在CFST 中传播的波频选为50 kHz~250 kHz。随着波频的增加,扩散系数线性减小,耗散系数线性增大,ATME呈现下降的趋势。表明频率是影响能量扩散的一个重要参数,为其他结构超声检测频率选取提供参考。

(2) 对剥离损伤尺寸面积占比为0%~28%,涵盖了工程中出现剥离损伤尺寸的大多数情况。在最佳激励频率150 kHz 下,基于归一化D及σ 值分别提出剥离损伤判定模型,该模型可对应D及σ 值求解得出CFST 剥离损伤面积占比,实现定量检测。

(3) 为了验证判定模型准确性,在试验相同条件下,利用有限元软件COMSOL 仿真,对剥离损伤面积不同的模型进行了有限元分析,从时域信号、波速及应力波传播理论3 方面验证了仿真模型的准确性。在此基础上,再通过归一化仿真D及σ 值与通过判定模型求解得到的D及σ 值对比,验证了所提出损伤判定模型的准确性。

本文初步采用超声能量扩散方法对CFST 界面剥离损伤定量检测,根据扩散及耗散系数提出的剥离损伤判定模型。随着传感器工艺、性能提升,后续在实际工程中的应用尚待进一步探索。

猜你喜欢
扩散系数钢管尺寸
表观扩散系数值与肝细胞癌分级的相关性以及相关性与肿瘤大小关系的分析
微型钢管桩在基坑支护工程中的应用
时变因素影响下混凝土中氯离子扩散计算方法
CIIE Shows Positive Energy of Chinese Economy
浅探输变电钢管结构的连接方法
ACS6000中压传动系统在钢管轧制中的应用
D90:全尺寸硬派SUV
定位于材料基因组计划的镍基高温合金互扩散系数矩阵的高通量测定
非肿块型乳腺癌的MR表观扩散系数及肿瘤大小与Ki-67表达的相关性研究
佳石选赏