Anisotropic Heat Transfer Through Woven Carbon Fibers

Composite materials such as carbon-fiber-reinforced polymer (CFRP) have outstanding properties. It is a lightweight material with high stiffness and high temperature tolerance and therefore used in aerospace industry, civil engineering, and also for high-end sports goods.

Carbon crystals form flat ribbons that are bundled together and woven in different structures as required by the application area. The bundles have anisotropic material properties. For thermal properties this means that the thermal conductivity along the fiber axis is much higher than perpendicular to it.

Figure 1: Model geometry: fibers embedded in an epoxy matrix (hidden) with infinite element domain.

In COMSOL Multiphysics, implementing anisotropic properties is straightforward in the global coordinate system, which is set by default. However, in the present case, an anisotropic thermal conductivity needs to be defined along woven fibers where the global coordinate system is inapplicable. The curvilinear coordinate system provides the possibility to create a coordinate system following the curves of a geometry in which anisotropic material properties or anisotropic physics can be defined.

This tutorial shows how to use the curvilinear coordinates interface and how to apply it to define anisotropic thermal conductivity.

Model Definition

The model represents a cutout of a carbon-fiber-reinforced polymer. The geometry used in this case is shown in Figure 1. The fiber bundles have circular cross-sections and are embedded in a matrix made of epoxy.

The “infinite domains” truncate the geometry to model a few fibers only. With the Heat Transfer Module you can assign them as Infinite Element Domains and thus suppress boundary effects. Without the Heat Transfer Module, the boundary conditions at the outer side affects the solution, in this case the maximum temperature. Increase the number of fibers to reduce these effects.

Material Properties

The material properties are summarized in Table 1.

Table 1: Material Properties.
Material Property	Epoxy	Carbon (Core)	Carbon (infinite domain)
Thermal conductivity	0.2 W/(m·K)	{60,4,4} W/(m·K)	60 W/(m·K)
Density	1200 kg/m³	1500 kg/m³	1500 kg/m³
Heat capacity at constant pressure	1000 J/(kg·K)	1000 J/(kg·K)	1000 J/(kg·K)

Note the syntax of the thermal conductivity for carbon (core). In the general case of an anisotropic thermal conductivity, it is a second order tensor. In the present case, the tensor is diagonal.

Note that the conductivity is higher in the fibers direction and lower in perpendicular direction. The coordinate system used for k must then provide an x-component following the shape of the fibers. The Curvilinear Coordinates interface provides appropriate tools to create such a base vector system.

Curvilinear Coordinates

Three predefined methods and a user-defined method are available to set up a curvilinear coordinate system. Further details can be found in the section The Curvilinear Coordinates Interface in the COMSOL Multiphysics Reference Manual. Here you use the diffusion method which solves Laplace’s equation resulting in a scalar potential. It is the same as solving the stationary heat transfer equation with temperature boundary conditions resulting in a temperature gradient and forming the first base vector of the new coordinate system. The second base vector is specified manually and the cross product of both forms the third base vector.

Figure 2 shows the base vector system for a single fiber.

Figure 2: Curvilinear coordinate system from diffusion method.

Alternatively the Flow Method is available, which results in a vector potential. This is equivalent to solving Stokes flow (also known as creeping flow) where the obtained velocity field forms the first base vector. The third option is to choose the Elasticity Method for solving an eigenvalue problem.

If the option is selected, the new curvilinear system is available as input for the list and thus providing new (x, y, and z) coordinates.

Boundary Conditions

For the curvilinear coordinates interface, the inlet and outlet boundaries define the direction of the first base vector. The heat transfer analogy consists of setting a high temperature at the inlet and a low temperature at the outlet. All other boundaries are thermally insulated walls.

For the heat transfer interface, a constant temperature boundary condition is set at the outermost walls. A boundary heat source described with a Gaussian pulse is applied to the center of the geometry and a convective cooling boundary condition is applied on both sides.

Infinite Elements

To truncate the geometry the Infinite Element Domain feature can be used. Boundary conditions applied to these elements can be imagined as boundary conditions at an infinite distance of the modeling domain. So it does not affect the solution of this particular problem. This works by scaling the width of the domain to be much larger than the original geometry.

Results and Discussion

From the interface a new coordinate system is obtained as shown in Figure 2. The temperature distribution on the surface shows a high temperature at the center where the maximum of the Gaussian function is located and decreases with increasing distance from the center. The temperature drop of 293 K, as specified in the boundary conditions, occurs mainly in the Infinite Element Domains.

Figure 3: Temperature distribution on the surface.

Figure 4 shows clearly that the heat spreads preferentially along the fiber axis.

Figure 4: Temperature at the center plane and fiber structure (gray).

Notes About the COMSOL Implementation

This tutorial demonstrates how to use the interface for defining anisotropic thermal conductivity. Hence, the instructions focus on this part and start with loading the file carbon_fiber_geom.mph. The steps needed to create this file are quite complex. This document does not go into details but provides a short summary instead.

The geometry sequence calls geometry subsequences depending on the global parameter q. The subsequences define different cross-sections and can be found under the node. Call the elliptical cross-section with q = 1 and the rectangular cross-section with q = 2.

All subsequent geometry features are based on these subsequences. Inside some of the features, a selection of geometric entities is created automatically by selecting the option. Instead of selecting objects manually, these selections are used as input in the following geometry node. This approach ensures that all geometry operations adapt and produce the desired geometry automatically, even if a geometry parameter changes.

Selections are also used on the finalized geometry to ensure that physical properties are assigned to the intended entities. These selections are defined under . They are used to automatically set up boundary and domain conditions as well as the mesh. The resulting model is consistent for any choice of parameters. The extra time needed to set up this kind of geometry sequence and to define selections is regained through an accelerated physics modeling and meshing process.

Heat_Transfer_Module/Tutorials,_Conduction/carbon_fibers_infinite_elements

Modeling Instructions

Root

Start with loading the model file that contains the geometry and selections used throughout the modeling process.

Application Libraries

From the menu, choose .

In the window, select in the tree.

Click

Root

Click the

button in the toolbar.

Component 1 (comp1)

Add the interface for the fibers.

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Curvilinear Coordinates (cc)

In the window for , locate the section.

From the list, choose .

Locate the section. Select the check box.

According to the section, the second basis vector is specified manually. The y direction feels natural.

Coordinate System Settings 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Diffusion Method 1

In the toolbar, click

and choose .

is the default boundary condition where the normal component of the vector field is zero. The direction of the first basis vector is specified with inlet and outlet boundary conditions.

In the window for , locate the section.

From the list, choose .

Inlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Diffusion Method 1

In the window, click .

Outlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Mesh 1

Now, build a suitable mesh manually. Start with meshing the fibers.

Free Triangular 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Distribution 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type 2.

Click

Swept 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Distribution 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type 8.

Click

Free Tetrahedral 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Domain 12 only.

Swept 2

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Distribution 1

Right-click and choose .

In the window for , locate the section.

In the text field, type 3.

Mesh the remaining parts with a free swept mesh.

Click

The final mesh looks like this.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Add a stationary study to compute the new coordinate system with the .

Find the subsection. In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Study 1

Step 1: Stationary

In the toolbar, click

Results

Vector Field (cc)

The default plots show the coordinate system with streamlines for the vector field. To create the plot shown in Figure 2, add a selection to the dataset. The plot group will then use this subset of the whole geometry only.

Study 1/Solution 1 (sol1)

In the window, expand the node, then click .

Selection

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click

In the dialog box, type 20, 34, 54, 68, 88, 102, 122, 136 in the text field.

Click .

Coordinate System Volume 1

In the window, expand the node, then click .

In the window for , locate the section.

Find the subsection. In the text field, type 16.

Find the subsection. In the text field, type 2.

In the toolbar, click

Click the

button in the toolbar.

Now, add the interface to the component.

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Find the subsection. In the table, clear the check box for .

Click in the window toolbar.

In the toolbar, click

to close the window.

Heat Transfer in Solids (ht)

Click the

button in the toolbar.

Definitions

Infinite Element Domain 1 (ie1)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

In order to apply a heat source right in the center of the model, the feature is used. It computes the center of mass, which is automatically the center of the geometry.

Mass Properties 1 (mass1)

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Define a heat source via a local variable, which is defined on the boundary only. Use the mass properties variable for the center of mass to apply the source term exactly in the center.

Variables 1

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Q_in	1e5[W/m^2]exp(-5e6[1/m^2]((x-mass1.CMX)^2+(z-mass1.CMZ)^2))	W/m²	Boundary heat source

In the next section, you define the materials according to the Material Properties section.

Materials

Epoxy

In the toolbar, click

In the window for , type Epoxy in the text field.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	0.2	W/(m·K)	Basic
Density	rho	1200	kg/m³	Basic
Heat capacity at constant pressure	Cp	1000	J/(kg·K)	Basic

Carbon

In the toolbar, click

In the window for , type Carbon in the text field.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Thermal conductivity	{k11, k22, k33} ; kij = 0	{60, 4, 4}	W/(m·K)	Basic
Density	rho	1500	kg/m³	Basic
Heat capacity at constant pressure	Cp	1000	J/(kg·K)	Basic

Carbon (Infinite Element Domain)

In the toolbar, click

In the window for , type Carbon (Infinite Element Domain) in the text field.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	60	W/(m·K)	Basic
Density	rho	1500	kg/m³	Basic
Heat capacity at constant pressure	Cp	1000	J/(kg·K)	Basic

Add a second node for the fibers and choose the curvilinear system as reference system. This way the thermal conductivity is high along the fiber axis and low perpendicular to it.

Heat Transfer in Solids (ht)

In the window, under click .

Solid 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Set up the boundary conditions, consisting in a heat source, a convective heat flux accounting for cooling, and a fixed temperature at the very outer boundaries of the infinite domain.

Boundary Heat Source 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Qb text field, type Q_in.

Heat Flux 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the h text field, type 10.

Temperature 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

The first study was used to compute the curvilinear system. Add a second study to solve for the heat transfer only. Refer to in the section by selecting the solution as input for the variables not solved in this second study. This way the new coordinate system, which is initially unknown by , can be used for the heat transfer computation.

Add Study

In the toolbar, click