选择性特征热力学模块

1、1.热力学选择项目1.1 介绍FLAC3D 的热选择在材料和发展里允许瞬态热传导的模拟及热引起的替代和压力。这种选择有如下内容具体特征。1.2.3.4.5.6.7.8.两个热物质模型是到的：非线性传导和无效的模型。如同在 FLAC3D 的标准版本里，不同的区域可能有不同的模型和特性。任何机械模型可以被用于热模型。.温度，流，对流和绝热应变的边界条件可能被规定。.热源可能作为点源或者体积源入材料中。这些来源可以随着时间成指数.明确和暗示解题算法都可提供。.热选择为单机械压力和连合作准备并通过热膨胀系数的计算。.温度可以被通过 Fish 进入用户确定与温度有关的特性。这个部分包含热公式和 FLAC

2、3D 输入命令为说明热分析。提供适合解决热问题的建议，并提供几个证实问题。这些例子为建立 FLAC3D 模型热分析和连接热应力或者热易分析起到向导作用。1.2 数学模型的涉及在 FLAC3D 里的热传导的变量是温度和 3 个组成热流部分。这些变量因能量平衡方程式和 Fourier 理论推倒的热运移法律相互关联。Fourier 定律代能量平衡均衡弯曲到给定边界和初始条件的几何方程和热传导的微分方程。热体积应变被引入机械组成中形成热机械耦合。1.2. 1定义作为一次表达式，符号 ai 指示矢量a在一个指定系统里参考坐标中的组成部分 I， Aij是丈量A.的组成部分(i，j)。此外，f，i量变量或者

3、一个矢量组成部分.)。与 xi 有关的 f 的部分。 (f 可能是一个标求和只在指标 i，j 和 k 上使用，这起到 1，2，3 与空间尺寸有关的组成部分作用。当使用在基本方，指标带有任何方位。SI用来说明变量的参数和程度。 关于变换请参阅第 1.6 部分其他制。 下列尺度数目在瞬态热传导的描述过特性长度：起到作用。固体的体积L （1.1）c表面可以传递的热量热扩散：热传导系数。 密度。为体积的比热。特性时间：1.2.2 热平衡方程：能量平衡：通过面积上的热流量， 为通过体积的热流量，体积的热能量，通常，温度变化可能由在能量及参数可能被以表示：和体积应变方面的变化引起， ，以及热组成规律涉和是

4、物性参数，T 是温度。在 FLAC3D 里，认为和。 密度，假说是应力变化影响温度中作用可以忽略，类似的有效的假定。因此，的机械问题与固体和液体有关把公式（1.6）带入到公式（1.4）中可以得到：考虑到对于固体和液体，在恒定的压力作用基本上相等。 从而，Cv 和 Cp 被可互换使用。1.2.3 传递规律积基本不变，比热注意到那 体积Fourier s 定律是确定温度热流量之间的关系，对于稳定各向同性的固体可以表示为：T 温度，k 热传导系数。1.2.4 边界条件和初始条件把公式 (1.8)带入公式 (1.7)中到为热传导产生微分方程。 最初 条件符合规定的温度场。 边界条件一般被表示在里 温度

5、的条件或者热流矢量与边界垂直的组成部分。 在这个版本里 在 FLAC3D 中，4 类条件被考虑，相当于：(1)给温度； (2)给流与边界垂直的组成部分； (3)对流的边界； (4)在 FLAC3D 里，给定边界(绝热应变)边界。与边界垂直朝外表方向的流量， 对流的热传导系数， 表明温度， 流体表明的温度，1.2.5 热机械耦合解决办法热应力问题需要对的压力劳累关系的公式更改， 减去由于温度变化而引起的热应力变化。因为热变化不会导致变形。热应力的变化和温度变化的关系：线性热膨胀系数，基，以机械为基础运动速度和应力关系的比率的微分方程，不变适合热机械耦合模型。1.2.6 热流耦合引起的毛力通过考虑

6、由于热膨胀应力孔引起孔隙压力的的变化而吧热传输应用到水的计算中。运移 形式(看见公式(1.11)在流体力学的相互作用过)流体能量变化量， 热膨胀系数，在 FLAC3D 里，由于温度引起孔隙压力的变化。在没有变形的情况下由于温度引起的孔隙压力的变化量。对于理想多孔材料，这个系数是体积热膨胀系数的梯度。孔隙率。对于不可压缩材料，孔隙压力为：流体莫数。在连续方考虑总应力变化1.3 数学模型FLAC3D 中基于组成的区域介质接近四面体，平衡方程式 (1.7),以及 Fourier s 法律(1.8)使用一个有限差分，并依据能量方程控制数值。这个机械公式导出定律。通过以温度，热流矢量和温度梯度，速度，压

7、力，疲劳矢量迭代获得。根据显式和隐式，普通的微分方程组可以通过离散方式来解答。原理推倒如下。1.3. 1 有限差分 近按照惯例，n 表示四面体的四个节点数目，和上角标的(f)涉及价值达在表面 f 上的相关变量。在一个四面体内认为是线的温度变化； 温度梯度表示为，L 面法线方向， 外表面积， 热体积流量。1.3. 2 能量平衡方程能量平衡方程表示为：使用机械节点公式。首先考虑了单个的四面体，到考虑热流量和体力的情况下为：每一个单元上建立能量平衡方程，在每一个单元上，从某个上流出热量等于从其余个节点上流量的和，认为是。在(1.20) 与温度梯度有关四面体热流矢量的组成部分。(1.8)温度梯度有关的

8、运移方程。(1.16). 四面体节点温度梯度。为了节省计算时间，在 FLAC3D 里采用采用区域公式。在每一个节点 n，总流量 是属于该节点通过全部四面体的区域和除 2 得到。区域矩阵 M是在区域温度 T n(n = 1，8)时的流量。为节省时间这些矩阵是有 36 个项， 每 10 步计算和更新。温度矩阵，依次，全球节的流量通过多次迭代获得。对于表达式来说，是所有节点求和。应用公式. (1.19)带入公式(1.24)到,节点温度能量。总体积流量。在公式1.25中，表示在n点的能量平衡方程。能量输出平衡。在FLAC3D对于一个普通的微分方程的系统使用显式和隐式有限差分格式过应用。1.3.3 显式

9、有限差分公式解决。下面在显式有限差分公式，温度在一个节点被假设为线性变。在(1.25)左边向前使用有限差分表示，并且在时间 t 在热平衡输出。从最初温度场开始，温度随时间变化不断更新，如果温度不被确定，使用表示为：只有时间步限定到一定范围内才能保证数值稳定性。1.3.3.1 稳定标准为了得到稳定性标准，T0。考虑给定某些区域给定给在公式。 (1.23),温度从一种零初对于对流边界，为：在 n 节点在全球对流的时期温度系数，通过热时间迭代，新的温度为：为了保证那些计算被重复适合连续)， 在上述关系里的 T0 的系数一定是确定的。要求时间遵循：在 FLAC3D 使用的时间步在最小节价值的这右边在方

10、面 方程 (1.33),乘以 0.8 的安全系数。为了评价有关参数的影响，记住如下内容是有用的 关键的时间步的代表。假如最小的四面体特性长度，到：标准时间步，是特性参数，m 依据几何形状是常量，在公式(1.34)里需要一个要在整个四面体的前面扩散，符合的特性时间的度量标准的时间步。如果 Lc，使用表示相似特性长度，需要估计的整个时间。不考虑对流情况。给定特性参数，那么在给定四面体时间是和边长成比例。如果数目很大，这种现实迭代不被采用。另一方面，这种第一个命令的方法的优势是给定的计算的时间步通常是小到足以满足温度随时间准确波动。1.3.4 隐式迭代模拟在多个层都是当瞬态，必须在尺寸方面限制到某一

11、特定值才保证在极小的时间系统是稳定的。隐式的公式消除这个限制，但是包含在每个时间步解决联立方在 FLAC3D 的这隐式公式使用这曲柄方法，假设节点温度随间隔 t 改变而变化。在此，引入中心差分等于在平衡热传递在 t 时刻和 t t 时间传递量的平均值。并公式如下：热量选择由定义（1.23）可以得到：和由公式（1.27）,最后的表达式为：把公式（1.40）带入（1.36）用公式（1.38）可以得到：定义已知的矩阵A和垂直 bt为可以把公式（1.42）写成如下形式：这样等式也许被写为每个结被介入在栅格。 收效的系统等式被解决在FLAC3D 运用Jacobi.s 重申方法。 在这种方法, 温度增加在

12、叠代 r + 1 和结 n 被计算使用再现联系网格中的每一个节点可以用这样的等式来表示，在FLAC3D中应用Jacobis iterative方法解决系统等式。应用这种方法，应用下边的关系计算在(1.43) (1.46)中 A 使用指标 j到初始近为：因此初始接近使用相似差分格式。是显式格式/。在 FLAC3D 内，考虑 3500 迭代，那标准形式：在隐式计算时间步大小要准确联系到隐式的稳定性。必须A满足如下公式才能保证柄方法是全部时间 t 内的稳定 。根据公式Eq. (1.43)定义的对于和公式定义的 必须遵循时间间隔 。如果雅各比无法获得，命令就会出息错误。此时必须降低时间步 ，才能保证正

13、常运行。此外，虽然隐含法准确的二次命令，但是一些选择方面需要合适的时间步。因此它的值必须比节点温度波动的波长小。在 FLAC3D 里选择这种隐式的方法，一定考虑计算时间和计算机器二。在这种隐含法方面，在每时间步需要至少迭代 3 次。计算的数量被要求 一个反复的台阶大约等于在明确计划里为一时间步需要的数值。 此外，中间值必须被在反复的程序里，需要额外的跟明确计划相互分配。那些不便，隐式计算中希望获得较大的时间步。1.3.5 热应力耦合瞬时热传递中热传输随时热应力计算耦合。耦合朝一个方向出现，即，温度可以导致压力变化，机械变化引起物体温度变化。在稳定的机械运动中能量的改变量被忽略。在公式(1.10

14、) 三角形区域应力改变。上述以为恒温三角形区域内在1.3.6 热孔隙压力耦合在公式 (1.34)由于温度变化引起毛栅栏点。力变化。结合公式(1.11),(1.12)以及(1.15)计算总应力方面的相应变化到：对于不可压缩固体颗粒，水特性也能通过温度变化而引起特性变化。1.4单一热问题和耦合热问题在Fl中可以解决单一热问题和机械热耦合问题。耦合变形在1.4.2 和 1.4.3部分。第一步在命令热的配置热计算界面。然后，必须选择明确和暗示热算法。选择确切算法，在热的任何阶段暗示的方式可能被开动，使用那些命令SET thermal implicit on or SET thermal implici

15、t off“.。在确定明确方式后，热时间步将被自动计算，但是更小的时间步可以通过命令SET thermal dt来选择。必须选定时间步的大小。这通过命令SET thermal dt，热模型和特性必须指定。指定密度，最初和边界条件，流体边界条件。1.4. 1单一热分析“SET mech off”命令确定单一热计算。“STEP”确定热计算步。当到达规定时间步，计算自动结束。“SET thermal age or SET thermal step”用来确定热计算时间步的范围。 “PROPERTY conductivity and PROPERTY”确定材料热物理参数。稳态热传导问题可以用 SOLVE

16、 命令。在某节点热不平衡传导中热传导梯度小于10-5范围内，计算才能正常运行。SOLVE ratio (or SET ratio)限制稳态热传导比率。STEP or SOLVE确定时间步和总计算时间。PRthermal确定时间步和其它热物理状态。1.4.2热机械耦合FLAC3D的热机械分析，热选择可能与机械计算相结合。 在热机械计算里，热计算可以选择瞬时和稳态热传输和热解决办法。在公式(1.10)中温度变化对体积的影响。选择热膨胀系数。在FLAC3D里机械模型里可以进行热机械分析，热模型中可以机械分析。M给定材料热区域。th null分析两个阶段要考虑一热机械分析。 第一个涉及与热相关的不同的

17、时标和机械过程。因为时间步相当于时间为信息从一个节点赋值到下一步，公式用来能比较时标。通常，与机械过程相关的t时间步具体形式：体积模数；剪切模量 密度标准长度。热机械时间步比率为：热扩散系数，对于非金属，取值一般为，对于岩石和固体， 密度取值为，根据式(1.55),热机械时间步比率一般为。近似为这比率保持非常大， 即使Lc 实际有1毫米，在FLAC3D内假设机械效应与扩散影响相比能可以瞬间出现。同样在考虑热传递中步考虑任何机械运动。在FLAC3D里第2 个问题关系到热机械耦合。当热计算不会影响机械变化时，耦合只出现朝一个方向，温度变化引起热应力变化。在大多数模型化过，最初在耦合的分析开始之前，

18、首先取得机械状况符合一个平衡状态。应用命令(SET thermal on mech off).获得介质弹性和热机械性能。SET thermal ond SET mech on表示热计算关闭，机械计算正在运行。这样系统循环进行运算。对非线性的机械模型(即，塑性体)来说，热变化一定 在更近的时间间隔传递遵循机械模数属性。 通常，热时间步很小的情况下， 在每步热计算中包含若干机械计算。在大的热时间范围，如果系统分析法热平衡，步扰乱介质的机械状态情况下，可以采用几个热时间步。通过热和机械的方式之间交替数字模拟控制。SET thermal on mech on控制热机械运行。SET mechRatio控

19、制机械平衡比率。SETmech substep n slave确定控制机械模数。SET thermal substep m控制热模数。如果，每热时间步，机械模件只需要一种附属措施到达均衡， 然后热附属台阶的的数量是时间机械附属步的二倍，但是它从未超过数m。无论附属台阶机械时间步连续或中止，系统恢复的原先的设置。下列信息在将打印到屏幕上当命令循环时。总时间步，模数运行，现在运行步，最大不平衡机械力率，最大热平衡利率，现在运行时间模数。这信息可以被保证计算正常进行。在第3 种方法，(STEP)命令中(SET thermal on mech on).命令确定热机械同时运行。每步热运行将会同时机械运行

20、。在选择了热时间步后，必须保证在热时间步足够机械运行的平衡。1. .1部分将叙述变量选择。1.4.3热孔隙应力耦合保证FLAC3D也能用来解决由于温度变化引起体孔隙压力的耦合问题。这个耦合也可以考虑机械压力。那些耦合因为这温度变化对体积变化的影响引起饱和的基体体积变化。 和积的热膨胀系数，可以通过令PROP RTY u thc来确定。CONFIG flui therma ).令确定热流体耦合题。 ET hermal n f uidof or SET thermal oluid on确定是否运行热流体耦合。只有确定小的时间步，才能保证热和流体时间步彼此协调。提供两个简单的例子说明热流体耦合的机械

21、流体反应。两个例子都是饱和态的例子，在第一个例子里，无拘束。第二个都限定。在两个例子中温度都因热源而增变化。对无拘束的验来说，体积增大，毛压力减少。 总平均的压力持 在零点约束。拘束的情况的体应变：孔隙压力常系数 常模数线性热胀系数 饱和区热膨胀系数排干水体积模数温度变化。在被限制的实例中，体积变化是零，毛在这个情况里增力由于温度增加增加。因为边孔隙应力变密度 (体积模) 1(K 0. 7两个例参如下：剪切模数 G) .孔隙率 (n 0. iot模数 ( ) 4i t线热系数 () 1胀系数 t) 3非排热膨胀系) 5热传导k)比热 (Cp) 1温度变化 ()T ) 0.4对于约束从公式(1.

22、56) 和 (1.57),分别选2.230和-0.918。例子1titleUndrained thermo-poro-elastic expan config fl thermal; - fish constants -def setuptest - unconfinedtestc_biom c_bioc c_poro c_cond c_cv =c_thex=4.1.; biot modulus; biot coefficient (default value)0.5 ; porosity1.; thermal conductivity specific heat; linear therma

23、l expan1. ;= 3.coefficientc_u_thc = 5. ; undrained thermal coefficient (beta)c_dens = 1.; drymedium density modulusmodulusc_k = 0.67 ; bulkc_g = 0.4 ; endsetupsheargen zone brick size11 1 p1 1 0 0 p2-0 1 0p3 0 0 1; -m prmechanical melu c_ksh c_g c_densflow minidensity; - fluid-mfl_isoset fl biot on;

24、 default value: poros ini biot_mod c_biom ini ftens -1e100.5 biot_c 1; - thermal m-mth_isoprop cond c_condapply vsource 1.sp c_cv thexp c_thexrange z -1 2u_thcc_u_thc; - fish function - def ana_soldelta_t = gp_temp(gp_head)v_strain_t = 3.0 * c_thex * delta_tterm_1 = term_2 =v_strainc_bioc * c_biom *

25、 c_u_thc * delta_t + c_k * c_k + c_bioc2 * c_biom= term_1 / term_2v_strain_tdelta_p = c_biom * (- c_bioc * v_strain + c_u_thc * delta_t)pn t = gp_n ea r( 1 . 0 , 1 . 0 , 1 . 0 )e rr _vo l = 100 . 0 * ( 6 . 0 * gp_xd i s p ( pn t )- v_ s t r a i n )/v_ s t r a ie rr _pp = 100 . 0 * ( gp_pp ( pn t )-

26、e nd;-s e tti ng s- s e tf l o ff s e tt h e rma l d t 8 e - 4;-t e s t- s e t mec h o fft h e rma l on s t e p 500; p a u s ed e lt a _p )/d e lt a _ps e t mec h on t h e rma l o s t e p 500a n a _ s o ll u i do ff p r p r p r r e t2 titl ev_ s t r a i n d e lt a _pe rr _vo le rr _ppUnd r a i n e d

27、 t h e rmo - po r o - e l a s ti ct e s t- c on f i n e d t e s t c on f i g f lt h e rma l;-f i s h c on s t a n t s- d e fs e t upc _b i om c _b i o c c _po r o c _ c ond c _ c v =c _ t h e x=4 .; b i o t modu l u s1 .; b i o tc o e ff i c i e n t( d e f a u lt v a l u e ) 0 . 5 ; po r o s it y1 .

28、;t h e rma lc ondu c ti v it y1 .;= 3 .s p ec i f i c h ea t;li n ea rt h e rma le xp a nc o e ff i c i e n tc _u_ t h c = 5 .; und r a i n e d t h e rma lc o e ff i c i e n t( b e t a )c _d e n s c _k =c _g = e nds e t up= 1 .0 . 67 ; d r ybu l kme d i um d e n s it y modu l u smodu l u s0 . 4 ;s h

29、 ea rg e n z on e b r i c k s i ze11 1 p1 1 0 0 p2 0 1 0 p3 0- 0 1;- m p rmec h a n i ca l me lu c _ks h c _gc _d e n si n if i xd e n s it yxy zf l u i df l _ i s o;- mf l ow m- set fl biot on; default value: poros 0.5 biot_c 1 ini biot_mod c_biomini ftens -1e10; - thermal m-mth_isoprop cond c_cond

30、apply vsource 1.sp c_cv thexp c_thex u_thc c_u_thcrange z -1 2; - fish function - def ana_soldelta_t = gp_temp(gp_head)v_strain_t = 3.0 *c_thex * delta_t c_u_thc * delta_t* v_strain_t + c_bioc * delta_p)delta_p = delta_sxxerr_sxx =c_biom *= - (c_k100.0 * (z_sxx(zone_head) - delta_sxx) / delta_sxxerr

31、_pp = 100.0 end; - settings set fl offset thermal dt; - test -* (gp_pp(gp_head) - delta_p) / delta_p-8e-4set mech off thermal on step 500set mech on thermal off step 500ana_solpr pr prretdelta_sxx delta_perr_sxx err_pp1.53FlD命令下列命令被提供运转热问题。 注意应用热命令可以对热进行选择。在不同时间内用CONFIG thermal给定热命令。命令流是确定热边界条件和范围。如

32、果没有特殊指定，命令适用于整个模型。 可选择的关键字可以 并赋予其数值。热边界条件使用两种关键字类型。数值赋予根据类型确定。热源 是指定点和某网格范围。当一个新来源用于一个栅栏点时， 新来源强度替换原有源强度。来描述热源能量衰竭。热源 是体积热源和网格范围内能量。新来源强度替换原有源强度。来描述热源能量衰竭。温度， 介质发生对流。 热传导系数。对流的边界条件限定范围内使用。 历史关键不适合对流。初始流量，流体在指定范围内的流量。这命令用来指定对(v0，的恒定流)或者在之外(v0)节点以外的热边界，新流量强度替换原有强度。下面例子。def decay decay=exp(deconst*(tht

33、ime-thini) endset thini=0.0 deconst=-1.0apply flux=10 hist=decay这个命令运行特殊点的热分析，在fl中必须给定阶段才能进行热分析。某节点温度给定，温度就是固定不变。给定温度允许变化。在给定点的温度。在给定范围内给定初始温度。这个命令是某个范围内热传导。各向同性热传导区域内没有热传导，节点热分析。确定温度范围。确定材料属性。节点热变化。节点温度热计算给定区域内材料属性。M不同热特性模型。下面列举关键字确切说明特性。确定材料密度非线性热传递系数, k比热, Cv(3)(4)性膨胀系数 t非排水热系数, 确定热结构耦合参数。是否运行流体计

34、算。CONFIG fl id命令给定就会运行流体计算，否则相反。是否运行结构计算。On命令给定就会运行流体否则相反。在耦合计算中 定耦合结构计算步。结构计算部分会自动调运到热计算中，参照1.4.2部分。在热传递中给定热限制时间。给定热时间步，时间步必须特殊给定。FLAC3D自动计算热时间步对于没有明确说明。可以允许应用不同的时间步。F AC3D中如果用户选值太大时，为了数值稳定性，时间步将被降低。命令运行，计算将不恢复到用户选择的价热模型里开通或断开。 默认断开CONFIGthermal命令给定，热模型计算运行，否则只有流体计算。确定热流比率。为全部节点最大不平衡模型热流除以平均全部栅栏点热流

35、大小。比率低于估计值时计算将停止热时间步的限定在热结构耦合中时间步将被调整。给定，热计算时间步给定，时热结构耦合部分时间步。这命令自动控制热时间步和耦合的热结构计算。 计算执行直到限制条件，命令结束。热时间限制计算机运行时间1.5.2Fish变量下列标量变量时用来帮助热分析。热时间步热流量比率加热时间一个FISH功能修改下列FLAC3D 节点变量：节点温度通过Fish功能可以获得热特性变量。参照1.51部分热物理参数。1.6热分析所以的热物理必须给定，不能转换变量。1.1，1.2给定热物理参数。可以被通用替换如五用替换 K，oF 替换 R 替代关系如下：一些1.7 实例验证在FLAC3D中，几

36、个证明问题证明了热量模型。在“OptionsTherma ” 目录中存有这些实例的数据文。1.7平面层中的导厚度L = 1 m 的平面板料最初是在恒温 0，一边的温度达到 100同时另一边是 0。最后达到恒温流和不变的温度分。Cr19 5 年对瞬时温度进行了解1 板表面一边的T2 板表面另一边的度L 板的度 T间Z 穿过板的距离和k等同于k/(Cp),k是温度传到率；是密度Cp是比。当 T2=0 时，解这个例子的热传导率是 1.6W/m，比热是 0.2J/Kg，材料的最大密度是 1000Kg/m3,温度 T是 100 。分析解答利用编的 FISH 功能直接比较大量的有选择的热量段，这段时间的分

37、析和大量的温度结果保存在表格。在 FLAC3D 模型中，板料被定义为 5 个区域。位置在 Z=0 的表面的恒温边界为 100，在Z0 的位置表面的恒温边界为 0。模型网格划分如图 1. 。包含FLAC3D 数据文件的例子1.3为这个问题使用明确公式化获得解答。在例子1.4 包含数据文件中使用了不确定的公式。为了精确解答在三热量段的分析和数值温度，在表1.2中做了比较。不确定表达表1.3中。在两个图中绘制了标准化的温度(T/T1) 和标准距离(z/L)，在表2，4和6中给出了温度的分析解，在表1，3和5给出了FLAC3D解。三热量时间是1.455, 7.273 和72.73 秒为确定解, 1.4

38、55， 11.45和71.45秒为不确定解。在所有情况下问题的解答达到了平衡热量状态，就以上两种解决公式，分析和数值温度之间在稳定状态的区别小于0.1。mth_isoprop conducti 1.6 spec_heat 0.2inifix fixdensity 1000tt100 range z -.001 .0010 range z 0.999 1.1 fish constants -; -def consc_cond = 1.6 ; conductivity c_dens = 1000. ; densityc_sph = 0.2length = 1. t1 = 100. ; specif

39、ic heat; wall thicknesswall temperature, face 1t a bn = - 1 t a b e = 0ov e r l = 1 .d = c _ c ond/l e ng t h/( c _d e n s * c _ s ph ) ov e r l * ov e r lp ip i; ma x numb e r o ft e rms- e x ac ts o l u ti ondo l 2 = d*/*t op =p i 2 = n_ma x2 .p i= 100t e p s = 1 . e - 5 ;sma ll v a l u e c omp a

40、r e d t o 1 e ndd e f num_ s o lt a bn = t a bn +2* do l 2 p i 2t _h a t t p2 =pn t = t h time t _h a t *gp_h ea dl oop wh il e pn t# nu ll r a d = s q r t ( gp_x( pn t ) 2 + gp_yi fr a d 1 . e - 4 t h e n( pn t ) 2 )*ov e r lx = gp_ z( pn t ) * ov e r lt a b l e ( t a bn , x ) = gp_ t emp ( pn t )/

41、t 1 e nd_ i fpn t = gp_n e x t ( pn t ) e nd_ l oope ndd e f a n a _ s o lt a b e = t a b e +2* do l 2 p i 2t _h a t t p2 =pn t = t h time t _h a t *gp_h ea dl oop wh il e pn t# nu ll r a d = s q r t ( gp_x( pn t ) 2 i fr a d 1 . e - 4 t h e nx = gp_ z( pn t ) * ov e r ln = 0+ gp_y( pn t ) 2 )*ov e

42、r ln it = 0t s um = 0 . 0t s umo = 0 . 0 c onv e r g e = 0l oop wh il e n n_ma xn = n + 1f n = t e rmt s umf l o a t ( n )= s i n ( p i *x* f n ) * e xp (- t p2* f n* f n )/f n= t s umo + t e rmif tsum = tsumo nit = n table(tabe,x) = converge = 1n = n_max elsetsumo = tsum end_if end_loopif converge

43、= 0then1. - x - top * tsumthenii = out( not converged: x= + string(x) exitend_if end_ifpnt = gp_next(pnt) end_loopend+ t = + string(thtime); - settings set mech off set thermal on; - test - cyc 20cons num_solana_sol-plot add plo show cyc 80 num_sol ana_sol plot mod plo show cyc 900 num_sol ana_sol p

44、lot modplot showtable1 2 line black red black red black red black red13 4line15 6lineplo ploplocreatetempadd tab 2 4 6 1 lin 3 lin 5 lin bla bla bla blu gre cyanshowsave con_exp.sav例1.4平面薄板传导隐式解法titleConduction in a plane sheet: implicit method config thermalgen zone brick size 1 1 25 p1 0.1 0.0.p20

45、.0.10.p30.0.1.; - thermal m-mth_isoprop conducti 1.6spec_heat 0.2inifix fixdensity 1000tt100 range z-.001 .0010 range z 0.999 1.1fish constants -; -def consc_cond = 1.6 ; conductivity c_dens = 1000. ; densityc_sph = 0.2length = 1.; specific heat; wall thicknesswall temperature, face 1t1 = tabntabe10

46、0. ;= -1= 0overl = 1. d = c_cond dol2 = d *top = 2. / length/ (c_dens * c_sph) overl * overlpip i 2 = p i * p in_ma x = 100 ; ma x numb e r o ft e rms- e x ac ts o l u ti on t e p s = 1 . e - 5 ;sma ll v a l u e c omp a r e d t o 1e ndd e f num_ s o lt a bn = t a bn +2* do l 2 p i 2t _h a t t p2 =pn

47、 t = t h time t _h a t *gp_h ea dl oop wh il e pn t# nu ll r a d = s q r t ( gp_x( pn t ) 2 + gp_yi fr a d 1 . e - 4 t h e n( pn t ) 2 )*ov e r lx = gp_ z( pn t ) * ov e r lt a b l e ( t a bn , x ) = gp_ t emp ( pn t )/t 1 e nd_ i fpn t = gp_n e x t ( pn t ) e nd_ l oope ndd e f a n a _ s o lt a b e

48、 = t a b e +2* do l 2 p i 2t _h a t t p2 =pn t = t h time t _h a t *gp_h ea dl oop wh il e pn t# nu ll r a d = s q r t ( gp_x( pn t ) 2 i fr a d 1 . e - 4 t h e nx = gp_ z( pn t ) * ov e r l n = 0n it = 0t s um = 0 . 0t s umo = 0 . 0 c onv e r g e = 0l oop wh il e n n_ma xn = n + 1+ gp_y( pn t ) 2 )

49、*ov e r lf n = t e rmt s umf l o a t ( n )= s i n ( p i *x* f n ) * e xp (- t p2* f n* f n )= t s umo + t e rm/f ni ft s um = t s umo t h e n n it = nt a b l e ( t a b e , x ) = 1 .- x -t op * t s um c onv e r g e = 1n = n_ma x e l s et s umo = t s umend_if end_loopif converge = 0 thenii = out( not

50、converged: x= + string(x) + t = + string(thtime) exitend_if end_ifpnt = gp_next(pnt) end_loopend; - settings - set mech offset thermal on; - test -; - start with explicit method - set out lpt1cyc 20 cons num_solana_solplot add table 1 2 line black red black plo show; - then switch to implicit - set

51、thermal implicit on dt=1.e-1 cyc 100num_solana_solred black red black redplot mod 1 plot show cyc 600 num_sol ana_sol plot mod 1plot show3 4 line5 6 lineplo ploplocreatetempadd tab 2 4 6 1 lin 3 lin 5 lin blashowbla bla blu gre cyansave con_imp.sav1.7.2 空心圆筒的受热最初在恒温0一个无限长度的空心圆筒，在圆筒内径加热到100外径保持在0情况下，

52、研究在达到热平衡状态下，确定温度和热应力。Nowacki (1962)在稳定热量状态的条件下，根据温度和半径在切向和轴向的应力，对这个问题给予了解答。T 温度；r 径向应力r 圆筒中心到边的距离；a 内径；b 外径Ta 内径温度 ；t 切向应力a 轴向应力m3kat/2G=K-2/3 G K是体积弹性模量G 是剪应力模量和是线性热量膨胀系数温度和应力的分析解，在FLAC3D数据文件中以编程的方式，应用FISH功能实现。在表格中分析解和数值解可以直接对照。以下是这些例子的性质参数：内径a外径b原料性质密度1.0m2.0m2000kg/m3比热（Cp）880.0J/Kg.C热导率(K) 4.2W/

53、m.C-6线膨胀系数（at）5.410 /剪应力模量（G）28.0Gpa体积模量(K) 48.0 Gpa圆筒的四分之一在FLAC3D中被模拟，图1.4显示了FLAC3D中的分格情况。模型的内径恒温边界条件为100，外径为0。FLAC3D模型可以在耦合和非耦合热结构下运行。用热弹性分析非耦合是非常有效的。通过运行耦合和非耦合模型使这问题得到证明。先运行这非耦合模型：热计算首先执行到达热平衡状态；然后热感应力学应力被计算。在例子1.5中显示了FLAC3D数据文件。例1.5空心热管titleHeating of a hollow cylinder - thermal & mechanical cal

54、culations config thermaldef gcons ; geometrical constants c_b = 2.n1 = 10y1 = 1. / float(n1)eps y1m y1pend=1.e-4y1y1- eps+ epsgcons gen zonedimencsp1 c_b0 0 p2 0 y1 0 p3 0 0 c_b &n1 1 12-1 1 1 1 size; - mechanical mm p rf i x f i x f i xf i xe lu 48 e 9 r a ng e r a ng e r a ng e r a ng et h e rma l

55、s h 28 e 9z x yyz- . 0001- . 0001- . 0001y1m y1p. 0001. 0001. 0001x yy; - r a ng e r a ng e mmc yy- e nd1 0 - 1 0 e nd2 0 2 0 r 1 . 0001n amei nn ame o t h_ i s oe nd1 0 - 1 0e nd2 0 20 r1 . 9999no tp r op c ondu c ti 4 . 2t h e xp 5 . 4 e - 6s p ec _h ea t880i n i f i xf i xd e n s it y 2000tt100 r

56、 a ng ei n0 r a ng es e tti ng sou t- y s i s;- ; un c oup l e ds e t mec h o fft h e rma l ons o l v es e t mec h on s o l v es e tt h e rma l;t h e rma l o ff on;c oup l e dy s i son t h e rma l on r a ti o 1 e - 3s e t s e t s e ts e tmec hmec ht h e rma lr a ti o 1 e - 3mec hs ub s t e p 10000 a

57、 u t o;s l a v es o l v e-f i s hc on s t a n t s- d e f c _gc _kc on s= 28 e 9= 48 e 9;s h ea r modu l u s; bu l k modu l u sc _ a l = 5 . 4 e - 6 ;c o e ff i c i e n t o ft h e rma le xp ant 1 = 100 .;bound a r y t emp e r a t u r e l n ( c _b )( c _b * c _b - 1 . )( c _k - c _g * 2 ./ 3 . ) / ( c

58、 _k + c _go c 1o c 2o c 3=1 .1 .0 . 5/*/ 3 . )c _mmu t a b1t a b2 t a b3 t a b4t a b5= c _g * ( 3 . * c _k * c _ a l )/( c _k + 4 . * c _g /=12345;nume r i ca ly ti ca l nume r i ca ly ti ca lnume r i ca lt emp e r a t u r e t emp e r a t u r er a d i a ls t r e ss r a d i a ls t r e ss t a ng e n t

59、i a ls t r e ss t a b6 t a b7 t a b8 e ndc on s=678;y ti ca l nume r i ca ly ti ca lt a ng e n ti a ls t r e ss a x i a la x i a ls t r e ss s t r e ss d e f num_ s o lt pn t = gp_h ea dl oop wh il e pn t # nu ll r a d = s q r t ( gp_y( pn t ) 2 + i fr a d 1 . e - 4 t h e ngp_ z( pn t ) 2 )x = gp_x(

60、 pn t )t a b l e ( t a b1 , x ) = gp_ t emp ( pn t ) e nd_ i fpn t = gp_n e x t ( pn t ) e nd_ l oope ndd e f a n a _ s o lt nn = 0pn t = gp_h ea dl oop wh il e pn t # nu ll /t 1r a d = s q r t ( gp_y( pn t ) 2 + gp_ zi fr a d 1 . e - 4 t h e n( pn t ) 2 )x = gp_x( pn t )t a b l e ( t a b2 , x ) = l