# Heat Transfer in Orthotropic Material Plate due to Laser Beam

Find the temperature distribution on a square plate made of orthotropic material. A laser beam hits the middle of the plate and generates a heat flux. This example shows how to

• Specify an orthotropic thermal conductivity that also depends on the temperature using a conductivity matrix and a function handle.

• Define a heat flux pulse of a Gaussian shape using a function handle. The example shows a function modeling a beam with the steady amplitude that is switched off after a particular time, as well as a function modeling the beam with its amplitude gradually decaying over time.

### Thermal Analysis Setup

Create and plot a geometry representing a 40-by-40-by-3 cm plate.

```L = 0.4; % m W = 0.4; % m H = 0.03; % m g = multicuboid(L,W,H); figure pdegplot(g,FaceLabels="on",FaceAlpha=0.3)```

Create a finite element analysis model with the plate geometry for transient thermal analysis.

`model = femodel(AnalysisType="thermalTransient",Geometry=g);`

Specify the initial temperature of the plate, 300 K.

`model.CellIC = cellIC(Temperature=300);`

Generate a mesh.

`model = generateMesh(model);`

Specify the heat flux boundary condition on top of the plate. First, specify a heat flux as a short pulse created by the laser beam by using the `heatFluxSteady` function. See Heat Flux Functions.

`model.FaceLoad(2) = faceLoad(Heat = @heatFluxSteady);`

Plot the heat flux returned by `heatFluxSteady` at the initial time t = 0 s.

```a = linspace(-L/2,L/2); b = linspace(-W/2,W/2); [x,y] = meshgrid(a,b); location.x = x; location.y = y; state.time = 0; flux = heatFluxSteady(location,state); figure contourf(x,y,flux,100,LineStyle="none") colorbar xlabel("x") ylabel("y")```

### Isotropic Thermal Conductivity

First, assume that the plate is made of an isotropic material with the thermal conductivity 50 $W/\left(m\cdot K\right)$.

`K = 50;`

Specify the material properties. The mass density of the plate is 1960 $kg/{m}^{3}$, and the specific heat is 710 $J/\left(kg\cdot K\right)$.

```model.MaterialProperties = ... materialProperties(ThermalConductivity=K, ... MassDensity=1960, ... SpecificHeat=710);```

Solve the thermal problem for the times between 0 and 2 seconds.

```t = 0:0.1:2; Rt_iso = solve(model,t);```

Plot the results using the Visualize PDE Results Live Editor task. First, create a new live script by clicking the New Live Script button in the File section on the Home tab.

On the Live Editor tab, select Task > Visualize PDE Results. This action inserts the task into your script.

To plot the temperature of the plate:

1. In the Select results section of the task, select `Rt_iso` from the drop-down list.

2. In the Specify data parameters section of the task, set Type to Temperature.

3. Select the Animate checkbox to see how the temperature changes over time.

``` ```

### Orthotropic Thermal Conductivity

For orthotropic thermal conductivity, create a vector of thermal conductivity components ${k}_{x}$,${k}_{y}$, and ${k}_{z}$.

```kx = 50; % W/(m*K) ky = 50; % W/(m*K) kz = 500; % W/(m*K) K = [kx;ky;kz];```

Specify material properties of the plate.

```model.MaterialProperties = ... materialProperties(ThermalConductivity=K, ... % W/(m*K) MassDensity=1960, ... % kg/(m^3) SpecificHeat=710); % J/(kg*K)```

Solve the thermal problem.

`Rt_ortho = solve(model,t);`

Plot the results using the Visualize PDE Results Live Editor task.

1. In the Select results section of the task, select `Rt_ortho` from the drop-down list.

2. In the Specify data parameters section of the task, set Type to Temperature.

3. Select the Animate checkbox to see how the temperature changes over time.

``` ```

### Orthotropic Thermal Conductivity Depending on Temperature

For orthotropic thermal conductivity that also depends on temperature, create a vector of thermal conductivity components ${k}_{x}$,${k}_{y}$, and ${k}_{z}$ and use it in a function handle to account for temperature dependency.

```kx = 50; % W/(m*K) ky = 50; % W/(m*K) kz = 500; % W/(m*K) K = @(location,state) [kx;ky;kz] - 0.005*state.u;```

Specify material properties of the plate.

```model.MaterialProperties = ... materialProperties(ThermalConductivity=K, ... % W/(m*K) MassDensity=1960, ... % kg/(m^3) SpecificHeat=710); % J/(kg*K)```

Solve the thermal problem.

`Rt_ortho_t = solve(model,t);`

Plot the results using the Visualize PDE Results Live Editor task.

1. In the Select results section of the task, select `Rt_ortho_t` from the drop-down list.

2. In the Specify data parameters section of the task, set Type to Temperature.

3. Select the Animate checkbox to see how the temperature changes over time.

``` ```

### Thermal Pulse with Amplitude Decay

Now specify a heat flux as a pulse with the amplitude gradually decaying over time by using the `heatFluxTimeDependent` function. See Heat Flux Functions.

`model.FaceLoad(2) = faceLoad(Heat = @heatFluxTimeDependent);`

Visualize the heat flux for all time steps by creating an animation. Keep a fixed vertical scale by scaling all plots to use the maximum amplitude 350000000$W/{m}^{2}$ as a colorbar upper limit.

```state.time = t; flux = heatFluxTimeDependent(location,state); figure for index = 1:length(t) contourf(x,y,flux(:,:,index),100,LineStyle="none") colorbar clim([0 350000000]) xlabel("x") ylabel("y") M(index) = getframe; end```

To play the animation, use the `movie(M)` command.

Solve the thermal problem.

`Rt = solve(model,t);`

Plot the results using the Visualize PDE Results Live Editor task.

1. In the Select results section of the task, select `Rt` from the drop-down list.

2. In the Specify data parameters section of the task, set Type to Temperature.

3. Select the Animate checkbox to see how the temperature changes over time.

``` ```

### Heat Flux Functions

Create a function representing a Gaussian heat flux with the center at (`x0`,`y0`) and the maximum amplitude 350000000$W/{m}^{2}$. The heat flux has a steady amplitude for 0.03 seconds, then it switches off.

```function fluxValSteady = heatFluxSteady(location,state) A = 350000000; x0 = 0; y0 = 0; sx = 0.1; % standard deviation along x sy = 0.05; % standard deviation along y fluxValSteady = A*exp(-(((location.x-x0).^2/(2*sx^2))+ ... ((location.y-y0).^2/(2*sy^2)))); % Beam switches off after 0.03 seconds if state.time > 0.03 fluxVal = 0*location.x; end end```

Create a function representing a Gaussian heat flux with the center at (`x0`,`y0`) and with the amplitude 350000000$W/{m}^{2}$ amplitude gradually decaying over time.

```function fluxValTime = heatFluxTimeDependent(location,state) t = state.time; A = 350000000; x0 = 0; y0 = 0; sx = 0.1; % standard deviation along x sy = 0.05; % standard deviation along y for index = 1:length(t) fluxValTime(:,:,index) = A*exp(-(((location.x-x0).^2/(2*sx^2)) + ... ((location.y-y0).^2/(2*sy^2)))).*exp(-t(index)); end end```