daily pastebin goal
42%
SHARE
TWEET

Untitled

a guest Mar 21st, 2019 57 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. Needs["NDSolve`FEM`"];
  2. Needs["DifferentialEquations`NDSolveProblems`"];
  3. Needs["DifferentialEquations`NDSolveUtilities`"];
  4.      
  5. Lr = 2*10^-2; (*dimension of computational domain in r-direction*)
  6. Lz = 10^-2;   (*dimension of computational domain in z-direction*)
  7. mesh = ToElementMesh[FullRegion[2], {{0, Lr}, {0, Lz}},MaxCellMeasure -> {"Length" -> Lr/50}, "MeshOrder" -> 1]
  8. mesh["Wireframe"]
  9.      
  10. lambda = 22;         (*heat conductivity*)
  11. density = 7200;      (*density*)
  12. Cs = 700;            (*specific heat capacity of solid*)
  13. Cl = 780;            (*specific heat capacity of liquid*)      
  14. LatHeat = 272*10^3;  (*latent heat of fusion*)
  15. Tliq = 1812;         (*melting temperature*)
  16. MeltRange = 100;     (*melting range*)
  17. To = 300;            (*initial temperature*)      
  18. SPow = 1000;         (*source power*)
  19. R = Lr/4;            (*radius of heat source spot*)
  20. a = Log[100]/R^2;            
  21. qo = (SPow*a)/Pi;
  22. q[r_] := qo*Exp[-r^2*a]; (*heat flux distribution*)        
  23. tau = 10^-3;         (*time step size*)
  24. ProcDur = 0.2;       (*process duration*)
  25.      
  26. Heviside[x_, delta_] := Module[{res},                                                                
  27.                                res = Piecewise[
  28.  
  29.                                                {                                                                  
  30.                                                 {0, Abs[x] < -delta},                                                                      
  31.                                                 {0.5*(1 + x/delta +  1/Pi*Sin[(Pi*x)/delta]), Abs[x] <= delta},                                                                                        
  32.                                                 {1, x > delta}                                                                      
  33.                                                }
  34.  
  35.                                               ];
  36.                                              res
  37.                               ]
  38.      
  39. HevisideDeriv[x_, delta_] := Module[{res},                                                                      
  40.                                     res = Piecewise[                                                                      
  41.                                                    {
  42.  
  43.                                                     {0, Abs[x] > delta},
  44.  
  45.                                                     {1/(2*delta)*(1 + Cos[(Pi*x)/delta]), Abs[x] <= delta}                                                                      
  46.                                                    }                                                                      
  47.                                                    ];                                                                      
  48.                                     res                                                                      
  49.                                   ]
  50.      
  51. EffectHeatCapac[tempr_] := Module[{phase},                                                                      
  52.                                   phase = Heviside[tempr - Tliq, MeltRange/2];
  53.                                   Cs*(1 - phase) + Cl*phase +LatHeat*HevisideDeriv[tempr - Tliq, 0.5*MeltRange]                                                                      
  54.                                  ]
  55.      
  56. ts = AbsoluteTime[];
  57. xlast = Table[{To}, {methodData["DegreesOfFreedom"]}];
  58. TemprField = ElementMeshInterpolation[{mesh}, xlast];
  59. NumTimeStep = Floor[ProcDur/tau];
  60. For[i = 1, i <= NumTimeStep, i++,
  61.  
  62.    (*
  63.     (*Setting of PDE coefficients for linear problem*)
  64.       pdeCoefficients=InitializePDECoefficients[vd,sd,"ConvectionCoefficients"->     {{{{-lambda/r, 0}}}},
  65. "DiffusionCoefficients" -> {{-lambda*IdentityMatrix[2]}},
  66. "DampingCoefficients" -> {{Cs*density}}];    
  67.    *)
  68.  
  69. (*Setting of PDE coefficients for nonlinear problem*)
  70.  
  71.  pdeCoefficients =
  72.  InitializePDECoefficients[vd, sd,
  73.  "ConvectionCoefficients" -> {{   {{-(lambda/r), 0}}  }},
  74.  "DiffusionCoefficients" -> {{-lambda*IdentityMatrix[2]}},
  75.  "DampingCoefficients" -> {{EffectHeatCapac[TemprField[r, z]]*
  76.  density}}];
  77.  
  78.  discretePDE = DiscretizePDE[pdeCoefficients, methodData, sd];
  79.  {load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
  80.  DeployBoundaryConditions[{load, stiffness, damping},
  81.  discreteBCs];
  82.  
  83.  A = damping/tau + stiffness;
  84.  b = load + damping.xlast/tau;
  85.  
  86.  x = LinearSolve[A,b,Method -> {"Krylov", Method -> "BiCGSTAB",
  87.  "Preconditioner" -> "ILU0"}];
  88.  TemprField = ElementMeshInterpolation[{mesh}, x];
  89.  xlast = x;            
  90.  ]
  91. te = AbsoluteTime[];
  92. te - ts
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top