conceptual understanding of quadratic regression It's clear to me how to interpret the coefficients of a quadratic regression:
data <- data.frame(hours=c(6, 9, 12, 14, 30, 35, 40, 47, 51, 55, 60),
                   happiness=c(14, 28, 50, 70, 89, 94, 90, 75, 59, 44, 27))

data$hours2 <- data$hours^2

quadraticModel <- lm(happiness ~ hours + hours2, data=data)

summary(quadraticModel)



           Estimate Std. Error t value Pr(>|t|)    
(Intercept) -18.25364    6.18507  -2.951   0.0184 *  
hours         6.74436    0.48551  13.891 6.98e-07 ***
hours2       -0.10120    0.00746 -13.565 8.38e-07 ***

happinessPredict <- predict(quadraticModel,list(hours=hourValues, hours2=hourValues^2))

plot(data$hours, data$happiness, pch=16)
lines(hourValues, happinessPredict, col='blue')

however, what isn't clear is why this works. Both hours and hours2 increase ever more positively. How does squaring hours and add it to the model allow to capture the quadratic trend?
Is there anyone who could provide me with a non-mathematical explanation for this?

The individual associations of your hours and hours2 with happiness are extremely weak in your example, and nothing completely "non-mathematical" can explain this. Maybe the following plot can help illustrate how multiple regression allows the predictor hours2 to improve on predictions based solely on hours.

The values are circles. The dashed black line shows the linear association of happiness with hours alone. Not very good, not even "statistically significant" (p = 0.53 for the hours coefficient).
The solid black line shows the full model. You might think of this as starting with a linear extrapolation of the values near hours = 0, shown in the blue line. You might then think of the (squared) hours2 term as providing a non-linear correction to that extrapolation. Subtract the red curve from the blue line and you get the full model.
Code in R:
plot(happiness~hours,data,bty="n",xlim=c(0,60),ylim=c(0,300))
abline(lm(happiness~hours,data),lty=2)
abline(-18.2536,6.7444,col="blue") # "extrapolation" from 'hours' near 0
curve(.1012*x^2,from=0,to=60,add=TRUE,col="red") # non-linear "correction"
curve(-18.2536+6.7444*x-0.1012*x^2,from=0,to=60,add=TRUE) # full model
legend("topleft",bty="n",
        legend="black dashed, linear 'hours' alone
               \nblack solid, full model
               \nblue, 'hours' component, full model
               \nred, negative of 'hours2' component, full model")

