CN112395676A - 含阻尼柔性结构抗线性爆炸荷载设计动力系数方法 - Google Patents

含阻尼柔性结构抗线性爆炸荷载设计动力系数方法 Download PDF

Info

Publication number
CN112395676A
CN112395676A CN202011299957.2A CN202011299957A CN112395676A CN 112395676 A CN112395676 A CN 112395676A CN 202011299957 A CN202011299957 A CN 202011299957A CN 112395676 A CN112395676 A CN 112395676A
Authority
CN
China
Prior art keywords
vibration
elastic
load
displacement
plastic
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Withdrawn
Application number
CN202011299957.2A
Other languages
English (en)
Inventor
张冰冰
耿少波
范志强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Taiyuan Institute of Technology
Original Assignee
Taiyuan Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Taiyuan Institute of Technology filed Critical Taiyuan Institute of Technology
Priority to CN202011299957.2A priority Critical patent/CN112395676A/zh
Publication of CN112395676A publication Critical patent/CN112395676A/zh
Withdrawn legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Architecture (AREA)
  • Civil Engineering (AREA)
  • Operations Research (AREA)
  • Structural Engineering (AREA)
  • Evolutionary Computation (AREA)
  • Buildings Adapted To Withstand Abnormal External Influences (AREA)

Abstract

本发明涉及一种含阻尼柔性结构抗线性爆炸荷载设计动力系数方法，属于抗爆设计技术领域，具体的柔性结构指的是：在线性爆炸荷载作用时长ti范围内，结构振动未达到弹性振动最大值，爆炸荷载卸载后，依靠惯性力在tT时刻达到弹性位移最大值yT后，进一步发展为塑性振动，在某一时刻tm，达到了结构总的弹塑性位移最大值ym，根据爆炸对建筑结构的作用过程，将该过程分为弹性阶段强迫振动、弹性阶段自由振动和塑性阶段自由振动；通过对上述不同的阶段的计算，充分考虑结构类型及结构阻尼对爆炸荷载动力系数的影响，尽量使设计的建筑结构更符合实际，在降低成本的同时，满足抗爆要求。

Description

含阻尼柔性结构抗线性爆炸荷载设计动力系数方法
技术领域
本发明涉及一种含阻尼柔性结构抗线性爆炸荷载设计动力系数方法，属于抗爆设计技术领域。
背景技术
目前，现有的建筑物在进行抗爆设计时，均需要考虑抗爆结构构件的抗爆能力。常规爆炸施加在结构上的爆炸荷载作用时长ti很短，我国及国外人防结构进行抗爆设计时，按等冲量线性荷载处理，进一步将该爆炸动荷载超压峰值Δpm与结构弹塑性阶段抗力动力系数kh相乘后，作为静载进行结构抗爆设计数值。其中规范给出的动力系数kh公式，未考虑结构阻尼的影响，也未考虑进行抗爆设计的结构类型，导致对实际结构精准化设计缺乏有力支撑，也造成了建筑结构在抗爆设计时，设计偏保守，在部分应用环境中，这种设计，将造成建设成本的增加。
发明内容
为解决现有技术存在的技术问题，本发明提供了一种含阻尼柔性结构抗线性爆炸荷载设计动力系数方法。
为实现上述目的，本发明所采用的技术方案为含阻尼柔性结构抗线性爆炸荷载设计动力系数方法，所述的柔性结构指的是：在爆炸荷载作用时长ti范围内，结构振动未达到弹性振动最大值，爆炸荷载卸载后，依靠惯性力在tT时刻达到弹性位移最大值yT后，进一步发展为塑性振动，在某一时刻tm，达到了结构总的弹塑性位移最大值ym，
根据爆炸对建筑结构的作用过程，将该过程分为弹性阶段强迫振动、弹性阶段自由振动和塑性阶段自由振动三个阶段；
a、弹性阶段强迫振动
在弹性阶段且在荷载作用时长范围0＜t≤ti内，动力等效体系的运动微分方程为：
Figure BDA0002786490650000021
其中，t为柔性结构爆炸作用下的时间参数，ti为爆炸荷载作用时长，Me为弹性阶段等效结构质量，Ce为弹性阶段等效结构阻尼，Ke为弹性阶段等效结构刚度，
Figure BDA0002786490650000022
为柔性结构等效体系振动加速度，
Figure BDA0002786490650000023
为柔性结构等效体系振动速度，y为柔性结构等效体系振动位移，ΔPe(t)为柔性结构承受的随时间t变化的爆炸动荷载，等效结构系数计算公式分别为：
Figure BDA0002786490650000024
其中，m为真实结构每延米质量，l为真实结构跨长，ξ为真实结构阻尼比，K为真实结构刚度，kM为弹性阶段质量变换系数，kL为弹性阶段荷载变换系数；由于爆炸冲击荷载持续时间非常短，可简化为等冲量的线性荷载，我国防护工程规范推荐采用的爆炸荷载为：
Figure BDA0002786490650000025
其中，ti为爆炸荷载作用时长，Δpm为爆炸荷载超压峰值，结构承受爆炸荷载之前初始位移、初速度均为0，求解该微分方程后，可确定此阶段位移和速度表达式为：
Figure BDA0002786490650000031
Figure BDA0002786490650000032
其中，无阻尼自振频率ω、含阻尼自振频率ωd、爆炸荷载超压峰值Δpm作为静载时对应的静位移yst各参数计算如下：
Figure BDA0002786490650000033
在爆炸荷载作用结束卸载的ti时刻，对应的位移和速度为
Figure BDA0002786490650000034
Figure BDA0002786490650000035
b、弹性阶段自由振动
由于进行设计的结构类型为柔性结构抗爆设计，当爆炸荷载消去，结构仍为弹性状态，此时结构为无外荷载、以位移yi及速度vi为初始条件的含阻尼弹性阶段自由振动，即当ti＜t≤tT时，等效体系的运动微分方程为
Figure BDA0002786490650000036
tT为柔性结构完成弹性振动，即将进入塑性振动的临界时刻，求解该方程后，位移及速度解答为：
Figure BDA0002786490650000037
Figure BDA0002786490650000041
其中将公式(7)、(8)代入公式(10)、(11)后，且令
Figure BDA0002786490650000042
则在tT时刻，结构弹性振动达到最大位移，此时位移、速度分别为：
Figure BDA0002786490650000043
Figure BDA0002786490650000044
c、塑性阶段自由振动
当结构振动时间大于tT时刻，为无外荷载、以yT及vT为初始条件的含阻尼塑性阶段自由振动，在tm时刻，结构振动达到最大位移，即当tT＜t≤tm时，动力等效体系的振动方程为
Figure BDA0002786490650000045
其中，me为塑性阶段等效结构质量，ce为塑性阶段等效结构阻尼，qm为结构最大抗力，其计算公式为：
Figure BDA0002786490650000046
其中，km为塑性阶段质量变换系数，kl为塑性阶段荷载变换系数，解(16)方程，求出此阶段位移和速度解为：
Figure BDA0002786490650000051
Figure BDA0002786490650000052
d、弹塑性阶段基于动力系数的延性比
当结构振动到最大位移ym时，对应的时刻为tm，此时速度vm＝0，代入(18)式，则：
Figure BDA0002786490650000053
将tm带入到(17)中得出结构弹塑性振动最大位移为：
Figure BDA0002786490650000054
如果令
Figure BDA0002786490650000055
将(21)、(16)代入(20)后，转变为：
Figure BDA0002786490650000056
由弹塑性阶段的抗爆设计，可知抗力动力系数kh和延性比β分别为
Figure BDA0002786490650000061
将(13)、(22)带入延性比β公式(23)中，可到
Figure BDA0002786490650000062
与现有技术相比，本发明具有以下技术效果：本发明根据实际情况，充分考虑结构的类型及结构阻尼对爆炸荷载动力系数的影响，尽量使设计的建筑结构更符合实际，在降低成本的同时，满足抗爆要求。并且通过该方法能够实现实际结构的精准化设计，也为抗爆设计奠定基础。
具体实施方式
为了使本发明所要解决的技术问题、技术方案及有益效果更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
含阻尼柔性结构抗线性爆炸荷载设计动力系数方法，所述的柔性结构指的是：在爆炸荷载作用时长ti范围内，结构振动未达到弹性振动最大值，爆炸荷载卸载后，依靠惯性力在tT时刻达到弹性位移最大值yT后，进一步发展为塑性振动，在某一时刻tm，达到了结构总的弹塑性位移最大值ym，
根据爆炸对建筑结构的作用过程，将该过程分为弹性阶段强迫振动、弹性阶段自由振动和塑性阶段自由振动三个阶段；
a、弹性阶段强迫振动
在弹性阶段且在荷载作用时长范围0＜t≤ti内，动力等效体系的运动微分方程为：
Figure BDA0002786490650000071
其中，t为柔性结构爆炸作用下的时间参数，ti为爆炸荷载作用时长，Me为弹性阶段等效结构质量，Ce为弹性阶段等效结构阻尼，Ke为弹性阶段等效结构刚度，
Figure BDA0002786490650000072
为柔性结构等效体系振动加速度，
Figure BDA0002786490650000073
为柔性结构等效体系振动速度，y为柔性结构等效体系振动位移，ΔPe(t)为柔性结构承受的随时间t变化的爆炸动荷载，等效结构系数计算公式分别为：
Figure BDA0002786490650000074
其中，m为真实结构每延米质量，l为真实结构跨长，ξ为真实结构阻尼比，K为真实结构刚度，kM为弹性阶段质量变换系数，kL为弹性阶段荷载变换系数；由于爆炸冲击荷载持续时间非常短，可简化为等冲量的线性荷载，我国防护工程规范推荐采用的爆炸荷载为：
Figure BDA0002786490650000075
其中，ti为爆炸荷载作用时长，Δpm为爆炸荷载超压峰值，结构承受爆炸荷载之前初始位移、初速度均为0，求解该微分方程后，可确定此阶段位移和速度表达式为：
Figure BDA0002786490650000076
Figure BDA0002786490650000081
其中，无阻尼自振频率ω、含阻尼自振频率ωd、爆炸荷载超压峰值Δpm作为静载时对应的静位移yst各参数计算如下：
Figure BDA0002786490650000082
在爆炸荷载作用结束卸载的ti时刻，对应的位移和速度为
Figure BDA0002786490650000083
Figure BDA0002786490650000084
b、弹性阶段自由振动
由于进行设计的结构类型为柔性结构抗爆设计，当爆炸荷载消去，结构仍为弹性状态，此时结构为无外荷载、以位移yi及速度vi为初始条件的含阻尼弹性阶段自由振动，即当ti＜t≤tT时，等效体系的运动微分方程为
Figure BDA0002786490650000085
tT为柔性结构完成弹性振动，即将进入塑性振动的临界时刻，求解该方程后，位移及速度解答为：
Figure BDA0002786490650000086
Figure BDA0002786490650000091
其中将公式(7)、(8)代入公式(10)、(11)后，且令
Figure BDA0002786490650000092
则在tT时刻，结构弹性振动达到最大位移，此时位移、速度分别为：
Figure BDA0002786490650000093
Figure BDA0002786490650000094
c、塑性阶段自由振动
当结构振动时间大于tT时刻，为无外荷载、以yT及vT为初始条件的含阻尼塑性阶段自由振动，在tm时刻，结构振动达到最大位移，即当tT＜t≤tm时，动力等效体系的振动方程为
Figure BDA0002786490650000095
其中，me为塑性阶段等效结构质量，ce为塑性阶段等效结构阻尼，qm为结构最大抗力，其计算公式为：
Figure BDA0002786490650000096
其中，km为塑性阶段质量变换系数，kl为塑性阶段荷载变换系数，解(16)方程，求出此阶段位移和速度解为：
Figure BDA0002786490650000101
Figure BDA0002786490650000102
d、弹塑性阶段基于动力系数的延性比
当结构振动到最大位移ym时，对应的时刻为tm，此时速度vm＝0，代入(18)式，则：
Figure BDA0002786490650000103
将tm带入到(17)中得出结构弹塑性振动最大位移为：
Figure BDA0002786490650000104
如果令
Figure BDA0002786490650000105
将(21)、(16)代入(20)后，转变为：
Figure BDA0002786490650000106
由弹塑性阶段的抗爆设计，可知抗力动力系数kh和延性比β分别为
Figure BDA0002786490650000111
将(13)、(22)带入延性比β公式(23)中，可到
Figure BDA0002786490650000112
通过上述方法，结构实际抗爆设计进行如下举例。
1、一种人防结构进行抗爆设计时，要求结构设计为柔性结构，该结构自振圆频率ω与爆炸动荷载作用时长ti乘积ωti为1.0，延性比β为1.60，阻尼比ξ为1％，km-l与km-l取值分别为0.66、0.78(《地下防护结构》方秦，柳锦春编著，ISBN 9787508470009中国水利水电出版社，2010年出版)，其它参数分别由其它公式计算得出，由柔性结构定义对θT(即ωtT)进行进行θT>θi(即ωtT>ωti)限定后假设初值后代入公式(24)，由延性比β的具体数值，采用迭代方法求出抗力动力系数kh为0.32。
2、一种人防结构进行抗爆设计时，要求结构设计为柔性结构，该结构自振圆频率ω与爆炸动荷载作用时长ti乘积ωti为0.8，延性比β为1.40，阻尼比ξ为5％，km-l与km-l取值分别为0.66、0.78(《地下防护结构》方秦，柳锦春编著，ISBN 9787508470009中国水利水电出版社，2010年出版)，其它参数分别由其它公式计算得出，由柔性结构定义对θT(即ωtT)进行θT>θi(即ωtT>ωti)限定后假设初值后代入公式(24)，由延性比β的具体数值，采用迭代方法求出抗力动力系数kh为0.27。
3、一种人防结构进行抗爆设计时，要求结构设计为柔性结构，该结构自振圆频率ω与爆炸动荷载作用时长ti乘积ωti为1.0，延性比β为1.20，阻尼比ξ为10％，km-l与km-l取值分别为0.66、0.78(《地下防护结构》方秦，柳锦春编著，ISBN 9787508470009中国水利水电出版社，2010年出版)，其它参数分别由其它公式计算得出，由柔性结构定义对θT(即ωtT)进行进行θT>θi(即ωtT>ωti)限定后假设初值后代入公式(24)，由延性比β的具体数值，采用迭代方法求出抗力动力系数kh为0.35。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包在本发明范围内。

Claims (1)

1.含阻尼柔性结构抗线性爆炸荷载设计动力系数方法，其特征在于：所述的柔性结构指的是：在爆炸荷载作用时长ti范围内，结构振动未达到弹性振动最大值，爆炸荷载卸载后，依靠惯性力在tT时刻达到弹性位移最大值yT后，进一步发展为塑性振动，在某一时刻tm，达到了结构总的弹塑性位移最大值ym，
根据爆炸对建筑结构的作用过程，将该过程分为弹性阶段强迫振动、弹性阶段自由振动和塑性阶段自由振动三个阶段；
a、弹性阶段强迫振动
在弹性阶段且在荷载作用时长范围0＜t≤ti内，动力等效体系的运动微分方程为：
Figure FDA0002786490640000011
其中，t为柔性结构爆炸作用下的时间参数，ti为爆炸荷载作用时长，Me为弹性阶段等效结构质量，Ce为弹性阶段等效结构阻尼，Ke为弹性阶段等效结构刚度，
Figure FDA0002786490640000012
为柔性结构等效体系振动加速度，
Figure FDA0002786490640000013
为柔性结构等效体系振动速度，y为柔性结构等效体系振动位移，ΔPe(t)为柔性结构承受的随时间t变化的爆炸动荷载，等效结构系数计算公式分别为：
Figure FDA0002786490640000014
其中，m为真实结构每延米质量，l为真实结构跨长，ξ为真实结构阻尼比，K为真实结构刚度，kM为弹性阶段质量变换系数，kL为弹性阶段荷载变换系数；由于爆炸冲击荷载持续时间非常短，可简化为等冲量的线性荷载，我国防护工程规范推荐采用的爆炸荷载为：
Figure FDA0002786490640000021
其中，ti为爆炸荷载作用时长，Δpm为爆炸荷载超压峰值，结构承受爆炸荷载之前初始位移、初速度均为0，求解该微分方程后，可确定此阶段位移和速度表达式为：
Figure FDA0002786490640000022
Figure FDA0002786490640000023
其中，无阻尼自振频率ω、含阻尼自振频率ωd、爆炸荷载超压峰值Δpm作为静载时对应的静位移yst各参数计算如下：
Figure FDA0002786490640000024
在爆炸荷载作用结束卸载的ti时刻，对应的位移和速度为
Figure FDA0002786490640000025
Figure FDA0002786490640000026
b、弹性阶段自由振动
由于进行设计的结构类型为柔性结构抗爆设计，当爆炸荷载消去，结构仍为弹性状态，此时结构为无外荷载、以位移yi及速度vi为初始条件的含阻尼弹性阶段自由振动，即当ti＜t≤tT时，等效体系的运动微分方程为
Figure FDA0002786490640000031
tT为柔性结构完成弹性振动，即将进入塑性振动的临界时刻，求解该方程后，位移及速度解答为：
Figure FDA0002786490640000032
Figure FDA0002786490640000033
其中将公式(7)、(8)代入公式(10)、(11)后，且令
Figure FDA0002786490640000034
则在tT时刻，结构弹性振动达到最大位移，此时位移、速度分别为：
Figure FDA0002786490640000035
Figure FDA0002786490640000036
c、塑性阶段自由振动
当结构振动时间大于tT时刻，为无外荷载、以yT及vT为初始条件的含阻尼塑性阶段自由振动，在tm时刻，结构振动达到最大位移，即当tT＜t≤tm时，动力等效体系的振动方程为
Figure FDA0002786490640000041
其中，me为塑性阶段等效结构质量，ce为塑性阶段等效结构阻尼，qm为结构最大抗力，其计算公式为：
Figure FDA0002786490640000042
其中，km为塑性阶段质量变换系数，kl为塑性阶段荷载变换系数，解(16)方程，求出此阶段位移和速度解为：
Figure FDA0002786490640000043
Figure FDA0002786490640000044
d、弹塑性阶段基于动力系数的延性比
当结构振动到最大位移ym时，对应的时刻为tm，此时速度vm＝0，代入(18)式，则：
Figure FDA0002786490640000045
将tm带入到(17)中得出结构弹塑性振动最大位移为：
Figure FDA0002786490640000046
如果令
Figure FDA0002786490640000047
将(21)、(16)代入(20)后，转变为：
Figure FDA0002786490640000051
由弹塑性阶段的抗爆设计，可知抗力动力系数kh和延性比β分别为
Figure FDA0002786490640000052
将(13)、(22)带入延性比β公式(23)中，可到
Figure FDA0002786490640000053
CN202011299957.2A 2020-11-19 2020-11-19 含阻尼柔性结构抗线性爆炸荷载设计动力系数方法 Withdrawn CN112395676A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011299957.2A CN112395676A (zh) 2020-11-19 2020-11-19 含阻尼柔性结构抗线性爆炸荷载设计动力系数方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011299957.2A CN112395676A (zh) 2020-11-19 2020-11-19 含阻尼柔性结构抗线性爆炸荷载设计动力系数方法

Publications (1)

Publication Number Publication Date
CN112395676A true CN112395676A (zh) 2021-02-23

Family

ID=74607424

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011299957.2A Withdrawn CN112395676A (zh) 2020-11-19 2020-11-19 含阻尼柔性结构抗线性爆炸荷载设计动力系数方法

Country Status (1)

Country Link
CN (1) CN112395676A (zh)

Similar Documents

Publication Publication Date Title
CN103473386B (zh) 一种确定水平运动下击暴流风剖面的方法
CN108416092B (zh) 起波配筋rc梁爆炸作用等效静载荷确定方法
CN109359265B (zh) 抽水蓄能调压室水流惯性时间常数允许值确定方法及装置
CN112395678B (zh) 含阻尼柔性刚性临界结构抗线性爆炸设计动力系数方法
CN112364425A (zh) 含阻尼刚性结构抗线性爆炸荷载设计动力系数方法
CN112395676A (zh) 含阻尼柔性结构抗线性爆炸荷载设计动力系数方法
CN114329736B (zh) 爆炸荷载作用下低阻尼刚性梁构件残余变形的求解方法
CN114329735A (zh) 爆炸荷载作用下中阻尼柔性梁构件残余变形的求解方法
CN112417559A (zh) 含阻尼柔性结构抗指数型爆炸荷载设计动力系数方法
CN112270032A (zh) 含阻尼柔性刚性临界结构抗指数型爆炸荷载设计动力系数方法
CN112417561A (zh) 含阻尼刚性结构抗指数型爆炸荷载设计动力系数方法
Forrai et al. Structural control technology: system identification and control of flexible structures
CN110990910A (zh) 时程激励下线性耗能结构响应的快速迭代法
CN114091207B (zh) 一种冲击载荷下的隔冲器时变可靠性分析方法
CN114329737B (zh) 爆炸荷载作用下中阻尼刚性梁构件残余变形的求解方法
Lu et al. An analytical solution for dynamic response of water barrier subjected to strong shock waves caused by an underwater explosion to dams
CN114329733A (zh) 爆炸荷载作用下高阻尼柔性梁构件残余变形的求解方法
CN114329734A (zh) 爆炸荷载作用下低阻尼柔性梁构件残余变形的求解方法
Wang et al. Nonlinear dynamic model of a turbine blade considering vibration and crack coupling
Shiau et al. Dynamic response of a spinning Timoshenko beam with general boundary conditions under a moving skew force using global assumed mode method
Li et al. Vibration Control of a Jacket Offshore Wind Turbine Under Earthquake Wind and Wave Loads by Tuned Mass Damper
Zhang et al. Optimization and evaluation of tuned inerter-based dampers for mitigating coupled responses of offshore wind turbines
Tombari et al. A Novel Reduced Column Approach for the Mitigation of Earthquake-Induced Vibrations of Wind Turbines
Song et al. Comparative analysis on dynamic response of wind turbine tower in different sites
CN116305606A (zh) 一种调谐惯容电涡流阻尼器的参数优化设计方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
WW01 Invention patent application withdrawn after publication

Application publication date: 20210223

WW01 Invention patent application withdrawn after publication