# Nonlinear finite elements/Homework 4/Solutions

## Problem 1: Axially Loaded Bar with Nonlinear Modulus

Given:

Bar of uniform cross-section under axial body load and compressive axial end load (${\displaystyle R}$ ). Nonlinear material with constitutive relation

${\displaystyle \sigma =E_{0}(1-\varepsilon )\varepsilon ,\qquad \varepsilon ={\cfrac {du}{dx}}~.}$

Finite element model uses quadratic elements with shape functions

{\displaystyle {\begin{aligned}N_{1}(x)&={\cfrac {2(x_{2}-x)(x_{3}-x)}{h^{2}}}\\N_{2}(x)&={\cfrac {4(x-x_{1})(x_{3}-x)}{h^{2}}}\\N_{3}(x)&={\cfrac {2(x-x_{1})(x-x_{2})}{h^{2}}}\end{aligned}}}

where ${\displaystyle h=x_{3}-x_{1}}$ .

### Solution

#### Part 1

Plot the shape functions for an element. Show that they add up to 1 at all points within the element.

A plot of the shape functions for an element is shown in Figure 1(a). The shape functions for a mesh are shown in Figure 1(b). If we add the shape functions we get

{\displaystyle {\begin{aligned}\sum _{i}N_{i}(x)&={\cfrac {2(x_{2}-x)(x_{3}-x)}{h^{2}}}+{\cfrac {4(x-x_{1})(x_{3}-x)}{h^{2}}}+{\cfrac {2(x-x_{1})(x-x_{2})}{h^{2}}}\\&={\cfrac {2}{h^{2}}}\left[(x_{2}-x)(x_{3}-x)+(x-x_{1})(x_{3}-x)+(x-x_{1})(x_{3}-x)+(x-x_{1})(x-x_{2})\right]\\&={\cfrac {2}{h^{2}}}\left[(x_{2}-x_{1})(x_{3}-x)+(x-x_{1})(x_{3}-x_{2})\right]\\&=\left({\cfrac {2}{h^{2}}}\right)\left({\cfrac {h}{2}}\right)(x_{3}-x_{1})=1\qquad \square \end{aligned}}}
 Figure 1 (a). Shape functions for an element. Figure 1 (b). Shape functions for a mesh.

#### Part 2

Plot the displacement, stress, strain, and the error norm when you use 10 elements to discretize the bar. Discuss and comment on your observations.

If we use an isoparametric form of the equations to simplify the integration of the stiffness matrix and load vector terms, the shape functions become

{\displaystyle {\begin{aligned}N_{1}(\xi )\,&=-{\frac {1}{2}}\left(1-\xi \right)\xi \\N_{2}(\xi )\,&=\,\left(1+\xi \right)\left(1-\xi \right)\\N_{3}(\xi )\,&={\frac {1}{2}}\left(1+\xi \right)\xi ~.\end{aligned}}}

The map between ${\displaystyle x}$  and ${\displaystyle \xi }$  is

${\displaystyle x(\xi )=x_{1}N_{1}(\xi )+x_{2}N_{2}(\xi )+x_{3}N_{3}(\xi );~\qquad x_{2}={\cfrac {x_{1}+x_{3}}{2}}~.}$

Therefore, the Jacobian of the deformation is

${\displaystyle J={\cfrac {dx}{d\xi }}={\cfrac {h}{2}}~.}$

The terms of the stiffness matrix are given by

{\displaystyle {\begin{aligned}K_{ij}&=\int _{x_{1}}^{x_{3}}AE_{0}\left(1-\sum _{k=1}^{3}u_{k}{\cfrac {dN_{k}(x)}{dx}}\right){\cfrac {dN_{i}(x)}{dx}}{\cfrac {dN_{j}(x)}{dx}}~dx\\&=\int _{-1}^{1}AE_{0}\left[1-\sum _{k=1}^{3}u_{k}\left(J^{-1}{\cfrac {dN_{k}(\xi )}{d\xi }}\right)\right]\left(J^{-1}{\cfrac {dN_{i}(\xi )}{d\xi }}\right)\left(J^{-1}{\cfrac {dN_{j}(\xi )}{d\xi }}\right)~J~d\xi ~.\end{aligned}}}

The element stiffness matrix is

${\displaystyle \mathbf {K} ={\cfrac {AE_{0}}{3h^{2}}}{\begin{bmatrix}7h+u_{3}+15u_{1}-16u_{2}&&-8(h+2u_{1}-2u_{2})&&h-u_{3}+u_{1}\\-8(h+2u_{1}-2u_{2})&&16(h-u_{3}+u_{1})&&-8(h-2u_{3}+2u_{2})\\h-u_{3}+u_{1}&&-8(h-2u_{3}+2u_{2})&&7h-15u_{3}-u_{1}+16u_{2}\end{bmatrix}}}$

The terms of the tangent stiffness matrix are given by

${\displaystyle T_{ij}=K_{ij}+\sum _{m=1}^{3}{\frac {\partial K_{im}}{\partial u_{j}}}u_{m}~.}$

Therefore, the element tangent stiffness matrix is

${\displaystyle \mathbf {T} ={\begin{bmatrix}7h+2u_{3}+30u_{1}-32u_{2}&&-8(h+4u_{1}-4u_{2})&&h-2u_{3}+2u_{1}\\-8(h+4u_{1}-4u_{2})&&16(h-2u_{3}+2u_{1})&&-8(h-4u_{3}+4u_{2})\\h-2u_{3}+2u_{1}&&-8(h-4u_{3}+4u_{2})&&7h-30u_{3}-2u_{1}+32u_{2}\end{bmatrix}}}$

The terms of the element load vector are given by

{\displaystyle {\begin{aligned}f_{i}&=\int _{x_{1}}^{x_{3}}axN_{i}(x)~dx+\left.N_{i}\left(AE{\cfrac {du}{dx}}\right)\right|_{x_{1}}^{x_{3}}\\&=\int _{-1}^{1}a\left(\sum _{k=1}^{3}x_{k}N_{k}(\xi )\right)N_{i}(\xi )~J~d\xi +\left.N_{i}(\xi )R_{i}\right|_{-1}^{1}~.\end{aligned}}}

Therefore, the element load vector is

${\displaystyle \mathbf {f} ={\cfrac {ah}{6}}{\begin{bmatrix}x_{1}\\2(x_{1}+x_{3})\\x_{3}\end{bmatrix}}+{\begin{bmatrix}-R_{1}\\0\\R_{3}\end{bmatrix}}}$

The stress-strain curve in compression is shown in Figure 2.

 Figure 2. Stress-strain curve in compression.

If we use a mesh containing ${\displaystyle 2}$  equally sized elements, we get the solutions shown in Figure 3.

 Figure 3(a). Displacement Figure 3(b). Stress
 Figure 3(a). Strain Figure 3(b). Error norm

If we use a mesh containing ${\displaystyle 10}$  equally sized elements, we get the solutions shown in Figure 4.

 Figure 4(a). Displacement Figure 4(b). Stress
 Figure 4(a). Strain Figure 4(b). Error norm

Observations:

1. The two-node and three-node elements give the same solution.
2. The solution using two elements shows that the stresses and strains are interpolated linearly within each element and are non longer constant.
3. There is no jump in the stress/strain at the nodes.
4. The solution converges quite rapidly - four iterations. This is due to the shape of the stress-strain curve in compression.

#### Part 2

Do the same for 100 elements. What are your observations on the effect of mesh refinement on the solution?

If we use a mesh containing ${\displaystyle 100}$  equally sized elements, we get the solutions shown in Figure 5.

 Figure 5(a). Displacement Figure 5(b). Stress
 Figure 5(a). Strain Figure 5(b). Error norm

Observations:

1. The solution is approximated quite accurately with 10 elements. Using 100 elements does not lead to any significant improvement in the solution.
2. The solution converges in four iterations unlike the case where a tensile load was applied.

#### Part 4

Now, increase the applied load. At what load does the solution fail to converge? Explain why the solution fails to converge at this load.

The solution converges at all loads unlike in the tension problem. If the load is increased to ${\displaystyle R=-20.1}$ , the elements invert and we get a negative Jacobian. Even though we get a solution, that solution is unphysical.

The Matlab code use to compute the results is shown below.

 function HW4Prob1_2 E0 = 10; R = -2; a = 1; W = 1; H = 1; A = 1; L = 1; tol = 1.0e-6; nElem = 10; h = L/nElem; nNode = 2*nElem+1; for i=1:nNode xnode(i) = (i-1)*h/2; end for i=1:nElem node1 = 2*(i-1)+1; elem(i,:) = [node1 node1+1 node1+2]; end u = zeros(nNode,1); condition = 1; condition_res = 1; count = 0; while (condition > tol) count = count + 1 K = zeros(nNode); T = zeros(nNode); f = zeros(nNode,1); for i=1:nElem node1 = elem(i,1); node2 = elem(i,2); node3 = elem(i,3); Ke = elementStiffness(A, E0, xnode, u, node1, node2, node3); Te = elementTangentStiffness(A, E0, xnode, u, node1, node2, node3); fe = elementLoad(a, xnode, node1, node2, node3); K(node1:node3,node1:node3) = K(node1:node3,node1:node3) + Ke; T(node1:node3,node1:node3) = T(node1:node3,node1:node3) + Te; f(node1:node3) = f(node1:node3) + fe; end f(nNode) = f(nNode) + R; u0 = 0; Kred = K(2:nNode,2:nNode); Tred = T(2:nNode,2:nNode); fred = f(2:nNode); ured = u(2:nNode); r = Kred*ured - fred; Tinv = inv(Tred); u_new = ured - Tinv*r; condition = norm(u_new - ured)/norm(u_new) u = [u0; u_new]; iter(count) = count; normres(count) = norm(r); normu(count) = condition; end figure; p01 = semilogy(iter, normres, 'ro-'); hold on; p02 = semilogy(iter, normu, 'bo-'); set([p01 p02], 'LineWidth', 2); xlabel('Iterations','FontName','palatino','FontSize',18); ylabel('Error Norm','FontName','palatino','FontSize',18); set(gca, 'LineWidth', 2, 'FontName','palatino','FontSize',18); legend([p01 p02],'|r|','|u_{n+1}-u_n|/|u_{n+1}|'); axis square; figure; p0 = plotDisp(E0, A, L, R, a); p = plot(xnode, u, 'ro-', 'LineWidth', 2); hold on; xlabel('x', 'FontName', 'palatino', 'FontSize', 18); ylabel('u(x)', 'FontName', 'palatino', 'FontSize', 18); set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18); legend([p0 p], 'Linear (Exact) ','Nonlinear (FEM)'); axis square; for i=1:nElem node1 = elem(i,1); node2 = elem(i,2); node3 = elem(i,3); [xloc(i,:), eps(i,:), sig(i,:)] = ... elementStrainStress(E0, xnode, u, node1, node2, node3); end figure; p0 = plotStress(E0, A, L, R, a); for i=1:nElem p1 = plot(xloc(i,:), sig(i,:), 'r-o','LineWidth',2); hold on; end xlabel('x', 'FontName', 'palatino', 'FontSize', 18); ylabel('\sigma(x)', 'FontName', 'palatino', 'FontSize', 18); set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18); legend([p0 p1], 'Linear (Exact) ','Nonlinear (FEM)'); axis square; figure; p0 = plotStrain(E0, A, L, R, a); for i=1:nElem p1 = plot(xloc(i,:), eps(i,:), 'r-o','LineWidth',2); hold on; end xlabel('x', 'FontName', 'palatino', 'FontSize', 18); ylabel('\epsilon(x)', 'FontName', 'palatino', 'FontSize', 18); set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18); legend([p0 p1], 'Linear (Exact) ','Nonlinear (FEM)'); axis square; function [Ke] = elementStiffness(A, E0, xnode, u, node1, node2, node3) x1 = xnode(node1); x2 = xnode(node2); x3 = xnode(node3); h = x3 - x1; u1 = u(node1); u2 = u(node2); u3 = u(node3); C1 = A*E0/(3*h^2); Ke(1,1) = 7*h + u3 + 15*u1 - 16*u2; Ke(1,2) = -8*(h + 2*u1 - 2*u2); Ke(1,3) = (h - u3 + u1); Ke(2,2) = 16*(h - u3 + u1); Ke(2,3) = -8*(h - 2*u3 + 2*u2); Ke(3,3) = (7*h - 15*u3 - u1 + 16*u2); Ke(2,1) = Ke(1,2); Ke(3,1) = Ke(1,3); Ke(3,2) = Ke(2,3); Ke = C1*Ke; function [Te] = elementTangentStiffness(A, E0, xnode, u, ... node1, node2, node3) x1 = xnode(node1); x2 = xnode(node2); x3 = xnode(node3); h = x3 - x1; u1 = u(node1); u2 = u(node2); u3 = u(node3); C1 = A*E0/(3*h^2); Te(1,1) = 7*h + 2*u3 + 30*u1 - 32*u2; Te(1,2) = -8*(h + 4*u1 - 4*u2); Te(1,3) = (h - 2*u3 + 2*u1); Te(2,2) = 16*(h - 2*u3 + 2*u1); Te(2,3) = -8*(h - 4*u3 + 4*u2); Te(3,3) = (7*h - 30*u3 - 2*u1 + 32*u2); Te(2,1) = Te(1,2); Te(3,1) = Te(1,3); Te(3,2) = Te(2,3); Te = C1*Te; function [fe] = elementLoad(a, xnode, node1, node2, node3) x1 = xnode(node1); x2 = xnode(node2); x3 = xnode(node3); h = x3 - x1; C1 = a*h/6; fe(1,1) = x1; fe(2,1) = 2*(x1+x3); fe(3,1) = x3; fe = C1*fe; function [x, eps, sig] = elementStrainStress(E0, xnode, u, ... node1, node2, node3) x1 = xnode(node1); x2 = xnode(node2); x3 = xnode(node3); h = x3 - x1; u1 = u(node1); u2 = u(node2); u3 = u(node3); ximin = -1; ximax = 1; nxi = 2; dxi = (ximax - ximin)/nxi; for i=1:nxi+1 xi = ximin + (i-1)*dxi; J = h/2; dN1_dxi = xi - 1/2; dN2_dxi = -2*xi; dN3_dxi = xi + 1/2; B = [dN1_dxi/J dN2_dxi/J dN3_dxi/J]; u = [u1; u2; u3]; eps(i) = B*u; sig(i) = E0*(1-eps(i))*eps(i); F = 1 + eps(i); if (F < 0) display('Material has inverted!'); end N1 = -0.5*(1-xi)*xi; N2 = (1+xi)*(1-xi); N3 = 0.5*(1+xi)*xi; x(i) = x1*N1 + x2*N2 + x3*N3; end function [p] = plotDisp(E, A, L, R, a) dx = 0.01; nseg = L/dx; for i=1:nseg+1 x(i) = (i-1)*dx; u(i) = 1/(6*A*E)*(-a*x(i)^3 + (6*R + 3*a*L^2)*x(i)); end p = plot(x, u, 'LineWidth', 3); hold on; xlabel('x', 'FontName', 'palatino', 'FontSize', 18); ylabel('u(x)', 'FontName', 'palatino', 'FontSize', 18); set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18); function [p] = plotStrain(E, A, L, R, a) dx = 0.01; nseg = L/dx; for i=1:nseg+1 x(i) = (i-1)*dx; eps(i) = 1/(2*A*E)*(-a*x(i)^2 + (2*R + a*L^2)); end p = plot(x, eps, 'LineWidth', 3); hold on; xlabel('x', 'FontName', 'palatino', 'FontSize', 18); ylabel('\epsilon(x)', 'FontName', 'palatino', 'FontSize', 18); set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18); function [p] = plotStress(E, A, L, R, a) dx = 0.01; nseg = L/dx; for i=1:nseg+1 x(i) = (i-1)*dx; sig(i) = 1/(2*A)*(-a*x(i)^2 + (2*R + a*L^2)); end p = plot(x, sig, 'LineWidth', 3); hold on; xlabel('x', 'FontName', 'palatino', 'FontSize', 18); ylabel('\sigma(x)', 'FontName', 'palatino', 'FontSize', 18); set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18); 

You could alternatively write your own Maple code or use ANSYS for these calculations depending on your interests.

## Problem 2: Nonlinear Thermal Conduction

Consider the model chip system shown in Figure 6.

 Figure 6. Schematic of a model chip.

Assume that heat is being generated in the silicon carbide region of the chip at the rate of 100 W/cm${\displaystyle ^{2}}$ . The cooling channels contain a fluid at a constant 250 K and (supposedly) can dissipate heat at the rate of 800 W/cm${\displaystyle 2}$ . The surroundings are at room temperature (300 K).

The thermal conductivity of the silicon carbide and the aluminum nitride are also given in Figure 6. The thermal conductivity of the epoxy resin is 1.5 W/m-K.

 Figure 7. Meshed model

#### Part 1

First run a linear simulation of the problem using room temperature thermal conductivities of the materials. This will act as a check of the model and setup and will help you determine any additional constraints that you will need. Do the cooling channels really dissipate heat at 800 W/cm${\displaystyle 2}$ ?

The following ANSYS input is used for this simulation. Due to symmetry, only half of the model will be meshed Figure 7. A constant temperature of 300K is applied around the edge Epoxy resin. Heat is being generated from areas made up by the Silicon Carbide at the rate of 100 W/cm${\displaystyle ^{2}}$ . Each line forming the cooling channels has constant temperature of 250K.

/prep7 et,1,55 mp,kxx,1,28.634 mp,kxx,2,84.873 mp,kxx,3,1.5 radius = 0.00005 ewidth = 0.0001 eheight = 0.0001 i=0 k,i+1,0,0 k,i+2,radius,0 k,i+3,ewidth,0 k,i+4,ewidth,eheight k,i+5,0,eheight k,i+6,0,radius local,11,1,0,0 k,i+7,radius,45 k,i+8,radius,-45 csys k,i+9,0,-radius k,i+10,0,-eheight k,i+11,ewidth,-eheight l,i+2,i+3,2 l,i+3,i+4,2 l,i+4,i+5,2 l,i+5,i+6,2 local,11,1,0,0 l,i+6,i+7,2 l,i+7,i+2,2 l,i+2,i+8,2 l,i+8,i+9,2 csys l,i+9,i+10,2 l,i+10,i+11,2 l,i+11,i+3,2 l,i+7,i+4,2,2 l,i+8,i+11,2 a,i+2,i+3,i+4,i+7 a,i+7,i+4,i+5,i+6 a,i+2,i+3,i+11,i+8 a,i+8,i+11,i+10,i+9 AGEN,7,all,,,7*ewidth,0,0,,1 asel,s,area,,5,28 ARSYM,X,all, , , ,1,0 asel,s,area,,29,52 AGEN,,all,,,7*ewidth*7,,,,,1 asel,all *do,i,1,37,7 blc4,ewidth*i,0,ewidth,eheight blc4,ewidth*(i+1),0,ewidth,eheight blc4,ewidth*(i+2),0,ewidth,eheight blc4,ewidth*(i+3),0,ewidth,eheight blc4,ewidth*(i+4),0,ewidth,eheight blc4,ewidth*i,-eheight,ewidth,eheight blc4,ewidth*(i+1),-eheight,ewidth,eheight blc4,ewidth*(i+2),-eheight,ewidth,eheight blc4,ewidth*(i+3),-eheight,ewidth,eheight blc4,ewidth*(i+4),-eheight,ewidth,eheight *enddo blc4,ewidth*43,0,ewidth,eheight blc4,ewidth*43,-eheight,ewidth,eheight blc4,ewidth*44,0,ewidth,eheight blc4,ewidth*44,-eheight,ewidth,eheight *do,i,1,6 blc4,(i-1)*0.0001,0.0001,0.0001,0.0001 *enddo *do,i,0,9 blc4,0.002+(0.0001*i),0.0001,0.0001,0.0001 *enddo *do,i,0,14 blc4,0.0005+(i*0.0001),0.0001,0.0001,0.0001 *enddo *do,i,0,14 blc4,0.003+(0.0001*i),0.0001,0.0001,0.0001 *enddo blc4,0.0044,-0.0001,0.0001,0.0001 blc4,0.0044,0,0.0001,0.0001 blc4,0.0044,0.0001,0.0001,0.0001 *do,i,1,51 blc4,(i-1)*ewidth,2*eheight,ewidth,0.0003 blc4,(i-1)*ewidth,-eheight,ewidth,-0.0003 *enddo *do,i,0,5 blc4,0.0045+(i*0.0001),eheight,ewidth,0.0001 blc4,0.0045+(i*0.0001),0,ewidth,0.0001 blc4,0.0045+(i*0.0001),-eheight,ewidth,0.0001 *enddo aovlap,all lesize,all,,,2 aatt,1 asel,s,loc,x,0,0.0045 asel,r,loc,y,-0.0001,0.0001 amesh,all asel,s,loc,x,0,0.0005 asel,a,loc,x,0.002,0.003 asel,r,loc,y,0.0001,0.0002 amesh,all asel,all aatt,2 asel,s,loc,x,0.0005,0.002 asel,a,loc,x,0.003,0.0045 asel,r,loc,y,0.0001,0.0002 amesh,all asel,all aatt,3 asel,s,loc,x,0,0.005 asel,r,loc,y,-0.0001,-0.0004 amesh,all asel,s,loc,x,0,0.005 asel,r,loc,y,0.0002,0.0005 amesh,all asel,s,loc,x,0.0045,0.005 asel,r,loc,y,-0.0001,0.0002 amesh,all asel,all ESEL,S,MAT,,1 /color,elem,12,all esel,all ESEL,S,MAT,,2 /color,elem,8,all esel,all ESEL,S,MAT,,3 /color,elem,4,all esel,all nsel,s,loc,x,0,0 dsym,symm,x nsel,all ASEL,S,MAT,,2 bfa,all,hgen,1000000 asel,all lsel,s,loc,x,0,0.005 lsel,r,loc,y,0.0005 DL,all, ,TEMP,300,1 lsel,s,loc,x,0,0.005 lsel,r,loc,y,-0.0004 DL,all, ,TEMP,300,1 lsel,s,loc,y,-0.0004,0.0005 lsel,r,loc,x,0.005 DL,all, ,TEMP,300,1 lsel,all LSEL,S,LENGTH,,0.3927E-04 DL,all, ,TEMP,250,1 lsel,all fini /solu solve fini /post1 plnsol,temp plnsol,tfsum fini 

Shown in Figure 8 and 9 below are the contour plots of temperature and the vector sum of heat flux. From Figure 9 we see that the cooling channels do not dissipate heat at the rate suggested in the problem statement.

 Figure 8. Temperature of the linear thermal conductivity model (K). Figure 9. Flux of the linear thermal conductivity model (W/m-K).

#### Part 2

Fit curves to the thermal conductivity versus temperature data for Silicon Carbide and Aluminum Nitride.

Below are the curve-fits for each material.

• Silicon Carbide:${\displaystyle k=2(10^{-8})T^{4}-2(10^{-5})T^{3}+0.0057T^{2}-0.7605T+87.873}$
• Aluminum Nitride:${\displaystyle k=-1(10^{-9})T^{4}+1(10^{-7})T^{3}+0.0003T^{2}-0.0941T+28.634}$

#### Part 3

Perform a nonlinear simulation using these temperature dependent properties to find the equilibrium temperature inside the silicon carbide region. State all simplifying assumptions that you have made.

The same code that was used in the previous problem is used here. However, the thermal conductivities are tabulated according to their temperature.

/prep7 et,1,55 MPTEMP,,,,,,,, MPTEMP,1,20 MPTEMP,2,60 MPTEMP,3,100 MPTEMP,4,140 MPTEMP,5,180 MPTEMP,6,220 MPTEMP,7,260 MPTEMP,8,300 MPDATA,KXX,1,,27 MPDATA,KXX,1,,24 MPDATA,KXX,1,,22 MPDATA,KXX,1,,21 MPDATA,KXX,1,,20 MPDATA,KXX,1,,20 MPDATA,KXX,1,,20 MPDATA,KXX,1,,20 MPDATA,KXX,2,,75 MPDATA,KXX,2,,59 MPDATA,KXX,2,,52 MPDATA,KXX,2,,49 MPDATA,KXX,2,,47 MPDATA,KXX,2,,45 MPDATA,KXX,2,,43 MPDATA,KXX,2,,41 mp,kxx,3,1.5 radius = 0.00005 ewidth = 0.0001 eheight = 0.0001 i=0 k,i+1,0,0 k,i+2,radius,0 k,i+3,ewidth,0 k,i+4,ewidth,eheight k,i+5,0,eheight k,i+6,0,radius local,11,1,0,0 k,i+7,radius,45 k,i+8,radius,-45 csys k,i+9,0,-radius k,i+10,0,-eheight k,i+11,ewidth,-eheight l,i+2,i+3,2 l,i+3,i+4,2 l,i+4,i+5,2 l,i+5,i+6,2 local,11,1,0,0 l,i+6,i+7,2 l,i+7,i+2,2 l,i+2,i+8,2 l,i+8,i+9,2 csys l,i+9,i+10,2 l,i+10,i+11,2 l,i+11,i+3,2 l,i+7,i+4,2,2 l,i+8,i+11,2 a,i+2,i+3,i+4,i+7 a,i+7,i+4,i+5,i+6 a,i+2,i+3,i+11,i+8 a,i+8,i+11,i+10,i+9 AGEN,7,all,,,7*ewidth,0,0,,1 asel,s,area,,5,28 ARSYM,X,all, , , ,1,0 asel,s,area,,29,52 AGEN,,all,,,7*ewidth*7,,,,,1 asel,all *do,i,1,37,7 blc4,ewidth*i,0,ewidth,eheight blc4,ewidth*(i+1),0,ewidth,eheight blc4,ewidth*(i+2),0,ewidth,eheight blc4,ewidth*(i+3),0,ewidth,eheight blc4,ewidth*(i+4),0,ewidth,eheight blc4,ewidth*i,-eheight,ewidth,eheight blc4,ewidth*(i+1),-eheight,ewidth,eheight blc4,ewidth*(i+2),-eheight,ewidth,eheight blc4,ewidth*(i+3),-eheight,ewidth,eheight blc4,ewidth*(i+4),-eheight,ewidth,eheight *enddo blc4,ewidth*43,0,ewidth,eheight blc4,ewidth*43,-eheight,ewidth,eheight blc4,ewidth*44,0,ewidth,eheight blc4,ewidth*44,-eheight,ewidth,eheight *do,i,1,6 blc4,(i-1)*0.0001,0.0001,0.0001,0.0001 *enddo *do,i,0,9 blc4,0.002+(0.0001*i),0.0001,0.0001,0.0001 *enddo *do,i,0,14 blc4,0.0005+(i*0.0001),0.0001,0.0001,0.0001 *enddo *do,i,0,14 blc4,0.003+(0.0001*i),0.0001,0.0001,0.0001 *enddo blc4,0.0044,-0.0001,0.0001,0.0001 blc4,0.0044,0,0.0001,0.0001 blc4,0.0044,0.0001,0.0001,0.0001 *do,i,1,51 blc4,(i-1)*ewidth,2*eheight,ewidth,0.0003 blc4,(i-1)*ewidth,-eheight,ewidth,-0.0003 *enddo *do,i,0,5 blc4,0.0045+(i*0.0001),eheight,ewidth,0.0001 blc4,0.0045+(i*0.0001),0,ewidth,0.0001 blc4,0.0045+(i*0.0001),-eheight,ewidth,0.0001 *enddo aovlap,all lesize,all,,,2 aatt,1 asel,s,loc,x,0,0.0045 asel,r,loc,y,-0.0001,0.0001 amesh,all asel,s,loc,x,0,0.0005 asel,a,loc,x,0.002,0.003 asel,r,loc,y,0.0001,0.0002 amesh,all asel,all aatt,2 asel,s,loc,x,0.0005,0.002 asel,a,loc,x,0.003,0.0045 asel,r,loc,y,0.0001,0.0002 amesh,all asel,all aatt,3 asel,s,loc,x,0,0.005 asel,r,loc,y,-0.0001,-0.0004 amesh,all asel,s,loc,x,0,0.005 asel,r,loc,y,0.0002,0.0005 amesh,all asel,s,loc,x,0.0045,0.005 asel,r,loc,y,-0.0001,0.0002 amesh,all asel,all ESEL,S,MAT,,1 /color,elem,12,all esel,all ESEL,S,MAT,,2 /color,elem,8,all esel,all ESEL,S,MAT,,3 /color,elem,4,all esel,all KBC,1 nsel,s,loc,x,0,0 dsym,symm,x nsel,all ASEL,S,MAT,,2 bfa,all,hgen,1000000 asel,all lsel,s,loc,x,0,0.005 lsel,r,loc,y,0.0005 DL,all, ,TEMP,300,1 lsel,s,loc,x,0,0.005 lsel,r,loc,y,-0.0004 DL,all, ,TEMP,300,1 lsel,s,loc,y,-0.0004,0.0005 lsel,r,loc,x,0.005 DL,all, ,TEMP,300,1 lsel,all LSEL,S,LENGTH,,0.3927E-04 DL,all, ,TEMP,250,1 lsel,all fini /solu solve fini /post1 plnsol,temp plnsol,tfsum fini 

Shown in Figure 10 and 11 below are the contour plots of temperature and the vector sum of heat flux.

 Figure 10. Temperature of the nonlinear thermal conductivity model (K). Figure 11. Flux of the nonlinear thermal conductivity model (W/m-K).

There are many approaches to create the mesh for this problem. You will have to find the way the work best for you. These two codes listed are just two ways that the meshes can be manipulated.