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 ( $R$ ). Nonlinear material with constitutive relation

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

Finite element model uses quadratic elements with shape functions

{\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 $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

{\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

{\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 $x$ and $\xi$ is

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

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

The terms of the stiffness matrix are given by

{\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

\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

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

Therefore, the element tangent stiffness matrix is

\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

{\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

\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 $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 $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:

The two-node and three-node elements give the same solution.
The solution using two elements shows that the stresses and strains are interpolated linearly within each element and are non longer constant.
There is no jump in the stress/strain at the nodes.
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 $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:

The solution is approximated quite accurately with 10 elements. Using 100 elements does not lead to any significant improvement in the solution.
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 $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 $^{2}$ . The cooling channels contain a fluid at a constant 250 K and (supposedly) can dissipate heat at the rate of 800 W/cm $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 $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 $^{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: $k=2(10^{-8})T^{4}-2(10^{-5})T^{3}+0.0057T^{2}-0.7605T+87.873$
Aluminum Nitride: $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.

Nonlinear finite elements/Homework 4