lanxicy.com

第一范文网 文档专家

第一范文网 文档专家

International Communications in Heat and Mass Transfer 36 (2009) 143–147

Contents lists available at ScienceDirect

International Communications in Heat and Mass Transfer

j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / i c h m t

Optimizing the exponent in the heat balance and re?ned integral methods☆

T.G. Myers

Department of Mathematical Sciences, Korean Advanced Institute of Science and Technology, 335 Gwahangno Yuseong-gu, Daejeon 305-701, Republic of Korea

a r t i c l e

i n f o

a b s t r a c t

The most contentious aspect of the Heat Balance Integral Method is the choice of power of the highest order term in the approximating function. In this paper we develop a new method, where the exponent is determined during the solution process. This is achieved by minimising the square of the difference of the terms in the heat equation. The solution requires no knowledge of an exact solution and generally produces signi?cantly better results than previous models. The method is illustrated by applying it to three standard thermal problems. ? 2008 Elsevier Ltd. All rights reserved.

Available online 6 December 2008 Keywords: Heat Balance Integral Method Re?ned integral method Heat ?ow Polynomial approximation

1. Introduction The Heat Balance Integral Method (HBIM), developed by Goodman [3], is a simple approximation method for solving the heat equation. However, since the heat equation is amenable to exact solution by a number of techniques the HBIM has had the greatest impact on Stefan problems where very few analytical solutions exist. Recently it has also been applied to modelling the temperature in a thermistor [6], the process of contact melting [11], the Korteweg-de-Vries equation [10], the solidi?cation of a molten material sprayed onto a substrate [9], the ignition time of wood [14] and the diffusion of solute in freezing salt solutions [4]. For a more complete list see [5]. Any improvement in the method can therefore have far reaching implications. For the basic problem of heating a semi-in?nite material occupying x ≥ 0 and initially at a constant temperature the HBIM involves three steps. Firstly, the heat penetration depth, δ(t) is introduced. For x ≥ δ the temperature change from the initial temperature is negligible. Secondly an approximating function, typically a polynomial, is introduced. This describes the temperature for 0 ≤ x ≤ δ(t). Finally, the heat equation is integrated over x a [0,δ] to produce what is termed the heat balance integral. The result of this process is a single ordinary differential equation for δ, which may often be solved analytically. The temperature, u(x,t), is then known through the approximating function which is in terms of δ(t). The greatest drawback with the method is the choice of approximating function. Goodman initially proposed a quadratic, but also brie?y mentioned a cubic. Wood [15] shows that for the quadratic there are six distinct formulations and the one employed by Goodman is generally only the third most accurate. Antic and Hill [1] employ a cubic pro?le in the study of diffusion in grain silos. Myers et al. [12] also

use a cubic to describe the melting of a ?nite thickness block. Mitchell and Myers [8] use a quartic approximation in a study of ablation. To ensure that the approximate solution merges smoothly with the constant far ?eld temperature it is standard to set ux(δ,t) = 0. It may also be shown that this leads to uxx(δ,t) = 0. Further, with the appropriate scaling one can set u(δ,t) = 0 and then the polynomial solutions reduce to xn u = an 1? : δ ?1?

The papers discussed so far involve the assumption that n is an integer. Braga et al. [2] take the next logical step and allow n to be a noninteger. For the ablation problem mentioned above they choose n based on matching the time to melting predicted by the exact solution and the approximate one (in fact their choice ensures that the temperature matches the exact solution at x = 0 for all time). After melting they change n to give better agreement with a numerical solution. Hristov [5] has carried out a detailed study of HBIM problems with a temperature of the form Eq. (1) and also provides a comprehensive review of HBIM studies. He introduces two alternative constraints to determine n based on different ways of matching to an exact solution. The key question then is what is the best choice of n? In most of the works mentioned above an exact or numerical solution is known and n is chosen based on matching to these solutions. This really makes the method redundant; if the exact or numerical solution is already known, why look for an approximate one? In order to measure the accuracy, without knowledge of the exact solution, Langford [7] introduced a function, E?t ? =

∫

δ @ u @ 2 u!2 ? dx 0: @ t @ x2 0

?2?

☆ Communicated by W.J. Minkowycz. E-mail address: tim.myers@ymail.com. 0735-1933/$ – see front matter ? 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.icheatmasstransfer.2008.10.013

If u is an exact solution then obviously E = 0. Approximate solutions will have E N 0. Taking the square of ut ? uxx prevents the cancelling of errors of opposite sign and also magni?es the importance of regions

144

T.G. Myers / International Communications in Heat and Mass Transfer 36 (2009) 143–147

where u does not closely satisfy the heat equation [7]. Langford focussed on the problem where the temperature at x = 0 is ?xed, taking n = 1, 2, 3 or 4, and found E = ent? 3/2, where en is a function of n alone. In the following work, rather than using Eq. (2) to indicate the error for a chosen n, we will minimise E to determine n. In doing so we prove that E ~ t? 3/2 for the ?xed temperature boundary condition (a result not shown in the short paper of Langford). For a ?xed ?ux condition we ?nd E ~ t? 1/2. For a cooling condition no such relation can be obtained. We illustrate the method on standard thermal problems, for which an exact solution is known, using the HBIM and related Re?ned Integral Method (RIM), see [13]. However, it is important to note that at no time do we use the exact solution to guide our approximate ones. Comparison of our results with the exact solution shows excellent agreement. 2. Application to standard thermal problems We now look at three standard thermal problems, each with an exact solution. In each case we assume a semi-in?nite media, initially at a constant temperature. Heat is then input at the boundary, either by specifying a new temperature there, a constant heat ?ux or a Newton cooling condition. Only for the ?rst problem do we show the non-dimensionalisation but obviously the scaling is similar in all subsequent cases. 2.1. Fixed temperature boundary condition Consider the problem where @u @2u =α 2 @t @x u?0; t ? = us ujxY∞ Yui u?x; 0? = ui ; ?3?

Substituting these into the expression for the error, Eq. (2), and evaluating the integral we ?nd h i n2 ?2n + 1? ?2n?1??n?1?2 ?δδt ?2n?3? + δ2 δ2 t n?2n?3? E?t ? = : ?8? 3 2 δ ?2n?3??4n ?1? This is the function we must minimise to determine n. At this stage we need to determine the heat penetration depth, δ. The HBIM requires integrating the heat equation over the region x a [0, δ]. Since u(δ, t) = ux(δ, t) = 0 this leads to d dt

∫

δ

0

udx = ?

@u @x

j

x=0

:

?9?

Substituting for u from Eq. (1) gives a single ordinary differential equation for δ with solution δ= p??????????????????????? 2n?n + 1?t ; ?10?

where δ(0) = 0. Note, δδt = n(n + 1) is independent of t and so the numerator of Eq. (8) is also independent of t and therefore, as pointed out by Langford [7], E = ent? 3/2 where en is a constant. Substituting for δ into Eq. (8) we ?nd ? E?t ? = 6n2 + 2n4 + 2n?7n3 ?1 ?p???????????????????? 2n?n + 1? t ?3=2 = en t ?3=2 ; ?11?

4?n + 1?2 ?2n?3??4n2 ?1?

which de?nes en. This formula shows singularities at n = ? 1, ±1/2 and n = 3/2. However for u to show a smooth transition to the constant solution at x = δ we assume n ≥ 2. The RIM formulation involves integrating the heat equation twice with respect to x

where u denotes the temperature, α the thermal diffusivity and us and ui are the temperature at x = 0 and initial temperature. We nondimensionalise so that ? = u?ui u us ?ui ?=x t ?= t x L τ ?4?

∫ ∫

0

δ

x

0

@u d dx = @t

∫

δ

0

@u @u ? @x @x

j

dx: x=0

?12?

and immediately drop the hats. Choosing τ = L2/α gives @u @2u = @ t @ x2 u?0; t ? = 1 ujxY∞ Y0 u?x; 0? = 0: ?5?

On the right hand side we see that ux(0, t) is a function of time and so may be treated as a constant in the integral. The other term integrates immediately. The left hand side may be integrated by parts, so δ

∫

δ

0

@u d? @t

∫ x @t dx = uj

δ 0

@u

x=δ

?u

j

x=0

?δ

@u @x

j

x=0

:

?13?

p???? ? This has the well-known exact solution u = erfc x=2 t ?. The HBIM and RIM involve choosing an approximating function that holds over a ?nite region δ(t), known as the heat penetration depth. For this example the heat penetration depth is de?ned by u(δ,t) =0 and ux(δ,t) =0. A further condition uxx(δ,t)=0 can be obtained by a smoothness argument or taking the total derivative (with respect to time) of u(δ,t)=0 [7,8,15]. If the approximating function is a polynomial of the form x x 2 x 3 u = a0 + a1 1? + a2 1? + a3 1? ; δ δ δ ?6?

Noting ut = uxx the ?rst term on the left hand side may be integrated and we arrive at the main equation for the RIM d dt

∫

δ

0

xu dx = u?0; t ?:

?14?

Substituting for u we ?nd an equation for the heat penetration depth, with solution δ= p??????????????????????????????? ?n + 1??n + 2?t : ?15?

then for a quadratic approximation a3 = 0 and the conditions u(δ,t) = 0 and ux(δ,t) = 0 determine a0 = a1 = 0 while the condition at x = 0 determines a2 = 1. For the cubic the smoothness condition also gives a2 = 0 and then we ?nd a3 = 1. In either case we are left with a function of the form of Eq. (1) where an = 1 (where n = 2 or 3 for the quadratic and cubic cases). With the temperature of Eq. (1) we ?nd @ u nxδt xn?1 = 2 1? @t δ δ @ 2 u n?n?1? xn?2 = 1? : δ @ x2 δ2 ?7?

This gives ? ?p???????????????????????????? n 10n5 ?39n4 + 34n3 + 27n2 ?20n?12 ?n + 1??n + 2? ?3=2 E?t ? = t = en t ?3=2 : 4?2n?3??4n2 ?1??n + 1?2 ?n + 2?2 ?16? Once δ is known we are left with the problem of determining n. This is achieved by minimising en. On Fig. 1 we show plots of en for the HBIM and RIM formulations for n a [0, 5]. For a sensible temperature pro?le we require n ≥ 2 and so ?nd a single minimum for each formulation. For the HBIM we ?nd a minimum of en ≈ 0.0169 for n ≈ 2.233. The minimum for the RIM formulation occurs for very similar n ≈ 2.219 with a marginally smaller error, en ≈ 0.0167. No other

T.G. Myers / International Communications in Heat and Mass Transfer 36 (2009) 143–147

145

3. Constant ?ux boundary condition We now change the problem by replacing the boundary condition at x = 0 with the constant ?ux condition @u ?0; t ? = ?1: @x This problem has the exact solution u=2

Fig. 1. Plot of en(n) for n a [0,5] for the constant temperature condition and the HBIM formulation (dashed line) and the RIM formulation (solid line).

?17?

r??? t ?x2 =?4t ? x e ?x erfc p?? π 2 t:

?18?

The approximate temperature now has the form u= δ xn 1? n δ ?19?

which leads to h i n?2n + 1? n?2n?1??n?1?2 + δt δ?2n?3??3n?2? + δ2 δ2 t ?5n?2??2n?3? δn2 ?4n2 ?1??2n?3? ?20? The HBIM formulation is still described by Eq. (9) but with the right hand side simpli?ed with ux(0,t) = ? 1. Substituting for u from Eq. (19) into Eq. (9) leads to

Fig. 2. Plot of u(x,t) at t = 1. The exact solution is shown as a solid line, the RIM solution with n = 2.219 as a dashed line and the HBIM solution with n = 2.233 as a dash–dot line.

E?t ? =

:

δ=

p???????????????????? n?n + 1?t : The RIM formulation is given by Eq. (14) and leads to

?21?

minima occur for larger n and in fact these curves asymptote to 0.088 for HBIM and 0.312 for the RIM as n→∞. Since en is independent of t these values of n will always provide the best approximation. On Fig. 2 we show a comparison of the exact solution against the approximate temperature from the HBIM and RIM formulations at t = 1. This means that for the HBIM we set n = 2.233 and use this value to calculate δ from Eq. (10), then u is given by Eq. (1). For the RIM we choose n = 2.219 and use Eq. (15) to calculate δ. The approximate curves show good agreement with the exact solution. This good agreement continues for all time (this has been con?rmed by testing the solutions for a wide range of times t a [10? 4,10]). For this problem we have found that n is close to 2, we therefore do not plot the standard Goodman solution with n = 2 since this is quite close to the curves already shown and confuses the ?gure. However, we will show it in the following section.

δ=

r?????????????????????????????????? 2?n + 1??n + 2?t : 3

?22?

In both cases δδt is independent of time and we can write an equation of the form E = ent? 1/2. A comparison of en for the HBIM and RIM formulations is shown in Fig. 3 for n a [1.7,5] (this range avoids the singularities and so highlights the required value of n). The minimum error and corresponding n values for HBIM and RIM are (en, n) = (0.0024, 3.584), (0.0029, 3.822) respectively. Note, these errors are an order of magnitude less than those in the previous example, we therefore expect even better agreement between the temperature pro?les. In Fig. 4 we compare the exact (solid line) and HBIM (dotted line) temperature predictions at time t = 1. The HBIM solution is almost indistinguishable from the exact solution. We do not show the RIM solution since it is virtually identical to the exact and HBIM solutions.

Fig. 3. Plot of en(n) for n a [1.7, 5] for a constant ?ux condition with the HBIM formulation (dash–dot line) and the RIM formulation (solid line).

Fig. 4. Comparison of exact (solid line), HBIM with n = 3.584 (dashed line) and HBIM with n = 2 (dash–dot) temperatures for constant ?ux condition at time t = 1.

146

T.G. Myers / International Communications in Heat and Mass Transfer 36 (2009) 143–147

Fig. 5. Plot of E(t) at t = 0.5, and t = 1 for n a [1.75, 5] for a cooling condition. The HBIM formulation is shown as a dashed line, the RIM formulation as a solid line.

Fig. 6. Plot of u(x, t) at t = 1. The exact solution is shown as a solid line, the RIM solution with n = 3.82, δ = 4.675 (dashed), the HBIM solution with n = 3.58, δ = 4.59 (dash–dot) and the HBIM solution (dotted) with n = 2, δ = 2.79.

So we conclude that either technique can be applied to the constant ?ux problem to obtain excellent agreement with the exact solution provided the choice n = 3.584 or n = 3.822 is made for the respective methods. The third curve is the standard HBIM solution of Goodman, with n = 2; it clearly gives poor agreement. An HBIM analysis with constant ?ux was carried by Braga et al. [2] whilst looking at the heating up phase before ablation starts. Their exponent, n = π / (4 ? π) ≈ 3.66 was chosen to match the time predicted by the exact solution when the surface x = 0 reaches the ablation temperature (in fact their solution matches the exact temperature at x = 0 for all time). Obviously this is very close to our value which was chosen without knowledge of the exact solution. Our HBIM and RIM formulations show 0.45% and 0.62% errors in the prediction of melting time. However, our general temperature pro?le will always be more accurate than when their value of n is used. Mitchell and Myers [8] choose n = 4 and use the HBIM, from the plot of en it is clear that this will show an even greater error. 4. Cooling condition Finally we replace the boundary condition at x = 0 with a Newton cooling condition @u ?0; t ? = u?0; t ??1: @x The exact solution is x x + 2t p?? : u = erfc p?? ?ex + t erfc 2 t 2 t The approximating function that satis?es Eq. (23) is u= δ xn 1? : n+δ δ The HBIM formulation gives " # d δ2 n?n + 1? = dt ?n + δ? n+δ with the implicit solution # 1 δ2 n+δ + nδ?n2 ln : t= n?n + 1? 2 n " ?25? ?24? ?23?

The RIM formulation gives # d δ3 δ ?1 + n??2 + n? = dt n + δ n+δ with the implicit solution t= ! 1 n+δ δ2 + nδ?n2 ln : ?n + 1??n + 2? n ?29? " ?28?

Obviously we have now lost the simple square root behaviour of δ. The expression for E is also rather cumbersome and will not be quoted. Instead for this case we will solve the problem numerically. This involves specifying a time t. For a given n either Eq. (27) or Eq. (29) is solved to determine δ, then δt may be calculated using the value of δ from Eq. (26) or Eq. (28). The temperature derivatives required for E are i @u δt xn?1 h x xn = ?n + δ? 1? n 1? + 2 @ t ?n + δ? δ δ δ @ 2 u n?n?1? xn?2 1? = : δ @ x2 δ?n + δ? ?30?

?31?

The integral E may now be evaluated numerically. This process is repeated over the speci?ed range of n. Note, for this situation we cannot write E = ent α (in fact from the analytical expression for E it is apparent that this form is not appropriate). The result is therefore qualitatively different to the previous examples and it is not certain that the position of the minimum of E is ?xed in time. For this reason, in Fig. 5a, b we plot E(t) at times t = 0.1, 0.5. In Fig. 5a the minima occur at n = 3.138, 3.313 for the HBIM and RIM. In Fig. 5b they occur at n = 2.845, 2.966. If we let t→0 we ?nd n→3.58, 3.82 (which are the values given by the constant ?ux analysis). The minimum error rapidly decreases with time. So n is indeed a function of time in this case. In Fig. 6 we compare temperatures at time t = 1. Since n is a decreasing function of t we show this result to indicate the level of accuracy achieved by using a constant n. We choose the n value for t ≈ 0, since this is when the greatest error will occur. Since the error

?26?

Table 1 Optimal values of n for different thermal problems Boundary condition HBIM n 2.235 3.584 3.58 RIM n 2.218 3.82 3.82 Typical error, HBIM en = 0.0169 en = 0.0024 E(t) ≈ 0.005 Typical error, RIM en = 0.0167 en = 0.0029 E(t) ≈ 0.003

?27?

Fixed temperature Constant ?ux Cooling condition

T.G. Myers / International Communications in Heat and Mass Transfer 36 (2009) 143–147

147

rapidly decreases with time minimising the error at t = 0 will lead to the smallest overall error. Consequently we take n = 3.58, 3.82 for the HBIM and RIM. We then calculate δ from Eq. (27) or Eq. (29) using these n values at the current time and ?nd E(1) = 0.005, 0.0035. Four solutions are shown, the exact solution the RIM with n = 3.82 and HBIM with n = 3.58, 2. Although not as accurate as in the previous example the ?rst two approximate solutions are still very good. The standard solution, with n = 2, provides a good approximation for small x but is very poor for larger values. 5. Discussion The aim of this work was to determine the optimal value of the exponent n for certain thermal problems. Where optimal has been de?ned in the sense of the Langford criteria of minimising the square of the heat equation. No doubt other criteria could be used, however Langford's method has the advantage over standard error measures in that it does not require knowledge of an exact solution. Further, since it involves the square of the solution, errors are magni?ed and so inaccurate solutions are punished more severely than accurate ones. We have demonstrated the method with three standard thermal problems and in all cases shown excellent agreement with exact solutions. We used problems with an exact solution simply for comparison. The method does not require any knowledge of these solutions. Although we have worked throughout with non-dimensional quantities the results hold for the equivalent dimensional systems. In Table 1 we summarise the cases and the appropriate value of n. This is the exponent of the polynomial that will give the best approximation (with respect to the error of Langford) to the solution of the heat equation with the speci?ed boundary condition. The error when using HBIM or RIM is similar in all cases. We have not carried out an exhaustive study of applications of the integral methods, since there are too many. However, we have illustrated the application of the new method to typical problems and

hope this is suf?cient to guide the user in other applications. In future we will extend the work to Stefan problems, which are, of course, where these methods are of the most practical use. References

[1] A. Antic, J.M. Hill, The double-diffusivity heat transfer model for grain stores incorporating microwave heating, Appl. Math. Model. 27 (2003) 629–647. [2] W.F. Braga, M.B.H. Mantelli, J.L.F. Azevedo, Approximate analytical solution for onedimensional ?nite ablation problem with constant time heat ?ux, AIAA Thermophys. Conference, June 2004. [3] T.R. Goodman, The heat-balance integral and its application to problems involving a change of phase, Trans. ASME 80 (1958) 335–342. [4] B.W. Grange, R. Viskanta, W.H. Stevenson, Diffusion of heat and solute during freezing of salt solutions, Int. J. Heat & Mass Trans. 19 (4) (1976) 373–384. [5] J. Hristov, The heat-balance integral method by a parabolic pro?le with unspeci?ed exponent: analysis and exercises, to appear in Thermal Sci. 13(2) 2009. [6] S. Kutluay, A.S. Wood, A. Esen, A heat balance integral solution of the thermistor problem with a modi?ed electrical conductivity, Appl. Math. Model. 30 (2006) 386–394. [7] D. Langford, The heat balance integral method, Int. J. Heat Mass Trans. 16 (1973) 2424–2428. [8] S.L. Mitchell, T.G. Myers, Heat balance integral method for one-dimensional ?nite ablation, AIAA, J. Thermophys. & Heat Trans. 22 (3) (2008) 508–514. [9] S.L. Mitchell, T.G. Myers, Approximate solution methods for one-dimensional solidi?cation from an incoming ?uid, App. Math. Comp. 202 (1) (2008) 311–326. [10] T.G. Myers and S.L. Mitchell, Application of the heat balance and re?ned integral methods to the Korteweg-de Vries equation, to appear in Thermal Sci. 13(2) 2009. [11] T.G. Myers, S.L. Mitchell, and G. Muchatibaya, Unsteady contact melting of a rectangular cross-section material on a ?at plate, Phys. Fluids 20 (2008) 103101. doi:10.1063/1.2990751. [12] T.G. Myers, S.L. Mitchell, G. Muchatibaya, M.Y. Myers, A cubic heat balance integral method for one-dimensional melting of a ?nite thickness layer, Int. J. Heat Mass Trans. 50 (2007) 5305–5317. [13] N. Sadoun, E.-K. Si-Ahmed, P. Colinet, On the re?ned integral method for the onephase Stefan problem with time-dependent boundary conditions, App. Math. Model. 30 (2006) 531–544. [14] M.J. Spearpoint, J.G. Quintiere, Predicting the piloted ignition of wood in the cone calorimeter using an integral model — effect of species, grain orientation and heat ?ux, Fire Safety J. 36 (4) (2001) 391–415. [15] A.S. Wood, A new look at the heat balance integral method, App. Math. Model. 25 (2001) 815–824.

相关文章:

- 热分配计法热计量改造_图文
- 5.12.4.2 系统原理概述 5.12.4.2.1 散热器热分配
*计法*是利用散热器热分配计...静态*平衡*阀安装在回水管道上,热量*计的积分*仪及远程模块需 安装在防护箱内,防护...

- 平衡计分卡7556239637
- 在线互动式文档分享平台,在这里,您可以和千万网友分享自己手中的文档,全文阅读其他用户的文档,同时,也可以利用分享文档获取的
*积分*下载文档

- 平衡计分卡
- 在线互动式文档分享平台,在这里,您可以和千万网友分享自己手中的文档,全文阅读其他用户的文档,同时,也可以利用分享文档获取的
*积分*下载文档

- 基于平衡计分卡法(BSC)的绩效考核设计
- 7 (一) 、应用
*平衡计分*卡法简历绩效目标指标体系???7 (二) 、应用*平衡计分*卡法于机关内部绩效管理和公务员队伍建设???8 五、改进和完善柳州市旅游局绩效考核...

- 辐射换热的计算
- 三、角系数的计算方法 直接
*积分法*求解角系数的方法 代数*积分法*几何分析法 1...相对性是在两物体处于*热平衡*时,净辐射换热量 为零的条件下导得的;完整性反映...

- 平衡计分卡在苏泊尔有限公司的应用研究
*平衡计分*卡在苏泊尔有限公司的应用研究_经管营销_专业资料。*平衡计分*卡在苏泊尔有限...再 加上管理者为实现个人目标利益最大化,常常会设法寻找时间短、见效快的项目...

- 九 溶解热的测定
- 在线互动式文档分享平台,在这里,您可以和千万网友分享自己手中的文档,全文阅读其他用户的文档,同时,也可以利用分享文档获取的
*积分*下载文档

- 双汇集团平衡计分卡应用研究 正文
- 双汇集团
*平衡计分*卡应用研究 正文_经管营销_专业资料。沈阳理工大学学士学位论文 ...用双汇集团董事长万隆的话说:“没有成功的IT系统,双汇做不了这么大, 也没法...

- 新课标2018版高考物理一轮复习第十二章热学专题十三热...
- 重力加速度大小取 g=10 m/s .求: 2 5 l (1)在大活塞与大圆筒底部接触前的瞬间,汽缸内封闭气体的温度; (2)缸内封闭的气体与缸外大气达到
*热平衡*时,缸内...

- 7.4 温度和温标
- 在线互动式文档分享平台,在这里,您可以和千万网友分享自己手中的文档,全文阅读其他用户的文档,同时,也可以利用分享文档获取的
*积分*下载文档

更多相关标签: