Skip to content

Add W alloy as material #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
N-Schwarz wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Add W alloy as material #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1f3f815
Add W alloy and improved namelist
Sep 18, 2025
5706bfe
add Makefile
Sep 18, 2025
c27a2e1
Zgeo for angle
Sep 18, 2025
392345f
Make T_init an input parameter
Sep 18, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 68 additions & 0 deletions tasks/loads_3d_walls/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
# Intel Fortran compiler
FC = ifort

# Base flags
RELEASEFLAGS = -O2 -fpp
DEBUGFLAGS = -O0 -g -check all -warn all -traceback -fpp

# OpenMP flag
OMPFLAG = -qopenmp

# Choose build type (default = release)
BUILD ?= release
ifeq ($(BUILD),debug)
FFLAGS = $(DEBUGFLAGS)
else
FFLAGS = $(RELEASEFLAGS)
endif

ifeq ($(WALL),W_alloy)
DEFINED = -DW_alloy
endif

ifeq ($(WALL),Be)
DEFINED = -DBE_WALL
endif

# Module source
MODSRC = mod_utils.f90
OBJMOD = mod_utils.o

# Executables and their sources
EXE1 = get_data_files
EXE2 = analyze_wall_T
EXE3 = calculate_T

OBJ1 = $(OBJMOD) get_data_files.o
OBJ2 = $(OBJMOD) analyze_wall_T.o
OBJ3 = $(OBJMOD) calculate_T.o

# Default target
all: $(EXE1) $(EXE2) $(EXE3)

# Dependencies: programs need the module
get_data_files.o: $(OBJMOD)
analyze_wall_T.o: $(OBJMOD)
calculate_T.o: $(OBJMOD)

# Module rule
$(OBJMOD): $(MODSRC)
$(FC) $(FFLAGS) $(DEFINED) -c $< -o $@

# Generic rule for objects
%.o: %.f90
$(FC) $(FFLAGS) $(DEFINED) -c $< -o $@

# Link executables
$(EXE1): $(OBJ1)
$(FC) $(FFLAGS) -o $@ $(OBJ1)

$(EXE2): $(OBJ2)
$(FC) $(FFLAGS) $(DEFINED) -o $@ $(OBJ2)

$(EXE3): $(OBJ3)
$(FC) $(FFLAGS) $(DEFINED) $(OMPFLAG) -o $@ $(OBJ3)

# Clean
clean:
rm -f *.o *.mod $(EXE1) $(EXE2) $(EXE3)
16 changes: 11 additions & 5 deletions tasks/loads_3d_walls/analyze_wall_T.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,18 @@ program analyze_wall_T

implicit none

#ifdef BE_WALL
real*8, parameter :: T_melt = 1556 ! T (K)
#elif defined(W_alloy)
real*8, parameter :: T_melt = 1800 ! T(K), approximate for W97Ni2Fe1
#else
real*8, parameter :: T_melt = 3695 ! T(K)
#endif

! Declarations
integer :: n_nodes, n_tri, i, i1,i2, i3, ierr
integer :: i_begin, i_end, i_step, i_jump_steps
real(8), parameter :: R_geo=2.8, T_melt=1556
real(8) :: R_geo=2.8d0, Z_geo=0.d0 ! R_geo, Z_geo for theta calculation
character(len=64) :: file_name, wall_f_name
real(8), allocatable :: nodes_xyz(:,:)
integer, allocatable :: indices(:,:)
Expand All @@ -18,7 +26,7 @@ program analyze_wall_T
real(8), dimension(3):: v21, v31, norm_tria
real(8), allocatable :: T_max(:), melt_duration(:), E_dep(:), tri_areas(:)

call read_namelist(i_begin, i_end, i_jump_steps, wall_f_name)
call read_namelist(i_begin, i_end, i_jump_steps, wall_f_name, R_geo=R_geo, Z_geo=Z_geo)

do i_step = i_begin, i_end, i_jump_steps

Expand Down Expand Up @@ -81,9 +89,7 @@ program analyze_wall_T
tri_R = sqrt(tri_x**2 +tri_y**2)

tri_phi = atan2(-tri_y,tri_x)
tri_theta = atan2(-(tri_Z-0.0d0),(tri_R-R_geo))

tri_theta = atan2(-(tri_Z-1.49),(tri_R-2.6)) ! UDP arc center
tri_theta = atan2(-(tri_Z-Z_geo),(tri_R-R_geo))

write(20, '(6ES14.6)') tri_phi, tri_theta, T_max(i), melt_duration(i), E_dep(i), tri_areas(i)

Expand Down
47 changes: 40 additions & 7 deletions tasks/loads_3d_walls/calculate_T.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,15 @@ module heat_diffusion
implicit none

#ifdef BE_WALL
real*8, parameter :: mat_density = 1.750d3 ! kg/m3
real*8, parameter :: mat_density = 1.750d3 ! kg/m3
real*8, parameter :: alpha_max = 6.75d-5
#elif defined(W_alloy)
!W97Ni2Fe1
real*8, parameter :: mat_density =18.5d3 !kg/m3
real*8, parameter :: alpha_max = 3.5d-5
#else
real*8, parameter :: mat_density = 19.254d3 ! kg/m3
real*8, parameter :: mat_density = 19.254d3 ! kg/m3
real*8, parameter :: alpha_max = 6.75d-5
#endif

contains
Expand Down Expand Up @@ -64,6 +70,12 @@ pure function spec_heat_capacity(T) result(cp)
cp = -1.338948d-11*TC**4 + 1.310600d-6*TC**3 - 3.143247d-3*TC**2 + 3.344647d0*TC + 1.741434d3
if (TC > 1283) cp = 3590
if (TC < 20) cp = 1807
#elif defined(W_alloy)
cp = 4.05260522e-09*T**3 -1.65961644e-05*T**2 + 4.79594887e-02*T + 1.25588588e+02
if (T>1800) then
TC = 1800.d0
cp = 4.05260522e-09*TC**3 -1.65961644e-05*TC**2 + 4.79594887e-02*TC + 1.25588588e+02
end if
#else
! DOI: 10.1007/978-94-007-7587-9_3
cp = 135.76 + 9.1159d-3*T + 2.3134d-9*T**3 - 6.5233d5 /T**2
Expand All @@ -89,6 +101,9 @@ pure function dspec_heat_capacity_dT(T) result(cp)
cp = -5.355792d-11*TC**3 + 3.9318d-6*TC**2 - 6.286494d-3*TC + 3.344647d0
if (TC > 1283) cp = 0
if (TC < 20) cp = 0
#elif defined(W_alloy)
cp = 1.2157815660000002e-08*T**2 -3.31923288e-05*T + 4.79594887e-02
if (T>1800) cp = 0
#else
! DOI: 10.1007/978-94-007-7587-9_3
cp = 9.1159d-3 + 6.9402d-9*T**2 + 13.0466d5 /T**3
Expand All @@ -113,7 +128,13 @@ pure function heat_conductivity(T) result(K)
K = 2.472945d-10*TC**4 - 7.354535d-7*TC**3 + 7.959136d-4*TC**2 - 4.470410d-1*TC + 2.073906d2
K = min(K,200.d0) !20 C
K = max(K,60.d0) !1283 C
if (TC > 1283) K = 84.7
if (TC > 1283) K = 84.7
#elif defined(W_alloy)
K = 4.48276523e-12*T**4 -1.68635680e-08*T**3 + 1.23066061e-05*T**2 + 1.91911470e-02*T + 7.08049586e+01
if (T>1800) then
TC = 1800.d0
K = 4.48276523e-12*TC**4 -1.68635680e-08*TC**3 + 1.23066061e-05*TC**2 + 1.91911470e-02*TC + 7.08049586e+01
end if
#else
! DOI: 10.1007/978-94-00 7-7587-9_3
K = 108.34 - 1.052d-2*T + 23419.9/T
Expand All @@ -139,6 +160,9 @@ pure function dheat_conductivity_dT(T) result(K)
K = 9.89178d-10*TC**3 - 22.063605d-7*TC**2 + 15.918272d-4*TC - 4.470410d-1
if (TC > 1283) K = 0
if (TC < 20) K = 0
#elif defined(W_alloy)
K = 1.793106092e-11*T**3 -5.0590704e-08*T**2 + 2.46132122e-05*T + 1.91911470e-02
if (T > 1900) K = 0
#else
! DOI: 10.1007/978-94-00 7-7587-9_3
K = -1.052d-2 - 23419.9/T**2
Expand Down Expand Up @@ -175,12 +199,21 @@ program calculate_T
integer, parameter :: nx=120
integer :: nt, i_time, ierr, i
real*8, parameter :: L=0.012d0
real*8, parameter :: stab_param = 0.01d0, T_init=473.d0, alpha_max = 6.75d-5
real(8) :: dx, dt, t_interval, time_now, time_before
real*8, parameter :: stab_param = 0.01d0
real(8) :: dx, dt, t_interval, time_now, time_before, frad=0.d0, T_init=473.d0
real(8), allocatable :: T_curr(:,:), q_perp(:)

call read_namelist(i_begin, i_end, i_jump_steps, wall_f_name)

#ifdef BE_WALL
write(*,*) 'Using Be wall'
#elif defined(W_alloy)
write(*,*) 'Using W alloy wall'
#else
write(*,*) 'Using W wall'
#endif

call read_namelist(i_begin, i_end, i_jump_steps, wall_f_name, frad=frad, T_init=T_init)

dx = L / real(nx - 1)
write(*,*) ' heat diffusivity = ', alpha_max, ' m2/s'
write(*,*) ' Radial grid width for T = ', dx
Expand Down Expand Up @@ -258,7 +291,7 @@ program calculate_T
! --- Get q_perp for wetted triangles
do i_tri=1, n_tri
if (ind_qperp(i_tri)>0) then
q_perp(ind_qperp(i_tri)) = q_heat_perp_3d(i_tri)*(1.d0-0.22d0)
q_perp(ind_qperp(i_tri)) = q_heat_perp_3d(i_tri)*(1.d0-frad)
endif
enddo

Expand Down
10 changes: 7 additions & 3 deletions tasks/loads_3d_walls/mod_utils.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -295,16 +295,20 @@ end subroutine read_data



subroutine read_namelist(i_begin, i_end, i_jump_steps, file_name)
subroutine read_namelist(i_begin, i_end, i_jump_steps, file_name, frad, R_geo, Z_geo, T_init)
!! Reads Namelist from given file.
integer, intent(inout) :: i_begin, i_end, i_jump_steps
integer :: fu, rc
logical :: file_exists
real*8, optional,intent(inout) :: frad !> radiation fraction (if simulation has no radiation,
!> but a radiation fraction is assumed in reality)
real*8, optional,intent(inout) :: R_geo, Z_geo !> for theta calculation
real*8, optional,intent(inout) :: T_init !> initial T for temperature calculation (K)
character(len=64),intent(inout) :: file_name
character(len=64) :: wall_f_name

! Namelist definition.
namelist /restart_inputs/ i_begin, i_end, i_jump_steps, wall_f_name
namelist /restart_inputs/ i_begin, i_end, i_jump_steps, wall_f_name, frad, R_geo, Z_geo, T_init

! Check whether file exists.
inquire (file='heat_load.nml', exist=file_exists)
Expand All @@ -327,4 +331,4 @@ end subroutine read_namelist



end module mod_utils
end module mod_utils