{ ELASTICITY.PDE This example shows the application of FlexPDE to a complex problem in thermo-elasticity. The equations of thermal diffusion and plane strain are solved simultaneously to give the thermally-induced stress and deformation in a laser application. A rod amplifier of square cross-section is imbedded in a large copper heat-sink. The rod is surrounded by a thin layer of compliant metal. Pump light is focussed on the exposed side of the rod. We wish to calculate the effect of the thermal load on the laser rod. The equations of Stress/Strain arise from the balance of forces in a material medium, expressed as dx(Sx) + dy(Txy) + Fx = 0 dx(Txy) + dy(Sy) + Fy = 0 where Sx and Sy are the stresses in the x- and y- directions, Txy is the shear stress, and Fx and Fy are the body forces in the x- and y- directions. The deformation of the material is described by the displacements, U and V, from which the strains are defined as ex = dx(U) ey = dy(V) gxy = dy(U) + dx(V). The eight quantities U,V,ex,ey,gxy,Sx,Sy and Txy are related through the constitutive relations of the material. In general, Sx = C11*ex + C12*ey + C13*gxy - b*Temp Sy = C12*ex + C22*ey + C23*gxy - b*Temp Txy = C13*ex + C23*ey + C33*gxy In orthotropic solids, we may take C13 = C23 = 0. Combining all these relations, we get the displacement equations: dx[C11*dx(U)+C12*dy(V)] + dy[C33*(dy(U)+dx(V))] + Fx = dx(b*Temp) dy[C12*dx(U)+C22*dy(V)] + dx[C33*(dy(U)+dx(V))] + Fy = dy(b*Temp) The "Plane-Strain" approximation is appropriate for the cross-section of a cylinder which is long in the Z-direction, and in which there is no Z-strain. The cylinder is loaded by surface tractions and body forces applied along the length of cylinder, and which are independent of Z. In this case, we may write C11 = G*(1-nu) C12 = G*nu b = G*alpha*(1+nu) C22 = G*(1-nu) C33 = G*(1-2*nu)/2 where G = E/[(1+nu)*(1-2*nu)] E is Young's Modulus nu is Poisson's Ratio and alpha is the thermal expansion coefficient. The displacement form of the stress equations (for uniform temperature and no body forces) is then (dividing out G): dx[(1-nu)*dx(U)+nu*dy(V)] + 0.5*(1-2*nu)*dy[dy(U)+dx(V)] = alpha*(1+nu)*dx(Temp) dy[nu*dx(U)+(1-nu)*dy(V)] + 0.5*(1-2*nu)*dx[dy(U)+dx(V)] = alpha*(1+nu)*dy(Temp) In order to quantify the "natural" (or "load") boundary condition mechanism, consider the stress equations in their original form: dx(Sx) + dy(Txy) = 0 dx(Txy) + dy(Sy) = 0 These can be written as div(P) = 0 div(Q) = 0 where P = [Sx,Txy] and Q = [Txy,Sy] The natural (or "load") boundary condition for the U-equation defines the outward surface-normal component of P, while the natural boundary condition for the V-equation defines the surface-normal component of Q. Thus, the natural boundary conditions for the U- and V- equations together define the surface load vector. On a free boundary, both of these vectors are zero, so a free boundary is simply specified by load(U) = 0 load(V) = 0. -- Submitted by Steve Sutton, Lawrence Livermore National Laboratory } TITLE "Thermo-Elastic Stress" COORDINATES cartesian2 { coordinate system, 1D,2D,3D, etc } VARIABLES { system variables } Tp Up Vp ! Wp { para problema 3D } SELECT { method controls } errlim = 1.0e-3 painted = on modes=2 { calcula 2 autovalores para o caso 2D; MUDAR PARA 3 NO CASO 3D ! ! ! } DEFINITIONS { parameter definitions } tf = 720 { tempo final de cálculo } g = 9.81 { gravidade, em m/s^2 } Text = 28 { temperatura externa; pode ser especificada uma função senoidal ou outra qualquer do tempo } { PROPRIEDADES DAS CAMADAS } T_lanc { declaração da temperatura de lançamento da camada - o valor será informado na respectiva REGION } tz { declaração do tempo de lançamento da camada - o valor será informado na respectiva REGION } { PROPRIEDADES DO CONTORNO DA CAMADAS } duracao_hc { declaração do tempo que o coeficiente de convecção fica com certo valor- o valor será informado no respectivo contorno da geometria da REGION } hc { declaração do coeficiente de convecção - o valor será informado no respectivo contorno da geometria da REGION } { AVANÇO DA HIDRATAÇÃO } GH { grau de hidratação em função do tempo e do tempo de início da hidratação } { PROPRIEDADES TÉRMICAS } k { declaração da condutividade térmica - o valor será informado na respectiva REGION } rho { declaração da massa específica - o valor será informado na respectiva REGION } c { declaração do calor específico - o valor será informado na respectiva REGION } Tad { declaração da temperatura adiabática - a função do tempo ou do grau de hidratação será informada na respectiva REGION } alpha { declaração do coeficiente de expansão térmica linear - será informado na respectiva REGION } { PROPRIEDADES ELÁSTICAS E VISCOELÁSTICAS (FLUÊNCIA) } E { declaração do módulo de elasticidade - a função do tempo ou do grau de hidratação será informada na respectiva REGION } nu { declaração do coeficiente de Poisson - o valor será informado na respectiva REGION } phi { declaração do coeficiente de fluência - a função do tempo ou do grau de hidratação será informada na respectiva REGION } { DEFINIÇÃO DE RELAÇÕES CONSTITUTIVAS PARA MATERIAL ISOTRÓPICO } lamb = E*nu/((1.+nu)*(1.-2.*nu)) { primeiro parâmetro de Lamé } mi = E/(2*(1+nu)) { segundo parâmetro de Lamé } { expressão do tensor de Hooke para o caso isotrópico } C11 = lamb+2*mi C22 = lamb+2*mi ! C33 = lamb+2*mi ! C23 = lamb ! C31 = lamb C12 = lamb ! C44 = mi ! C55 = mi C66 = mi { relações entre deformações e deslocamentos, para o caso ortotrópico, que inclui o caso isotrópico } ex = dx(Up) ey = dy(Vp) ! ez = dz(Wp) gxy = dy(Up) + dx(Vp) ! gyz = dz(Vp) + dy(Wp) ! gzx = dx(Wp) + dz(Up) { Lei de Hooke } Sxx = C11*ex + C12*ey {+ C13*ez} - (3*lamb + 2*mi)*alpha*(Tp-T_lanc) Syy = C12*ex + C22*ey {+ C23*ez} - (3*lamb + 2*mi)*alpha*(Tp-T_lanc) ! Szz = C13*ex + C23*ey + C33*ez - (3*lamb + 2*mi)*alpha*(Tp-T_lanc) ! Tyz = C44*gyz ! Tzx = C55*gzx Txy = C66*gxy S = TENSOR ( (Sxx,Txy{,Tzx}) , (Txy,Syy{,Tzx}) {, (Tzx,Tyz,Szz)}) u = vector(Up,Vp{,Wp}) { expressão do vetor dos deslocamentos } qx = dx(Tp) qy = dy(Tp) ! qz = dz(Tp) q = vector(qx,qy{,qz}) { Critério de Von Mises - OBS: válido apenas para 2D } p_ang= 0.5* arctan( 2*Txy/(Sxx-Syy) ) { Radians } Sxp0= 0.5*(Sxx+Syy)+ 0.5*(Sxx-Syy)*cos(2*p_ang)+ Txy*sin(2*p_ang) Syp0= 0.5*(Sxx+Syy)+ 0.5*(Sxx-Syy)*cos(2*(p_ang+pi/2)) + Txy*sin(2*(p_ang+pi/2)) { Test for highest value } p_angl= if Sxp0>Syp0 then p_ang else p_ang+ pi/2 Sxp= if Sxp0>Syp0 then Sxp0 else Syp0 Syp= if Sxp0>Syp0 then Syp0 else Sxp0 p_angle= p_angl * 180/pi { Degrees } mises{Sv=sqrt(3*J2)}= sqrt( 0.5*( (Sxx-Syy)^2 + (Syy{-Szz})^2 + ({Szz-}Sxx)^2)+ 3*(Txy^2{+Tyz^2+Tzx^2})) INITIAL VALUES Tp = T_lanc { a temperatura inicial será a temperatura de lançamento de cada camada } EQUATIONS { PDE's, one for each variable } { primeira lei de Newton - equação do movimento: div(S) + F = 0 } Up: dx(Sxx) + dy(Txy){ + dz(Tzx)} = 0 { the U-displacement equation } Vp: dx(Txy) + dy(Syy){ + dz(Tyz)} - rho*g = 0 { the V-displacement equation; CUIDADO AO CONSIDERAR O PROLEMA EM 3D E A VERTICAL FOR O EIXO Z ! ! ! } ! Wp: dx(Tzx) + dy(Tyz) + dz(Szz) = 0 { the W-displacement equation } Tp: div(k*grad(Tp)) + rho*c*dt(Tad)= rho*c*dt(Tp) { equação do calor } ! S: div(grad(S)) + lambda*S = 0 {CONSTRAINTS { Integral constraints } integral(Up)=0 { eliminate translations } integral(Vp)=0 ! integral(Wp)=0 integral(dx(Vp)-dy(Up)) = 0 { eliminate rotations } ! integral(dy(Wp) - dz(Vp)) = 0 ! integral(dz(Up) - dx(Wp)) = 0 } BOUNDARIES { The domain definition } REGION 1 { por convenção, referir-se-á à fundação mais o contorno do domínio em estudo } { as variáveis anteriormente declaradas, especificadas aqui para a fundação, na verdade valem para o domínio inteiro, uma vez que aqui este é definido } { nas demais regiões, as variáveis aqui definidas para todo o domínio, assumem os valores especificados em cada região } { propriedades da fundação} T_lanc = 28 { especificar a temperatura inicial da fundação } tz = 0 { a fundação tem tempo de lançamento igual a zero } { PROPRIEDADES DO CONTORNO DA CAMADAS } duracao_hc = 72 { especificação do tempo que o coeficiente de convecção fica com certo valor } hc = if t - tz <= duracao_hc then 300 else 12 { em geral, não se considera a convecção no contorno da fundação } { mas pode ser especificada uma troca de calor com o ambiente e, portanto, hc diferente de 0 } { AVANÇO DA HIDRATAÇÃO } GH = 0 { grau de hidratação em função do tempo e do tempo de início da hidratação } { PROPRIEDADES TÉRMICAS } k = 2.6 { especificação da condutividade térmica da fundação } rho = 2400 { especificação da massa específica da fundação } c = 0.26 { especificação do calor específico da fundação } Tad = 0 { a temperatura adiabática é sempre igual a zero } alpha = 1.e-6 { especificação do coeficiente de expansão térmica linear da fundação } { PROPRIEDADES ELÁSTICAS E VISCOELÁSTICAS (FLUÊNCIA) } E = 30e9 { especificação do módulo de elasticidade da fundação } nu = 0.20 { especificação do coeficiente de Poisson da fundação } phi = 0 { especificação do coeficiente de fluência da fundação } { Walk the domain boundary } start 'fora'(0,15) value(Up)=0 value(Vp)=0 {value(Wp)=0} line to (56,15) line to (56,44) line to (44,44) natural(Tp) = hc*(Tp-Text) line to (0,44) natural(Tp) = hc*(Tp-Text) line to (0,27) line to close REGION 2 "bloco de concreto" { a região 2 e as subsequentes referir-se-ão às camadas de concretagem } { uma camada deve ser apenas de um material; se uma camada tiver dois materiais, considerar duas camadas diferentes com o mesmo tempo de lançamento } { PROPRIEDADES DA CAMADA } T_lanc = 28 { especificação da temperatura inicial ou de lançamento da camada } tz = 0 { especificação do tempo de lançamento da camada } { PROPRIEDADES DO CONTORNO DA CAMADAS } duracao_hc = 0 { especificação do tempo que o coeficiente de convecção fica com certo valor } hc = if t - tz <= duracao_hc then 300 else 12 { especificação do coeficiente de convecção da camada } { AVANÇO DA HIDRATAÇÃO } GH = if t <= tz then 0 else ( t - tz) { grau de hidratação em função do tempo e do tempo de início da hidratação } { PROPRIEDADES TÉRMICAS } k = if t <= tz then 0 else 2.6 { especificação da condutividade térmica da camada } rho = if t <= tz then 0 else 2400 { especificação da massa específica da camada } c = if t <= tz then 0 else 0.26 { especificação do calor específico da fundação } Tad = if t <= tz then 0 else GH^2.512/(0.08605+0.03248*GH^2.512) { especificação da temperatura adiabática da camada } alpha = 1.e-6 { especificação do coeficiente de expansão térmica linear da fundação } { PROPRIEDADES ELÁSTICAS E VISCOELÁSTICAS (FLUÊNCIA) } E = if t <= tz then 0 else 25e9 { função( GH ) } { especificação do módulo de elasticidade da camada } nu = 0.20 { especificação do coeficiente de Poisson da camada } phi = if t <= tz then 0 else 0 { função( GH ) } { especificação do coeficiente de fluência da camada } { Walk the domain boundary } start(0,27) line to (44,27) line to (44,44) line to (0,44) line to close REGION 4 "concreto em torno do tubo" { For each material region } { PROPRIEDADES DA CAMADA } T_lanc = 28 { especificação da temperatura inicial ou de lançamento da camada } tz = 0 { especificação do tempo de lançamento da camada } { PROPRIEDADES DO CONTORNO DA CAMADAS } duracao_hc = 0 { especificação do tempo que o coeficiente de convecção fica com certo valor } hc = if t - tz <= duracao_hc then 300 else 12 { especificação do coeficiente de convecção da camada } { AVANÇO DA HIDRATAÇÃO } GH = if t <= tz then 0 else ( t - tz) { grau de hidratação em função do tempo e do tempo de início da hidratação } { PROPRIEDADES TÉRMICAS } k = if t <= tz then 0 else 2.6 { especificação da condutividade térmica da camada } rho = if t <= tz then 0 else 2400 { especificação da massa específica da camada } c = if t <= tz then 0 else 0.26 { especificação do calor específico da fundação } Tad = if t <= tz then 0 else GH^2.512/(0.08605+0.03248*GH^2.512) { especificação da temperatura adiabática da camada } alpha = 1.e-6 { especificação do coeficiente de expansão térmica linear da fundação } { PROPRIEDADES ELÁSTICAS E VISCOELÁSTICAS (FLUÊNCIA) } E = if t <= tz then 0 else 25e9 { função( GH ) } { especificação do módulo de elasticidade da camada } nu = 0.20 { especificação do coeficiente de Poisson da camada } phi = if t <= tz then 0 else 0 { função( GH ) } { especificação do coeficiente de fluência da camada } start 'esquerda' (15.5,33.5) arc (center=10,33.5) angle=360 start 'direita' (26.5,33.5) arc (center=32,33.5) angle=360 EXCLUDE start (15,33.5) natural(Tp) = hc*(Tp-Text) arc (center=10,33.5) angle=360 start (27,33.5) natural(Tp) = hc*(Tp-Text) arc (center=32,33.5) angle=360 TIME 0 TO tf { período de tempo considerado no cálculo } MONITORS { show progress } PLOTS { save result displays } for t = 0 by tf/100 to tf ! grid(x+Up,y+Vp{,z+Wp}) { mostra a malha utilizada } contour(Tp) as "Resultado do cálculo do campo de temperaturas" contour(Up) contour(Vp) vector(u) norm contour(Sxx) contour(Syy) contour(Txy) elevation(Up, Vp) on 'esquerda' elevation(Sxx) on 'esquerda' elevation(Syy) on 'esquerda' contour( p_angle) vector( cos(p_angl)*Sxp, sin(p_angl)*Sxp) norm notips as ' "Tension" ' contour( Sxp) contour( Syp) contour( Mises) END